Моделирование работы артезианской скважины для ОАО «Труновское»

Министерство образования и науки Российской Федерации

ГОУ Высшего профессионального образования

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

ДИПЛОМНАЯ РАБОТА № ______

Студента Котвицкого Сергея Григорьевича

Специальности 230401 «Прикладная математика»

Защищен

Тема Моделирование работы артезианской скважины для ОАО «Труновское»

Приказ о закреплении темы от « 27 » января 2011 г. №235/c

Чертежи листов

Пояснительная записка листов

Подпись лица, принявшего документы на кафедру

Министерство образования и науки Российской Федерации

ГОУ Высшего профессионального образования

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

Факультет Информационных технологий и телекоммуникаций

Кафедра Прикладная математика и компьютерные технологии

УТВЕРЖДАЮ

Зав. кафедрой

подпись, инициалы, фамилия

« « 2011 г.

ЗАДАНИЕ НА ДИПЛОМНУЮ РАБОТУ

Студент Котвицкий Сергей Григорьевич группа ПМ-061

1. Тема: Моделирование работы артезианской скважины для ОАО «Труновское»

утверждена приказом по СевКавГТУ № 235/с от «27»_января_2011 года.

2. Срок представления работы к защите « 20 « июня 2011г.

3. Исходные данные для научного исследования 1. Отчёт по производственной практике, 2. Материалы из литературы на тему дипломной работы

4. Содержание дипломной работы:

4.1. Задачи с подвижными границами: Общая информация о подземных водах, Задача Стефана о фазовом переходе, Задачи с подвижными границами, решаемые разложением в ряд.

4.2. Математическая модель работы артезианской скважины: Плоскорадиальный фильтрационный поток воды к артезианской скважине (постановка задачи), Приближенное аналитическое решение задачи о нестационарной плоскорадиальной фильтрации воды, Первая фаза истощения водоносного пласта, Вторая фаза истощения водоносного пласта, Численное решение задачи с помощью пакета Maple. Сравнение с результатами приближённого аналитического решения

4.3. Заключение_______________________________________________

4.4. Список использованной литературы:

5. Приложение: Приложение 1, Приложение 2

Руководитель работы Гоголева С. А. 

Консультанты по разделам:

«Информационные технологии и программирование» Ляхов В. Ф.

Задание к исполнению принял » 27 » января 2011 г.

подпись

ВВЕДЕНИЕ

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

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

1. ЗАДАЧИ С ПОДВИЖНЫМИ ГРАНИЦАМИ

1.1 Общая информация о подземных водах

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

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

По геолого-гидрогеологическим условиям могут быть выделены следующие основные типы геологических структур и образований, а также связанные с ними подземные воды [1]:

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

На долю речных долин и артезианских бассейнов платформ приходится более 60% всех разведанных и эксплуатируемых участков подземных вод.

В зависимости от условий залегания и гидродинамических особенностей подземные воды делят на верховодку, грунтовые и артезианские. В северных и северо-восточных районах России, находящихся в пределах зоны многолетнемерзлых пород, подземные воды делят на три типа: 1) надмерзлотные, залегающие над толщей многолетней мерзлоты, служащей для них водоупором; 2) межмерзлотные, заключенные внутри толщи многолетней мерзлоты; 3) подмерзлотные, находящиеся ниже толщи многолетней мерзлоты. Для водоснабжения используют, в основном подмерзлотные и межмерзлотные воды, к которым относятся подземные воды, проходящие в трещиноватых осадочных и изверженных породах Алданского района и аллювиальных отложениях речных долин в северных районах европейской части страны, в Сибири, а также в южных районах Восточной Сибири и Дальнего Востока [2].

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

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

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

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

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

Классификация источников подземных вод по величине их дебита приведена в табл. 1.1.

Таблица 1.1 — Классификация источников подземных вод по величине дебита

Класс	Название по дебиту	Дебит
		л/с	м/сут
I	Гигантские	10000	864000
II	Исполинские	10000-1000	864000-86400
III	Очень большие	1000-100	86400-8640
IV	Большие	100-10	8640-864
V	Значительные	10-1	864-86,4
VI	Малые	1-0,1	86,4-8,64
VII	Незначительные	0,1-0,01	8,64-0,864
VIII	Весьма незначительные	0,01	<0,864

Артезианская скважина — это буровая скважина, которая пробурена для эксплуатации подземных вод.

Артезианские водоносные горизонты залегают между двумя водоупорными слоями и надежно защищены от поверхностного загрязнения. В отличие от грунтовых вод они часто имеют отдаленную область питания — за несколько километров и даже за десятки и сотни километров. При вскрытии скважиной уровень артезианской воды всегда устанавливается значительно выше водоупорной кровли водоносного горизонта, а иногда артезианская вода сама изливается из скважины (фонтанирует). На тех участках, где артезианские воды получают питание, они приобретают характер или грунтовых со свободной поверхностью, или межпластовых грунтовых вод. Подземные воды всех перечисленных видов могут циркулировать в пустотах рыхлых зернистых или в трещинах скальных пород. В последнем случае подземные воды, относящиеся к любому из перечисленных видов, получают дополнительное название трещинных [3].

Работа артезианской скважины в общем виде описывается уравнением фильтрации:

(1.1.1)

Уравнение (1.1.1) есть дифференциальное уравнение движения жидкости в пористой среде. Упомянутое уравнение имеет вид «уравнения теплопроводности», интегрирование которого при различных начальных и граничных условиях рассматривается в каждом курсе математической физики[4]. В данном случае нас интересует уравнение (1.1.1) в совокупности с граничными условиями:

(1.1.2)

— (1.1.3)

То есть нас интересует дифференциальное уравнение параболического типа (1.1.1) с граничными условиями (1.1.2) и (1.1.3). Данная задача называется — задачей с подвижными границами. Ниже рассмотрим основные задачи такого типа с их решениями.

1.2 Задача Стефана о фазовом переходе

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

Пусть жидкая среда занимает полупространство и пусть при температура всех слоёв жидкости одинакова и равна , где — температура отвердевания жидкости. Без ограничения общности мы будем считать .

С момента времени на границе поддерживается постоянная температура ниже температуры кристаллизации . В этом случае при вблизи граничной поверхности возникает слой твердой фазы, толщина которого с течением времени увеличивается (рис. 1). Фронт кристаллизации в любой момент времени отделяет твёрдую фазу от жидкой, двигаясь с некоторой скоростью в направлении жидкой фазы. По постановке задачи [5].

Рисунок 1 — Процесс крисстализации

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

Ясно выделяя движущийся фронт кристаллизации, обозначим индексом “1” величины, относящиеся к твёрдой фазы, а индексом “2” — к жидкой фазе. Тогда, считая, что свойства среды при фазовом переходе изменяются скачком, запишем уравнение теплопроводности для двух фаз:

, , ; (1.2.1)

, , ; (1.2.2)

где и — коэффициенты температуропроводности твёрдой и жидкой фаз соответственно.

Учитывая, что в начальный момент времени существует только жидкая фаза, начальное условие для задачи запишем в виде

, (1.2.3)

Краевые условия задачи сформулируем следующим образом:

а) на границах области

при ; (1.2.4)

при ;

б) на фронте фазового перехода

; (1.2.5)

, (1.2.6)

где — скрытая теплота кристаллизации, отнесённая к единице массы твёрдой фазы.

С помощью автомодельной переменной (преобразование Больцмана) приведём уравнение (1.2.1) и (1.2.2) к обыкновенным дифференциальным уравнениям для функций и :

, (1.2.7)

Полагая , запишем уравнение (1.2.7) в виде

(1.2.8)

Интегрируя выражение (1.2.8), получаем

, (1.2.9)

Интегрируя (1.2.9) ещё один раз находим общее решение уравнений (1.2.7) для и :

,

где ; , .

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

Для малых имеет место разложение

,

а при больших справедлива асимптотическая формула

Возвращаясь к переменным и , запишем найденные решения уравнений (1.2.1) и (1.2.2) в виде

, , ; (1.2.10)

, , ;

Выполняя теперь граничные условия (1.2.4) находим

, (1.2.11)

При этом замечаем, что начальное условие (1.2.3) также будет выполнено.

Из условия (1.2.5) на фронте фазового перехода, т.е. при , следует, что

; (1.2.12)

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

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

(1.2.13)

и его скорость

(1.2.14)

которая уменьшается со временем, т.е. по мере утолщения слоя твёрдой фазы.

Подставляя (1.2.12) в соотношение (1.2.11), получаем

; (1.2.15)

Теперь из равенств (1.2.11) и (1.2.15) находим все четыре константы:

(1.2.16)

Чтобы определить константу , надо воспользоваться условием Стефана (1.2.6) на фронте фазового перехода. Так как

,

то с учётом формул (1.2.10), (1.2.14) и (1.2.16) условие (1.2.6) приводит к уравнению

. (1.2.17)

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

В случае, когда начальная температура всех слоёв жидкости равна температуре фазового перехода, т.е. , из равенств (1.2.16) следует, что

; ;

; ,

и трансцендентное уравнение (1.2.17) принимает более простой вид

, (1.2.18)

где .

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

; .

1.3 Задачи с подвижными границами, решаемые разложением в ряд

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

, , ; (1.3.1)

, (1.3.2)

, , (1.3.3)

а на границе задано одно из следующих граничных условий:

, (1.3.4)

, (1.3.5)

. (1.3.6)

Неопределённость постановки состоит в том, что пока не оговаривается, что функции , , , , , известны или их нужно искать. Предположим только, что указанные функции и условия (1.3.2)-(1.3.6) таковы, что задача (1.3.1)-(1.3.6) имеет единственное аналитическое решение[6].

Напомним, что коэффициенты степенного ряда по x, являющегося решением уравнения (1.3.1), полностью определяется, если на одной из границ заданы два независимых условия. Пусть на границе заданы два условия: (1.3.2) и (1.3.3). Следовательно, функцию T ищем в виде ряда, используя только внутренние граничные условия (1.3.2) и (1.3.3), т.е. решение строим независимо от граничного условия (1.3.4) или (1.3.5) и (1.3.6). Для этого функцию разлагаем в ряд Тейлора в точке :

. (1.3.7)

Условия (1.3.2) и (1.3.3) непосредственно дают первые члены ряда (1.3.7). Далее, дифференцируя условие (1.3.2) по и используя основное уравнение (1.3.1), получим:

, , , (1.3.8)

аналогично из условия (1.3.3) имеем

, . (1.3.9)

В свою очередь, дифференцируя по условия (1.3.8) и (1.3.9) и вычисляя последующие дифференциалы, определим все производные :

; (1.3.10)

, (1.3.11)

Здесь и означают соответствующие производные по порядка и .

Таким образом последовательно определяются все коэффициенты ряда (1.3.7). Если известны значения функций ,,, то из граничного условия (1.3.4) или из (1.3.5) и (1.3.6) можно найти граничные функции, т.е. решить обратную задачу, а если известна одна из функций, например , то подставляя (1.3.7) в одно из условий (1.3.4)-(1.3.6), можно найти её, т.е. решить прямую задачу. Аналогично можно определить и неизвестные функции или .

В задаче Стефана: , , и задано одно из условий (1.3.4)-(1.3.6); требуется определить границу и температурное поле. Здесь — теплота фазового перехода, — коэффициент теплопроводности.

Прежде чем определить границу , отметим, что в данном случае ряд (1.3.7) имеет вид

. (1.3.12)

Для нахождения функции можно предложить по крайней мере два способа.

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

. (1.3.13)

Предположим, что . Пусть теперь задано, например, условие (1.3.5), при этом, не нарушая общности, полагаем . Дифференцируя условие (1.3.5) и используя основное уравнение (1.3.1), можно легко убедиться, что

, . (1.3.14)

Подставляя соотношение (1.3.11) в (1.3.14) и учитывая, что при граница , получим

,

а из выражений

,

непосредственно вытекает

Таким образом рекуррентно определяются все коэффициенты ряда (1.3.13). Для граничного условия (1.3.5) первые члены этого ряда имеют вид

. (1.3.15)

а для граничного условия (1.3.16) находим

Вопрос сходимости (1.3.15), (1.3.16) остаётся пока открытым. Предварительный анализ показывает, что ряды (1.3.15), (1.3.16) сходятся при очень малых значениях .

В случае задания граничного условия (1.3.4) можно легко убедиться, что все коэффициенты ряда (1.3.13) равны нулю. Действительно, дифференцируя (1.3.4) и используя основное уравнение (1.3.1), получим

.

С другой стороны, соотношение (1.3.10) при даёт

.

Но поскольку при граница , следовательно, и т.д.

Итак получено нулевое решение, что, конечно, невозможно. В чем же дело? Подозрение падает на ряд (1.3.13).

Вероятно, функция в точке не разложима в ряд Маклорена. Действительно, если обратиться к точному решению задачи

. (1.3.17)

Отсюда ясна невыполнимость соотношения (1.3.13).

Второй способ имеет более общий характер. Сначала ряд (1.3.7) с коэффициентами (1.3.10), (1.3.11) преобразуем к виду, удобному для оценки коэффициентов, затем — в операторный вид, удобный для дальнейшего решения.

В силу предположения об аналитичности справедливо

.

Теперь рассмотрим выражение

,

Откуда следует цепочка равенств

. (1.3.18)

Аналогично получим

(1.3.19)

Последовательно понижая порядок дифференцирования на единицу в правых частях выражений (1.3.18) и (1.3.19), находим соответственно соотношения (1.3.10) и (1.3.11). Расписывая, в свою очередь, выражения (1.3.10) и (1.3.11), имеем

, (1.3.20)

(1.3.21)

Непосредственной проверкой при с помощью выражений (1.3.18), (1.3.19) убеждаемся в справедливости (1.3.20) и (1.3.21). Пусть для верны формулы (1.3.20) и (1.3.21), докажем справедливость (1.3.20) для . Подставляя (1.3.20) и (1.3.21) в (1.3.18), получим

.

Таким образом, выражение (1.3.7) можно переписать в виде

(1.3.22)

при этом, не нарушая общности, полагаем . Очевидно, если наложить определённые условия на производные функции , то вид выражения (1.3.22) удобен для оценки коэффициентов рядов, но трудно получить решение относительно . Преобразуем выражение (1.3.22) дальше. Для этого соберём те члены, которые содержат одинаковый сомножитель , в результате имеем

. (1.3.23)

Теперь для определения подставляем (1.3.23) в одно из условий (1.3.4)-(1.3.6). Например, условие (1.3.4) при даёт

, (1.3.24)

и, следовательно,

. (1.3.25)

Нетрудно убедиться, что выражение (1.3.25) верно при . Далее можно показать, что справедливо соотношение

, .

Здесь формально считаем, что двойной факториал от отрицательных чисел равен единице.

Последнее соотношение с учётом (1.3.17) даёт следующее выражение:

,

поскольку

.

Следовательно, условие (1.3.24) можно переписать в виде

(1.3.26)

Известно, что

.

Отсюда из соотношения (1.3.26) получим известное трансцедентное уравнение для определения :

.

Для того чтобы вычислить ряд (1.3.23), отдельно рассмотрим входящую в него сумму по . Обозначая , и преобразуя в необходимых случаях индексы суммирования, находим

.

Но поскольку , имеем

.

Следовательно, выражение (1.3.23) можно переписать в виде

. (1.3.27)

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

В общем случае, когда заданы условия (1.3.2), (1.3.3) из выражений (1.3.8), (1.3.9), используя соотношение (1.3.18), (1.3.19), можно найти подобными выкладками следующую формулу:

. (1.3.28)

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАБОТЫ АРТЕЗИАНСКОЙ СКВАЖИНЫ

2.1 Плоскорадиальный фильтрационный поток воды к артезианской скважине (постановка задачи)

Водоносный пласт представляет собой прямой круговой цилиндр радиуса и толщины . Кровля, подошва и боковая поверхность его непроницаемы. Центральная скважина радиуса полностью вскрывающая пласт открытым забоем (гидродинамически совершенная скважина). Работает на режиме постоянного массового расхода (дебита) воды скважины . Во всех горизонтальных плоскостях фильтрационный поток воды в пласте будет одинаковым — одномерным радиальным; его называют плоскорадиальным потоком[7].

В данном случае математическая постановка задачи запишется в виде

, , — (2.1)

уравнение фильтрации, где — радиальная пространственная переменная, — время; и -начальный и некоторый конечный момент времени; — давление; — коэффициент динамической вязкости воды; и — коэффициенты пористости и проницаемости пласта.

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

: , (2.2)

Граничные условия, которые следует выразить через соотношения, содержащие нужную нам функцию давления :

, (2.3)

, — (2.4)

Массовый расход воды (его алгебраическое значение) через любую цилиндрическую поверхность представляется как

,

где — радиальная скорость фильтрации согласно закону Дарси,

и учитывая, что , для дебита воды на скважине найдём.

Тогда краевое условие (2.3) перепишется так:

, ; (2.5)

Граничному условию (2.4) будет соответствовать соотношение

, . (2.6)

Приведём ещё выражение для средневзвешенного пластового давления, которым в дальнейшем воспользуемся для приближённого аналитического решения сформулированной задачи (2.1), (2.2), (2.5), (2.6). Вычитая из начальных запасов суммарный за время отбор воды из пласта, получаем запасы, оставшиеся в нём на текущий момент времени (их также иногда называют текущими запасами):

, (2.7)

Здесь

(2.8)

Средневзвешенное по пласту давление, а вычитаемое в квадратных скобках представляет собой относительный отбор воды из пласта.

Перейдем теперь к безразмерным переменным и величинам:

, , , , (2.9)

Тогда исходная система уравнений (2.1), (2.2), (2.5), (2.6) представится в безразмерном виде следующим образом:

,, ; (2.10)

: , ; (2.11)

: , ; (2.12)

: , . (2.13)

Из соотношений (2.7), (2.8) найдём выраженное в безразмерной форме уравнение —

(2.14)

где

(2.15)

и

(2.16)

характеризует отнесённые к начальным пластовым запасам отбор и текущие запасы воды соответственно.

2.2 Приближенное аналитическое решение задачи о нестационарной плоскорадиальной фильтрации воды

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

Заметим, что для нелинейных задач типа (2.10-2.13), а также и более общих задач, которым соответствуют константные (или нулевые) начальные условия, например, вида (2.11), характерна конечность скорости распространения передней границы (фронта) расширяющейся области возмущения. Это важное обстоятельство существенно отличает такие задачи от задач, связанных с линейными дифференциальными уравнениями в частных производных параболического типа, для которых, как известно, скорость распространения возмущений бесконечна.

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

2.2.1 Первая фаза истощения водоносного пласта

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

, , ; (2.17)

; (2.18)

, , (2.19)

: ; (2.20)

, ; (2.21)

: , ; (2.22)

, ; (2.23)

. (2.24)

, (2.25)

где

(2.26)

Приближённое решение этой начально-краевой задачи будем искать, полагая в (2.17)

(2.27)

и рассматривая уравнение

, , ; (2.28)

при условиях, соответствующих выше приведённым (значок тильда указывает на то, что отыскиваемое решения отличается от точного; в дальнейшем он опускается).

Последовательно интегрируя (2.28) по R, получаем

и

,

откуда с учётом (2.22-2.24), находим

, (2.29)

Функция , согласно (2.25), (2.26), определяется из уравнения

,

или

(2.30)

Можно заметить, что иногда при практических расчётах давление на внешнем непроницаемом контуре кругового пласта принимается равным средневзвешенному пластовому.

,

или

(2.31)

Здесь через обозначена зависимость, которая приближённо представляет удовлетворяющую уравнению (2.30) нужную нам функцию .Выполняя выкладки и учитывая, что (), а также предполагая, что с некоторого момента времени можно уже считать , в результате получим

, . (2.32)

Отсюда, с учётом условия , которое соответствует (2.18), найдём момент окончания первой и начала второй фазы процесса истощения:

(2.33)

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

, (2.34)

Где продолжительность периода первой фазы, согласно условию , определяется как

(2.35)

Заметим, что это значение следует и из (2.33) при разложении в степенной ряд и пренебрежении малыми членами первой степени.

Функции -(2.32) и — (2.34) фактически соответствуют случаю вырождения скважины в сток.

Вместо (2.29) рассматриваем следующие две зависимости на временных промежутках и соответственно:

(2.36)

(2.37)

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

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

В отличие от дебита скважины , соответствующего зависимости (2.36), из (2.37) для расхода воды на радиусе следует такое выражение:

, (2.38)

Поскольку — монотонная возрастающая функция, то и будет монотонно возрастать, все ближе подходя к . Потребовав, чтобы относительное расхождение между и было меньше 0.01% (достаточно высокая), и используя (2.34) при , имеем следующие соотношения для определения :

или (2.39)

(2.40)

С учётом того, что , из (2.39) и (2.40) находим

(2.41)

На временном промежутке функцию будем задавать в виде

,

или иначе, говоря квадрат ее — представлять фрагментом квадратичной параболы, проходящей через точки и и имеющей производную при , которая совпадает с ; последняя, согласно (2.34), на всём отрезке является константой, равной 16, а из (40). Тогда получим, что

, (2.42)

В таблицах 2.1, 2.2, 2.3 приведены результаты расчетов изменения во времени координаты фронта возмущения по формулам (2.42), (2.34) и значений давлений на скважине, полученных с использованием зависимостей (2.37), (2.36) для ( м, м) и .

Аналитические выражения для давления на скважине таковы:

, (2.43)

— (2.42)

(2.44)

— (2.34)

При получении второй формулы пренебрегали малой величиной по сравнению с 1.

Таблица 2.1 — Результаты расчётов для и при ( м, м) и .

10-10	0,000204	0,9995
10-9	0,000234	0,9956
10-8	0,000447	0,9747
10-7	0,001281	0,9298
10-6	0,004005	0,8748
10-5	0,012650	0,8176
	0,02	0,7947
10-4	0,04	0,7601
10-3	0,126491	0,7025
10-2	0,4	0,6450
	1	0,5991

Таблица 2.2 — Результаты расчётов для и при ( м, м) и .

10-10	0,000204	1
10-9	0,000234	0,9996
10-8	0,000447	0,9975
10-7	0,001281	0,9930
10-6	0,004005	0,9875
10-5	0,012650	0,9818
	0,02	0,9795
10-4	0,04	0,9760
10-3	0,126491	0,9703
10-2	0,4	0,9645
	1	0,9599

Таблица 2.3 — Результаты расчётов для и при ( м, м) и .

10-10	0,000801	1
10-9	0,000810	0,9997
10-8	0,000894	0,9971
10-7	0,001497	0,9812
10-6	0,004079	0,9403
10-5	0,012673	0,8863
	0,02	0,8638
10-4	0,04	0,8294
10-3	0,126491	0,7718
10-2	0,4	0,7143
	1	0,6685

2.2.2 Вторая фаза истощения водоносного пласта

На этом временном периоде рассматриваем следующую систему уравнений:

, , ; (2.45)

: , (2.46)

: , ; (2.47)

: , ; (2.48)

Принимая в уравнении фильтрации (2.45)

(2.49)

и присоединяя уравнение (2.14), вместо (2.45-2.48) запишем приближенную математическую постановку задачи в виде

, , ; (2.50)

: , (2.51)

:; : , ; (2.52)

,, (2.53)

Где определяется по (2.35), а задаётся зависимостью (2.37) при и :

, (2.54)

Интегрирование обыкновенного дифференциального уравнения по пространственной переменной (2.50), с учётом граничных условий (2.52), и, если мы будем пренебрегать по сравнению, с 1, даёт

(2.55)

Неизвестную функцию времени выразим через давление на непроницаемом граничном контуре — контурное давление :

(2.56)

Тогда из (2.55) и (2.56) имеем

, , (2.57)

-решение, которое удовлетворяет начальному условию (2.54), (2.51).

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

(2.58)

и

, , (2.59)

Предельный момент времени находится из условия

,

где

, (2.60)

Отсюда получаем

(2.61)

Значения давлений на скважине и на внешнем непроницаемом контуре пласта, рассчитанных по формулам (2.60) и (2.58) для и дано в таблице 2.2. В таблицах 2.4, 2.5, 2.6 вычисленное по (2.61) равно 3,7681, 80,293 и 4,352.

Таблица 2.4 — Результаты расчётов для и (вторая фаза) ( м, м) и .

	0,5991	1
0,1	0,5917	0,9926
0,5	0,5141	0,9150
1	0,4216	0,8225
2	0,2516	0,6525
3	0,1016	0,5025
3,5	0,0341	0,4350

Таблица 2.5 — Результаты расчётов для и (вторая фаза) ( м, м) и .

	0,9599	1
0,1	0,9592	0,9993
0,5	0,9512	0,9913
1	0,9413	0,9814
2	0,9216	0,9617
3	0,9021	0,9422
3,5	0,8924	0,9325
10	0,7712	0,8113
30	0,4512	0,4913
50	0,2112	0,2513
80	0,0012	0,0413

Таблица 2.6 — Результаты расчётов для и (вторая фаза) ( м, м) и .

	0,6685	1
0,1	0,661	0,9926
0,5	0,5835	0,9150
1	0,4910	0,8225
2	0,3210	0,6525
3	0,1710	0,5025
3,5	0,1035	0,4350
4	0,0410	0,3725
4,3	0,0002	0,3317

В ходе написания данной главы было получено приближенное аналитическое решение поставленной задачи(задача с подвижной границей), которое далее я сравню с численным решением, полученным при помощи математического пакета Maple. Были найдены функции для давлений на скважине и подвижной границы. Для первой фазы, характеризующейся падением пластового давления, распространяющимся от скважины к непроницаемой границе(пока не достигнет её) были получены следующие функции:

,

,

Данные функции действительны на промежутке , значение было получено для заданных условий при помощи (2.39) и (2.40). Функции

,

определены на промежутке . была определена с помощью формул

.

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

Для второй фазы, характеризующейся падением давления по всему пласту, были найдены следующие функции, определяющие падение давление на контуре и на скважине в определённый момент времени:

,

— момент истощения скважины определяется по формуле:

.

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

2.3 Численное решение задачи с помощью пакета Maple. Сравнение с результатами приближённого аналитического решения

Дифференциальные уравнения лежат в основе математического моделирования различных, в том числе физических, систем и устройств [8,9]. Рассмотрим основные возможности для решения дифференциальных уравнений в частных производных в среде Maple, а также сравним результаты, полученные с помощью данного математического пакета с моим приближенным аналитическим решением.

В Maple имеется функция pdsolve для решения дифференциальных уравнений в частных производных. Она может использоваться в следующих формах записи:

pdsolve(PDE, f, HINT, INTEGRATE, build)

pdsolve(PDE_systern, funcs, HINT, other_options)

pdsolve(PDE_system, conds, numeric, other_options)

pdsolve(PDE_system, conds, type=numeric, other_options)

Эта функция введена вместо устаревшей функции pdesolve. В функции pdsolve используются следующие параметры:

PDE — одиночное дифференциальное уравнение с частными производными;

PDEsystem — система дифференциальных уравнений с частными производными;

conds — начальные или граничные условия;

f — неопределенная функция или имя;

funcs — (опция) множество или список с неопределенными функциями или именами;

HINT — (опция) равенство в форме HINT = argument, где аргумент может быть символом ‘+’, ‘*’, любым алгебраическим выражением или строкой ‘strip’;

INTEGRATE — (опция) задает автоматическое интегрирование для множества ODEs (если PDE решается при разделении переменных;

build — опция, задающая попытку построения явного выражения для неопределенной функции, независимо от общности найденного решения;

numeric — ключевое слова, задающее решение в численном виде;

other_options — другие опции.

Для решения дифференциальных уравнений с частными производными и его визуализации в Maple служит специальный инструментальный пакет PDEtool:

> with(PDEtools);

[PDEplot, build, casesplit, charstrip, dchange, dcoeffs, declare, difforder, dpolyform, dsubs, mapde, separability, splitstrip, splitsys, undeclare].

Ввиду небольшого числа функций этого пакета приведем их определения: build(sol) — конструирует улучшенную форму решения, полученного функцией pdsolve;

casesplit(sys, o1, 02,…) — преобразует форму дифференциального уравнения;

charstrip(PDE.f) — находит характеристическую последовательность, дающую дифференциальное уравнение первого порядка;

dchange(tr,expr,o1,02,…) — выполняет замену переменных в математических выражениях или функциях;

dcoeff(expr,y(x)) — возвращает коэффициенты полиномиала дифференциального уравнения;

declare(expr) и др. — задает функцию для компактного ее отображения;

difforder(a.x) — возвращает порядок дифференциала в алгебраическом выражении а;

dpolyform(systno_Fn,opts) — возвращает полиномиальную форму для заданной системы sys не полиномиальных дифференциальных уравнений;

dsubs(deriv1=a,…fexpr) — выполняет дифференциальные подстановки в выражение ехрг;

mapde(PDE,intoff) — создает карту PDE в различных форматах into с опциональным заданием имени неизвестной функции f;

separability(PDE, F(x,y,…), ‘*’) — определяет условия разделения для сумм или произведений PDE;

splitstrip(PDE,f) — разделяет характеристическую последовательность на несоединенные поднаборы;

splitsys(sys.funcs) — разделяет наборы уравнений (алгебраические и дифференциальные) на несоединенные поднаборы;

undeclare(f(x),…) и др. — отменяет задание функции для компактного ее отображения.

Одна из важнейших функций пакета DEtools — DEtools[PDEplot] — служит для построения графиков решения систем с квазилинейными дифференциальными уравнениями первого порядка в частных производных. Эта функция используется в следующем виде:

PDEplot(pdiffeq, var, i_curve, srange, о)

PDEplot(pdiffeq, var, i_curve, srange, xrange, yrange, urange, o)

Здесь помимо упоминавшихся ранее параметров используются следующие:

pdiffeq — квазилинейные дифференциальные уравнения первого порядка (PDE),

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

* animate = true, false — включение (true) или выключение (false) режима анимации графиков;

* basechar = true, false, ONLY — устанавливает показ начального условия на плоскости (х,у);

* basecolor = b_color — устанавливает цвет базовых характеристик;

* ic_assumptions — задание (в виде равенств или неравенств) ограничений на начальные условия для первых производных;

* initcolor = i_color — инициализация цвета кривой начальных условий;

* numchar = integer — задает число отрезков кривых, которое не должно быть меньше 4 (по умолчанию 20);

* numsteps = [integer1,integer2] — задает число шагов интегрирования (по умолчанию [10,10]);

* obsrange = true, false— прекращение интегрирования (true) при выходе отображаемой переменной за заданные пределы или продолжение интегрирования (false) в любом случае;

* scenе=[х,у,и(х,у)] — вывод обозначений координатных осей.

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

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

скважина уравнение моделирование производная

Таблица 2.7 — Сравнение расчётов для и (первая фаза) ( м, м) и c результатами, полученными в среде Maple

		аналитическое	Maple
10-10	0,000204	0,9995	0,9953
10-9	0,000234	0,9956	0,9897
10-8	0,000447	0,9747	0,9718
10-7	0,001281	0,9298	0,9249
10-6	0,004005	0,8748	0,8703
10-5	0,012650	0,8176	0,8119
	0,02	0,7947	0,7916
10-4	0,04	0,7601	0,7576
10-3	0,126491	0,7025	0,6986
10-2	0,4	0,6450	0,6412
	1	0,5991	0,5958

Таблица 2.8 — Сравнение расчётов для и (вторая фаза) ( м, м) и c результатами, полученными в среде Maple

	аналитическое	Maple	аналитическое	Maple
	0,5991	0,5562	1	0,9998
0,1	0,5917	0,5509	0,9926	0,9985
0,5	0,5141	0,4866	0,9150	0,9519
1	0,4216	0,3919	0,8225	0,8870
2	0,2516	0,2621	0,6525	0,7072
3	0,1016	0,1064	0,5025	0,5375
3,5	0,0341	0,0365	0,4350	0,4726

Как видно из представленных данных численное и приближенное аналитическое решение дают похожие результаты. Максимальное относительное расхождение не превышает 7,75%.

ЗАКЛЮЧЕНИЕ

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

Далее было проведено сравнение результатов приближённого аналитического решения с численным решением, полученным в среде Maple для второй фазы. Приближенное аналитическое решение показало хорошую точность в сравнении с численным решением, максимальное относительное расхождение не превысило 7,75%.

СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ

1. Минкин Е.Л. Исследования и прогнозные расчеты для охраны подземных вод / Е.Л. Минкин. М.: Недра, 1972.

2. Справочник проектировщика. Водоснабжение населенных мест и промышленных предприятий / Под ред. И.А. Назарова. М.: Стройиздат, 1977.

3. Справочное руководство гидрогеолога / Под ред. В.М. Максимова. Л.: Недра, 1979.

4. Куштанова Г.Г. Подземная гидромеханика / Г.Г. Куштанова, М.Н. Овчинников. Казань, 2010.

5. Карлович М.Л. Дифференциальные уравнения математической физики / М. Л. Карлович, М.Ю. Иванович Издательство МГТУ имени Н. Э. Баумана, 2002.

6. Федоров Ф.М. Граничный метод решения прикладных задач математической физики / Ф.М. Федоров. Новосибирск: Наука Сибирская издательская фирма РАН, 2000.

7. Малых А.С. Специальный курс математического моделирования фильтрационных потоков / А.С. Малых, Г.П. Цыбульский. М., 2002.

8. Дьяконов В.П. Компьютерная математика. Теория и практика / В. П. Дьяконов. М.: Нолидж, 2001.

9. Иванов В.В. Методы вычислений на ЭВМ / В.В. Иванов. Киев: Наукова Думка, 1986.

10. Дьяконов В.П. Maple 9.5/10 в математике, физике и образовании / В.П. Дьяконов// Серия «Библиотека профессионала». М.: СОЛОН-Пресс, 2006. 720 с.

11. Лейбензон Л.С. Движение природных жидкостей и газов в пористой среде / Л.С. Лейбензон. Москва. 1947.

12. Басниев К.С. Подземная гидромеханика / К.С. Басниев, И.Н. Кочина. М.: Недра, 1993.

13. Пыхачев Г.Б. Подземная гидравлика / Г.Б. Пыхачев, Р.Г. Исаев. М.: Недра, 1973

14. Леонтьев Н.С. Основы теории фильтрации / Н.Е. Леонтьев. М.: Недра, 2009.

15. Щелкачёв В.Н. Подземная гидравлика / В.Н. Щелкачёв, Б.Б. Лапук. М.: РХД, 2001.

16. Великанов М.А. Движение подземных вод в крупнозернистых грунтах / М.А. Великанов. М., 1945.

17. Чарный И.А. Подземная гидромеханика / И.А. Чарный. ОГИЗ, Гостехидат, 1948.

18. Басниев К.С. Нефтегазовая гидромеханика / К.С. Басниев, Н.М. Дмитриев. Ижевск: РХД, 2005.

19. Евдокимова В.А. Сборник задач по подземной гидравлике / В.А. Евдокимова, И.Н. Кочина. М.: «Недра», 1979.

20. Баренблатт Г.И. Движение жидкостей и газов в природных пластах / Г. И. Баренблатт, В.М. Ентов, В.М. Рыжик. М.: Недра, 1984. 211 с.

21. Азиз Х. Математическое моделирование пластовых систем / Х. Азиз, Э. Сеттари М., 2004.

22. Яковлев В.П. Гидродинамический анализ недр / В.П. Яковлев. Баку: ОНТИ, 1937.

23.Вебстер А. Дифференциальные уравнения в частных производных математической физики / А. Вебстер, Г. Сеге. Государ. Технико-Теоретическое издательство, 1933.

24. Жукова Г.С. Уравнения в частных производных: примеры, задачи, методы решения / Г.С. Жукова, Е.М. Чечеткина. М., 2003.

25. Арсенин В.Я. Методы математической физики и специальные функции / В.Я. Арсенин. М.: Наука, 1972.

26. Владимиров В.С. Уравнения математической физики / В.С. Владимиров. М.: Наука, 1971.

27. Комеч А.И. Практическое решение уравнений математической физики / А.И. Комеч. М.: Наука, 1993.

ПРИЛОЖЕНИЕ 1

Листинг программы численного решения поставленной задачи при условиях ( м, м) и первая фаза в математическом пакете Maple:

> PDE := diff(P(r, t), t) = 1 / r * diff(r * diff(P(r, t), r), r);

R[c] := 2e-4: Q[c] := 0.1: x := Q[c] / 2 / R[c]: R[f] := 4*sqrt(t);

BC := {P(r, 0) = 1, D[1](P)(R[c], t) = x, D[1](P)(R[f], t) = 0};

SOL := pdsolve(PDE, BC, numeric, spacestep=0.001):

> P := SOL :- value():

tau := [1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 2.5e-5, 1e-4, 1e-3, 1e-2, 6.25e-2]:

printf(«t, cttRf(t), мtPc(t), Паtnn»):

for i in tau do

s := P(R[c], i):

v := 4*sqrt(i):

printf(«%0.2et%0.7ft%0.7fn», op(2, s[2]), v, op(2, s[3])):

od:

SOL :- plot3d(r=R[c]..1, t=1e-10..6.25e-2, axes = boxed);

SOL :- plot(r=R[c], t=1e-10..6.25e-2, axes = normal);

ПРИЛОЖЕНИЕ 2

Листинг программы численного решения поставленной задачи при условиях ( м, м) и вторая фаза в математическом пакете Maple:

> restart:

PDE := diff(P(r, t), t) = 1 / r * diff(r * diff(P(r, t), r), r);

R[c] := 2e-4: Q[c] := 0.1: x := Q[c] / 2 / R[c]:

BC := {P(r, 0) = 1, D[1](P)(R[c], t) = x, D[1](P)(1, t) = 0};

SOL := pdsolve(PDE, BC, numeric, timestep=1e-2, spacestep=0.001):

> tau := [6.25e-2, 0.1, 0.5, 1, 2, 3, 3.5]:

printf(«t, cttPc(t), ПаtPk(t), Паnn»):

for i in tau do

s := P(R[c], i):

d := P(1, i):

printf(«%0.2et%0.7ft%0.7fn», op(2, s[2]), op(2, s[3]), op(2, d[3])):

od:

SOL :- plot(r=R[c], t=6.25e-2..3,labels = [t, Pc], axes = normal);

SOL :- plot(r=1, t=6.25e-2..3,labels = [t, Pk], axes = normal);

SOL :- plot3d(r=R[c]..2e-4, t=1e-10..3, axes = normal);