Содержание к диссертации
Введение
1.Анализ состояния вопроса и задачи исследования 10
1.1.Параметры и режимы работы дегазационных скважин, пробуренных с поверхности .10
1.2. Обзор современных методов теории фильтрации, применительно к задачам подземной газодинамики 16
Задачи исследования 19
2. Описание модеж процесса дегазации подработанного горного массива скважинами, пробуренными с поверхности .20
2.1.Постановка задачи моделирования 20
2.2. Описание принципов организации пакета программ для расчета движения газа в подработанном горном массиве методом конечных элементов 28
3. Экспериментальные исследования газодинамических характеристик дегазируемого горного массива и идентификация модели 40
3.1.Характеристика объекта исследований и методика экспериментов 40
3.2. Результаты экспериментов и идентификация модели .45
4. Исследование и выбор оптимальных параметров сквазшн.пробуренных с поверхности 55
4.1.Расчет движения газа в подработанном горном массиве без учета дегазационных скважин 57
4.2. Расчет движения газа к скважине, пробуренной с поверхности 62
4.3.Оценка экономической эффективности рекомендуемых параметров дегазации шахт скважинами,пробуренными с поверхности 77
Заключение 83
Список использованной литературы, 86
Приложения
- Обзор современных методов теории фильтрации, применительно к задачам подземной газодинамики
- Описание принципов организации пакета программ для расчета движения газа в подработанном горном массиве методом конечных элементов
- Результаты экспериментов и идентификация модели
- Расчет движения газа к скважине, пробуренной с поверхности
Введение к работе
"Основными направлениями экономического и социального развития СССР на I98I-I985 годы и на период до 1990 года"предусматривается увеличение добычи угля к 1985 году до 770-800 млн.тонн в год. Основной прирост добычи долкен быть достигнут за счет повышения производительности труда. Выполнение этого задания Коммунистической партии Советского Союза требует широкого применения новой техники и прогрессивной технологии разработки угольных месторождений.
Многочисленные исследования технико-экономической эффективности систем разработки угояьных пластов показали, что в условиях комплексной механизации очистных работ наиболее эффективными являются столбовые системы с обратным порядком отработки выемочных полей и погашением выработок по мере подвигания забоя. Они имеют ряд преимуществ, основными из которых являются:
Обзор современных методов теории фильтрации, применительно к задачам подземной газодинамики
Основы теории фильтрации для задач подземной газодинамики были заложены акад.Л.С.Лейбензоном. Большой вклад в разработку теории подземной газодинамики внесли И.АЛарный,И.Я.Попубаринова--Кочина, акад.С.А.Христианович. Подробное изложение основ теории фильтрации газов содержится в классических работах [15,16,17] и более поздних [18,19,20J .
Несмотря на прогресс теории фильтрации в нашей стране и за рубежом, применение ее для решения технологических задач в горном деле ограничено Это объясняется сложностью численной реализации моделей, отражающих реальную технологическую ситуацию. С развитием средств вычислительной техники диапазон применения методов математического моделирования расширяется. Начиная с 50-х годов такие модели строились на основе метода конечных разностей (МКР), Метод состоит в замене производных в соответствующих уравнениях конечными разностями в узлах прямоугольной сетки. При этом задача сводится к системе линейных уравнений с сильно разреженной матрицей большой размерности [21] . В качестве типичных работ по применению МКР к задачам расчета фильтрации газов укажем [22,23,24] .
Но для реальных технологических задач МКР чрезмерно усложняется. Основные трудности при этом возникают не при решении разностных уравнений, а при их построении. Для нерегулярных областей требуются нерегулярные сетки, а также единообразные способы аппроксимации физических законов. То есть желательно использовать ЭВМ не только при решении дискретных задач, но и при их формулировке. Поэтому, начиная с 60-х годов все большую популярность приобретает метод конечных элементов (МКЭ). Идея метода заключается в аппроксимации искомой функции на конечном множестве подобластей, называемых конечными элементами. Аппроксимирующие функции определяются через значения искомой функции в конечном множестве точек рассматриваемой области, называемых узлами. При этом задача сводится к поиску множества узловых значений искомой функции, доставляющих минимум некоторой величине, связанной с физической сущностью задачи. Для задач фильтрации минимизируется функционал связанный (по Эйлеру) с соответствующим дифференциальным уравнением. Процесс минимизации сводится к решению системы линейных алгебраических уравнений, относительно узловых значений искомой функции.
По мнению одного из теоретиков метода f.Стренга: "...новая идея сделала возможным использовать в большей степени мощности вычислительных машин (при построении аппроксимаций, при решении полученных дискретных задач и графическом изображении результатов), чем любые другие способы научных вычислений".
Исторически метод возник, как развитие матричных методов расчета конструкций. Первое формальное изложение метода связано с именами Тэрнера, Клаффа, Мартина, Топпа (1956). В последующие годы число публикаций по теории метода и его приложений стремительно нарастает [25] . В обширной литературе по методу выделяют два направления: математическое и инженерное. Наиболее полно отражает современное состояние первого и второго направлений, соответственно, йниги [26,27] и [28,29] . Метод широко применяется для решения задач подземной гидрогазодинамики [ЗО,ЗІ] , в том числе, для расчета фильтрации газов [32,33] .
Преимущества МКЭ при решении задач подземной газодинамики и геомеханики следующие:- возможность обработки областей произвольной формы;- единообразная обработка сложных граничных условий;- возможность подробного учета характеристик среды (неоднородность, анизотропия).
Оборотной стороной достоинств МКЭ является большой объем квалифицированной работы по программированию. Поэтому для расширения круга пользователей метода, стали создаваться пакеты прикладных программ для решения классов технологических задач. Подробный обзор идей и методов такого подхода к задачам фильтрации в рамках МКР можно найти в [34] . Развитие МКЭ дало возможность су-щественно повысить степень автоматизации решения задач механики сплошной среды. Это привело к изданию широкого спектра универсальных пакетов на основе МКЭ. Из крупнейших коммерческих пакетов такого рода назовем: NASTRAN (300), США [35] ; A SKA (30), ФРГ; MODUL-F (ЮО)» Фр, (в скобках - размер пакета в тыс.операторов). Имеются данные о создании подобных отечественных пакетов [ Зб] .
Современные тенденции автоматизации решения задач фильтрации в рамках МКЭ можно разбить по функциональным признакам на три направления. I.Автоматическая подготовка сеток:- априорная [37] ;- адаптивная [38,39,40] .
При этом широко используются диалоговые и графические средства [41,42] . 2.Применение "фронтального" метода сборки и решения линейной системы [43,44] с предварительной минимизацией ширины фронта [45,4б] . Предварительная перенумерация узлов сетки с
Описание принципов организации пакета программ для расчета движения газа в подработанном горном массиве методом конечных элементов
Расчет установившегося движения газа в подработанном горном массиве сводится к решению уравнения (2.5)в двумерной области 3г
На границе X области граничные условия видаїїі Ї2. Участки границы; о, - поток через границу. Решение уравнения (2.5) с граничными условиями (2.9) и (2.10) эквивалентно минимизации функционала [28І
Метод конечных элементов заключается в минимизации X на множестве узловых значений ф1 конечно-элементной сетки [ 28] Эта задача сводится к решению линейной системы где D J - квадратная матрица, традиционно называемая "матрицей жесткости", содержит данные о свойствах среды и геометрии области; (ф j - искомый вектор узловых значений функции f ; /pi - вектор, включающий данные о граничных условиях, источниках и стоках газа, традиционно называется "вектором нагрузки". Вариант процедуры МКЭ, реализованный в настоящей работе, состоит из следующих шагов.I.Разбиение области на треугольные элементы с линейной аппроксимацией искомой функции, вычисление и сборка по элементам матрицы жесткости и вектора нагрузки. 3.Решение системы (2.12) методом квадратного корня.. В цилиндрических координатах уравнение (2.5) имеет вид ( в осесимметрическом случае):Граничные условия для уравнения (2.10) выражаются (2.9) и формулой
Здесь h,j- направляющие косинусы вектора нормали к границе. krv» K&- главные значения тензора проницаемости. Вариационная формулировка задачи (2.13) с условиями (2.9) и (2.14) связана с рассмотрением функционала [28] : где F - расстояние от оси симметрии до центра тяжестиэлемента. Основные трудности, связанные с реализацией МКЭ вызваны большим объемом информации, подлежащей обработке. Поэтому, кроме модулей, реализующих собственно МКЭ, в состав пакета включены подпрограммы подготовки и проверки исходных данных и обработки полученных результатов.
Пакет предназначен для совместного решения уравнения (2.5) в плоской или осесимметричной постановке в произвольной двумерной области с граничными условиями 1-го, 2-го и 3-го рода, и системы алгебраических уравнений (2.10),(2.11), описывающих дегазационную систему. Основные эксплуатационные характеристики пакета следующие.Вычислительная среда ОСЕСКоличество модулей ......... 20Язык программирования ФОРТРАНОбъем текста программ 1,5 тыс.оп.
Максимальный размер загрузочного модуля 50 Кбайт. Все подпрограммы пакета можно разделить на три функциональные части:
I.Препроцессор. Включает модули: - триангуляции произвольной двумерной области; - ввода и проверки данных о сетке, подготовленных вручную; - ввода и проверки граничных условий и других данных; - интерполяции входных полей (проницаемости и плотности функции источника). 2.Основная часть. Реализует процедуру МКЭ: - формирование матриц элементов; - сборка по элементам глобальной матрицы жесткости; - решение линейной системы. Решение системы выполняется стандартной программой, входящей в $SP ФОРТРАНа (МСНВ ). При этом верхняя полупента целиком находится в оперативной памяти, что позволяет повысить быстродействие. 3.Постпроцессор. Включает модули: - построения изолиний решения на АЦПУ; - вычисление характеристик дегазационной системы (давления, расходы, содержание СН , коэффициент эффективности дегазации и др.);- другие служебные подпрограммы.
Опыт эксплуатации пакета показал, что подготовку сетки удобно производить заранее (рис.2.3). Обычно триангуляция области выполняется автоматически подпрограммой GtRID , использующей для генерации сетки семейство четырехугольных зон с квадратичной аппроксимацией координат [58] . Эта подпрограмма позволяет моделировать двумерные области, границы которых описываются кривыми не выше второго порядка. В программе осуществляется нумерация узлов и вычисляется ширина ленты матрицы жесткости. Полученные данные заносятся на внешний носитель. Недостатком программы является то, что для областей, имеющих разрезы, нумерация узлов производится неоптимально, в смысле ширины ленты. Заметим, чторесурсы (память и время), необходимые для решения системы (2.12) возрастает как квадрат ширины ленты, т2 . Для таких случаев предусмотрена возможность ввода данных о сетке, построенной вручную. При этом возникает проблема проверки входных данных, так как сетка может содержать до 900 элементов и 576 узлов. То есть проверке подлежит порядка 10 чисел. Эту задачу выполняет программа CARD, которая проверяет:- количество введенных узлов (элементов);- порядок нумерации узлов в элементе;- номера узлов (элементов);- повторное введение данных об элементе.
Результаты экспериментов и идентификация модели
Тогда, при Н = 335 и hi = 210 м, из (3.1) имеем кар = 97 м. Таким образом, по результатам замеров расхода газа по глубине скважины определена средневзвешенная по газу высота залегания спутников. Для дальнейших расчетов принято к ср = 100 м.
Положим, что газовыделение из сближенных пластов эквивалентно газовыделению из некоторой плоскости, расположенной на высоте І7Ср. В дальнейшем эта плоскость упоминается как " эквивалентный источник". Интенсивность газовыделения из него в зависимости от расстояния от лавы определялось по формуле [ 50 ] :где fff)" интенсивность газовыделения, кг/ с-м2;
У - общее газовыделенив из подработанного горного массива, кг/с ; X - расстояние от лавы, м;извлекала дополнительно 3,5 м3/мин метана или 22,6 % от общего расхода его в скважине. Результаты измерений расхода (табл.3.5) по длине скважины показали, что приток газа происходит на высоте от /ц до 2 , что соответствует расположению сближенных пластов (рис 3.2). Пласты, расположенные выше /? газа не дают.
Xm - абсцисса максимума f ( ) » ы; L/\ - длина лавы, и Для оценки Хт применялась эмпирическая зависимость Г50] ;константа А в дальнейшем использовалась для определения интенсивности газовыделения по формуле (3.2).
Для оценки проницаемости массива решалась.задача фильтрации газа к скважине в осесимметричной постановке. Подробное обсуждение постановки и результатов решения содержится в разделе 4 настоящей работы. Рассмотрим алгоритм определения проницаемости на основании сравнения опытных и расчетных данных.
Пусть поступление газа в скважину происходит из цилиндрической области радиусом R . Тогда, в пределах области положимгде г у - цилиндрические координаты;
Kr Kg.- радиальная и вертикальная проницаемость пород; оС - коэффициент пропорциональности. Пусть вертикальная проницаемость, К , зависит лишь от расстояния до выработанного пространства, % : где К«(о)- проницаемость на верхней границе выработанного пространства; В - коэффициент пропорциональности. Характер движения газа в подработанном горном массиве зависит от изменения проницаемости по глубине. О реальном движении газа в массива можно судить по кривой притоков в скважину по глубине (см.рис.3.2). Проницаемость пород в зоне обрушения, по данным МакНИИ [бз] ;Тогда в соответствии с (3.7), вертикальная проницаемость массива определяется коэффициентом В Выберем В так, чтобы расчетный приток газа в скважину был наиболее близок к экспериментальному. Отклонение кривой оценивалось как где 00" оценка близости экспериментальной и расчетнойкривой; &р 0.д " РаЧ8Тная и экспериментальная кривая притока газа по длине скважины; - координата длины скважины; 1 - длина участка скважины, лежащего ниже эквивалентного источника. Выбор В , доставляющего минимум v(f ) осуществлялся решением осесимметричной задачи фильтрации газа к скважине. Путем численных экспериментов выбрано
При этом расчетная кривая отличается от экспериментальной, в смысле 0 , на 8% (рис.3.6).Подбор коэффициента анизотропии, oL , выполнялся следующим образом.Задавшись оС - І, в результате решения фильтрационной задачи совместно с системой уравнений, описывающих дегазационную систему, получаем рабочую точку на характеристике вакуум-насоса (рис.3.7). Если расчетная точка лежит выше экспериментальной (R ) то расчетная скважина "легче" натурной, и коэффициент рс следует уменьшить. Если модельная скважина "тяжелее" натурной (случай Rf на рис.3.7), следует увеличить ос . Путем численных экспериментов, по методу золотого сечения, выбрано значение oL =1,8. При этом экспериментальный расход по скважине отличается от расчетного на 10%. Таким образом
То есть, проницаемость вдоль напластования примерно в 1,8 раза выше проницаемости поперек слоев. Полученные зависимости в дальнейшем использовались для определения параметров дегазационной скважины пробуренной с поверхности. І.В условиях шахты "Стожковская" одна скважина, пробуренная с поверхности может обеспечить снижение метановыделения из подработанных пластов на 62%,2.Расход метана, извлекаемого скважиной, превышает величину снижения метановыделения в горные выработки на 29%, что подтверждает справедливость мнения о 3.Предложена методика определения исходных данных для математического моделирования процесса дегазации подработанного горного массива скважинами,пробуренными с поверхности, и выполнена идентификация модели.
Расчет движения газа к скважине, пробуренной с поверхности
Путем моделирования движения газа к скважине дается оценка оптимального, в смысле эффективности дегазации расстояния от за 63боя скважины до кровли вынимаемого пласта (величины породной пробки).
Пусть для дегазации подработанного горного массива пробурена вертикальная скважина. Вследствие низкой проницаемости вышележащих пород, ограничимся рассмотрением части массива, находящегося ниже эквивалентного источника метана РуР2 (см.рисЛ.I) с интенсивностью (f(f) (рисЛЛ).
Рассмотрим движение газа в цилиндрической области, образованной вращением фигуры ОЗАВС вокруг отрезка 03 (рис.4.5). Для постановки осесимметричной задачи предположим, что интенсивность потока газа через участок АВ постоянна, то есть
То есть, если скважина находится на расстоянии ?Хт от лавы, замена переменного потока из эквивалентного источника в цилиндрическую область радиуса R= Хт постоянным, дает относительную погрешность 0=0,6%,
Чтобы определить размер области по радиальной координате, необходимо указать расстояние от оси скважины, Rm , на котором скважина не будет оказывать влияния на давление газа в массиве. По имеющимся данным [57] , для условий Донбасса, радиус,обеспечивающий эффективность дегазации 50$ составляет 40-50 м. На расстоянии вдвое большем влияние скважины будет исчезающе мало. Учитывая также то, что задание завышенного максимального радиуса не ухудшит качества результатов расчета, принято /?Пч = 90 м. Задавая проницаемость массива, учтем, что для осевой симметрии необходимогде Кг? Ку " коэффициенты проницаемости. При этом, учитывается лишь зависимость проницаемости,пород от глубины залегания. Эта зависимость задавалась на основе современных представлений о зонах нарушенности подработанного горного массива (рисЛ.б). Методика и результаты определения проницаемости в условиях, типичных для разработки пологих пластов Донбасса, приведены в третьем разделе настоящей работы. При этих допущениях, движение газа в ОЗАВС (см.рис.4.5) описывается уравнением
Здесь U - безразмерная функция давления где г, у - цилиндрические координаты; Р(г,у)- давление газа; Р - аїмосферное давление; Рт - максимальное давление газа в области. Рассмотрим граничные условия на контуре ОЗАВС (см.рис.4.6). На участке ОС поддерживается внутришахтное давление, которое по сравнением с давлением газа в пассиве можно считать равным атмосферному Алгоритм расчета режима работы скважины. Так как разрежение зависит от притока газа в скважину, решение проводилось по ите рационной схеме. I.Задание исходного приближения разрежения в скважине Uj . 2.Решение задачи (4.14),(4.16-4.19) в области ОЗАВС. 3.Вычисление притока газа в скважину на основании полученного в п.2 поля скорости фильтрации. 4.Вычисление очередного приближения UN 5.В случае достижения желаемой точности по UN (обычно --5%), переход на п.7. 6.Переход на п.2. 7.Вычисление содержания СН в извлекаемой смеси и эффективности дегазации. Пункты 2-7 реализованы в виде подпрограмм пакета, описанного во втором разделе настоящей работы. Типичные результаты решения для различных значений величины природной пробки приведены на рис«4.6 и 4.7. На этих рисунках можно наблюдать механизм появления подсосов воздуха в дегазационную систему: по маре приближения забоя скважины, 3 » к верхней границе выработанного пространства, ОС, скважина попадает в зону повышенной проницаемости пород. При этом величина области разрежения (показана штриховкой) резко увеличивается. Момент соприкосновения нулевой изобары с границей выработанного пространства, ОС, является началом поступления воздуха в скважину. Величина 03, при которой начинаются подсосы зависит, прежде всего, от проницаемости пород, интенсивности и расположения источников метана, кроме того, от состава и конфигурации дегазационной системы. Описанный во второй главе настоящей работы пакет программ, позволяет исследовать зависимость эффективности дегазации от величины 03 в широком диапазоне указанных условий. Ниже, в качестве примера, приведены результаты такого исследования, для условий определенных в третьем разделе настоящей работы. Для оценки эффективности по формуле (4.1), необходимо знать количество газа, поступающего в выработанное пространство без дегазации и при дегазации. Считая, что газ поступает в исследуемую цилиндрическую область с участка площади эквивалентного источника, ограниченного окружностью радиуса R , имеем где R - радиус цилиндрической области; о - интенсивность истечения газа из эквивалентного источника на расстоянии Х$ от лавы. Для определения количества газа # , попадающего в выработанное пространство при наличии скважины, проинтегрируем нормальную к участку ОС (см.рис .4.6) компоненту вектора скорости фильтрации (рис.4.8) по отрезку границы, прилегающему к выработанному пространству где Р - плотность газа; f - нулевая точка функции Vn{r) (см.рис.4.8). Минус в (4.21) объясняется тем, что, при выборе положительного направления у вверх от выработанного пространства истечению газа в выработки соответствует отрицательное значение вертикальной компоненты скорости. На участке положительной Vn имеют место подсосы воздуха из выработанного пространства в скважину (участок ОО1 на рис.4.7). Тогда, масса воздуха &в , попадающего в скважину равна Здесь г =у C V ) і где п 2 - координаты начала и конца отрезка с/г Зная нормальную к скважине компоненту скорости фильтрации Vr , получаем зависимость расхода по скважине, 0-с , от глубины, Z , (рис Л.9): где І - расстояние вдоль скважины от ее забоя; V - радиальная компонента скорости; F - радиальная координата центра тяжести граничного элемента, соответствующего t рбщий дебит газа, QG , складывается из притоков вдоль скважины QLQ И массы газа, &3 , поступающего через забой скважины: где Lg - длина необсаженной части скважины. Для притока через забой, имеем где Rc - радиус скважины; Vy - вертикальная компонента скорости. Тогда, содержание метана, С » в смеси,извлекаемой скважиной, есть Кривая Q.(О,полученная описанным способом (см.рис.4.9),используется для верификации модели, путем сравнения с экспериментальной кривой , получение которой описано в разделе 3 настоящей работы. Таким образом, зная поле скорости фильтрации V&-,y), порченное в результате численного решения задачи (4.14),(4.16-4.19), можно вычислить эффективность дегазации п и содержание СН в добываемой смеси, С . Полученные значения п и с зависят от величины породной пробки, n . Из физических соображений ясно, что с углублением скважины, содержание метана должно падать,вследствие увеличения подсоса воздуха, а величина расхода газа - возрастать, по мере того, как скважина проникает в область высокой проницаемости. Совместное влияние этих тенденций приводит к тому, что до некоторой глубины эффективность скважины возрастает (так как преобладает увеличение дебита), а затем из-за снижения содержания метана, начинает падать. Таким образом, должна существовать оптимальная, в смысле эффективности величина h В результате серии численных экспериментов, заключавшихся в решении описанной выше задачи при разных значениях п , установлено, что для приведенных выше условий, типичных для пологих пластов Донбасса, оптимальная величина породной пробки составляет 20 -г 30 мощностей разрабатываемого пласта (рис.4.10) [б4] . Для определения оптимального расстояния между скважинами, поступим следующим образом. Скважина оказывает воздействие на поле давлений газа в области радиуса R (см.рис.4.6, 4.7 ). Физически это воздействие выражается в отклонении линий тока газа. Так как без скважин изобары практически горизонтальны (см.рис. 4.3), о зоне влияния скважины можно судить по отклонению радиальной составляющей скорости от нуля. Таким образом, по результатам решения задачи ( 4.14 - 4.19 ) можно определить радиус влияния скважины, R . Для условий, приведенных в третьем разделе настоящей работы, в результате численного моделирования движения газа к скважине, получено
