Исследование напряженно-деформированного состояния неоднородного массива горных пород при действии гравитационных сил Текст научной статьи по специальности «Физика»

© В.О. Каледин, Т.В. Бурнышева, А.Б. Цветков, 2002

УДК 539.3: 622.83

В.О. Каледин, Т.В. Бурнышева, А.Б. Цветков

ИССЛЕДОВАНИЕ НАПРЯЖЕННО-Г ДЕФОРМИРОВАННОГО СОСТОЯНИЯ I НЕОДНОРОДНОГО МАССИВА ГОРНЫХ ПОРОД ^“'ПРИ ДЕЙСТВИИ ГРАВИТАЦИОННЫХ СИЛ

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

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

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

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

Параметры напряженно-

деформированного состояния рассматриваемой кусочно однородной области должны удовлетворять следующим группам уравнений [2]: уравнениям равновесия

Біу а - р =0; (1)

физическому закону

а = Вг;

(2)

кинематическим соотношениям Коши

г = Уи;

(3)

уравнениям неразрывности

Ь(г) = 0; (4)

граничным условиям и условиям сопряжения в напряжениях

I *

а п|г = р, а и|^+с= а и|^.о (5)

и перемещениях

I *

и | Е = и , и|^+о=и|^-о. (6)

Введем матричные обозначения: а = [ах ау аі тху тхі туі Ї - вектор напряжений; г =[гх гу гЕ Уху Ухе УуЕ Ї - вектор деформаций; В - матрица упругих характеристик изотропного материала [3]

1 л 1-л л 1-л 0 0 0

л 1-л 1 1 1-л 0 0 0

Е(1 - л) л 1-л л 1-л 1 0 0 0

(1 + л)(1 - 2л) 0 0 0 1-2 л 2(1-л) 0 0

0 0 0 0 1-2л 2(1-л) 0

0 0 0 0 0 1-2 л

где Е - модуль упругости материал; л - коэффициент Пуассона.

Численное решение краевой задачи (1)-(6) будем получать методом конечных элементов (МКЭ) в форме метода перемещений. Как известно [3], такой подход дает хорошо обусловленную систему линейных уравнений, но имеет недостатки, вытекающие из невозможности точного выполнения граничных условий в напряжениях и условий сопряжения (5) на границах раздела горных пород. Кроме того, учет распределенных по объему сил тяжести при дискретизации задачи может вносить в решение дополнительную погрешность [4].

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

Поэтому в настоящей работе предпринята попытка, с одной стороны, уменьшить погрешность решения на первом

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

[5].

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

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

п=11 2;

(гТа-

' -єТа)dО.- ЩиТ PdО.,

О О

где г - вектор полной деформации, а - вектор напряже-

(7)

^ £ 0= (е0х ,80 у ,е0х> 70 XI ,У0 у1, Г 0ху ) - вектор начальной деформации, первый интеграл - энергия деформаций,

второй интеграл - работа массовых сил Р .

Кусочно-неоднородную расчетную область представим в виде конечно-элементной модели. Поскольку в пределах однородной подобласти геостатическое напряжение пропорционально координате z, то аппроксимация перемещений должна иметь порядок полинома не менее двух. Кроме того, неправильная форма включений требует значительной неравномерности сетки, что связано с неправильной формой конечных элементов. Отсюда вытекает требование постоянства якобиана преобразования, отображающего искаженный конечный элемент на элемент канонической формы. Таким образом, следует выбрать субпараметрический элемент с базисными функциями второго порядка. Этим требованиям отвечает десятиузловой тетраэдр с дополнительными узлами на серединах ребер (рис. 1). Как показали численные эксперименты, такой выбор гарантирует устойчивость результатов при весьма неравномерной сетке.

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

Рис. 2. Граница раздела двух сред

Рис. 1. Субпарамет-рический тетраэдр с десятью узлами

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

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

4е) ( 1 - л(е))

а(е) = ах

0х (е) (е)

л р

г(е) ( 1 - л(е))

а(е) =

0У .,(е) р(е)

г(е)

0і

р(е)

(8)

(е)

т(е) т(е)

_(е) = 1ху _(е) = 6уг _

0ху р(е) , Оуг р(е) , 0іх р(е)

(е) = Чх

(9)

где р

(е)_

геостатическое давление в центре элемента е ;

(е) (е) (е) (е) (е) (е)

а0х ,СТ0у >ст0г ,Т0XI ,Т0У1 ,Т0хУ - относительные компоненты напряжений; ^(е) - коэффициенты Пуассона.

Геостатическое давление в центрах элементов рассчитывалось по формуле:

р(Є) =Ё Рі ‘ Ні

(10)

і=1

где Ь - высота однородного слоя массива расположенного над расчетной точкой, р1 - плотность г-го однородного слоя. Методика расчета давлений более полно рассмотрена в статье [7].

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

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

На границе раздела двух сред непрерывными остаются шесть функций:

деформации в плоскости, касательной к поверхности раздела,

Ss = , 7= 7^ , (11)

и составляющие напряжений, действующие по нормали к ней:

СТ1 = т1 = т2

и п и п ’ ^п

т1 = г2

(12)

Остальные компоненты напряжений и деформаций разрывны.

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

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

Представим все составляющие вектора напряжений и вектора деформаций на границе раздела сред в произвольном базисе {х, у z} через непрерывные составляющие в базисе, связанном с границей {і, ґ, п}, где 5, t - единичные векторы, лежащие в касательной к границе раздела плоскости, и - единичный вектор нормали.

Имеем

о*

О п

Оп

У 5,

т,п Т 5п _

и

0{п5!}=

'Є* "

ЄҐ

Єп =

уя,

у ,п

у яп

Е цЕ

1-Ц2 1-Ц2

цЕ Е

1-Ц2 1-Ц2

0 0

0 0

0 0

0 0

5,}

' 1 0

0 1

0 0

2 Ц+Ц 1-Ц2 2 Ц + Ц 1-Ц2

0 0

0 0

0

0

Е

2(1+ц)

0

Ц+Ц

1-л2

Ц+Ц2

1-Ц2

0

1

0

0

2

00

00

00

00

(13)

1-3ц2-2ц3 (1-Ц2)Е 0 0

2(1+ц)

Е

0

2(1+ц)

Е

аҐ

Оп

У 5,

т,п

=WUíns,}

(14)

где через Р и W обозначены матрицы преобразования. Перейдем к базису {х,у^}:

р[хуг}_С . р1п5<} (15) и

Q{xУz} =С ' . ^п5,}

Рис. 3. Изопараметриче-ский шестигранник с восьмью узлами

где С и С' - матрицы перехода от произвольного базиса {х, у, z} к базису на границе {і, ґ, и}.

Разобьем трехмерную пространственную неоднородную область на конечные элементы так

чтобы каждый отдельный элемент содержал однородный материал. В качестве конечного элемента возьмем шестигранник, вершины которого выберем в качестве узлов (рис. 3).

В пределах этого элемента аппроксимируем каждую непрерывную компоненту вектора напряжений трилинейной функцией:

О

(ї,Ч,0 = [^Е N 2 Е ... Ж8Е ]>

“ху

(17)

где (£,ц,£) - базисные функции [6].

Перейдем к непрерывным компонентам тензоров в базисе {х, у, z}:

о(£,гі,С) = МСРи{5,и}. (18)

Таким образом, зная значения непрерывных функций и{5,п'} в узлах данного конечного элемента (КЭ), можно получить значения всех напряжений данного конечного элемента.

Используя первую методику, мы рассчитали несогласованные значения вектора напряжений внутри КЭ. Согласованные напряжения ст(£,^,£) на элементе определим, минимизируя невязку между согласованными и несогласованными напряжениями:

Ф(о) = |(о - Мо)2dV ^ тіп.

(19)

V

Возьмем производную от функционала по и приравняем ее к нулю. С учетом равенства (18) получим:

| PTCTNTNCPdV и{^”} = | РТСТЫТ а(!У . (20)

V V

Узловые значения согласованных напряжений определяются решением системы линейных алгебраических уравнений (20).

По предложенной методике рассчитывалось поле согласованных напряжений в неоднородном массиве с пластом в виде антиклинали (рис. 4). Расчетная область представляла собой вертикальный разрез глубиной 1000 м и шириной 20000 м. С учетом симметрии половина области (шириной 10000 м) разбивалась на 400 элементов.

сг

5

X

т

я,

т

5И

т

X

0

0

0

0

0

0

0

X

0

0

0

0

0

0

X

Г

яп

Исследовалось влияние коэффициента Пуассона пласта на напряжения в области. В первом случае задавались следующие характеристики пласта: Е = 0.45 МПа, G = 0.202 МПа, ц = 0.115, р = 2.8 г/см3. Коэффициент Пуассона вмещающей толщи был в 2 раза меньше коэффициента Пуассона пласта, модули Юнга и плотности совпадали. На рис. 5 представлены рассчитанные значения напряжений компоненты ох : точками - не-

0.00

-100.00

-200.00

-300.00

-400.00

-500.00

-600.00

-700.00

500.0

1000.0 1500.0

несогласованные

напряжения

согласованные

напряжения

согласованные, сплошной кривой - по уточненной методике.

Рис. 4. Вид половины поперечного сечения области с пластом в виде антиклинали

Рис. 5. Значения напряжений О х .

Рис. 6. Изолинии напряжений компоненты О х Рис. 7. Значения напряжений О х .

Рис. 8. Изолинии напряжений компоненты О х .

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

На рис. 6 представлены изолинии компоненты О х , градиент которой имеет (как видно из рисунка) горизонтальную составляющую. Горизонтальная составляющая антиградиента направлена в сторону зоны уменьшения напряжений, т.е. под свод антиклинали.

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

На верхней границе пласта относительная погрешность несогласованных напряжений о 2 по сравнению с согласованными напряжениями в узлах элементов равнялась 61 %, на нижней границе - 10 %.

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

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

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

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

В качестве примера использования основной методики, основанной на использовании классического метода конечных элементов (аппроксимация век-

z

1000 2000 3000 4000

1000 2000 3000 4000

1000 2000 3000 4000

тора перемещения многочленом второго порядка, КЭ - десятиузловой тетраэдр), рассмотрим задачу вычисления поля напряжений пласта в виде антиклинали (см. рис. 4.). На область действует распределенная сила в виде собственного веса. Исследовалось изменение областей относительного растяжения и сжатия в зависимости от варьируемой плотности. Расчет производился при следующих параметрах мате-

Рис. 9. Линии уровня О,г / p при Еі/Е2=50000/100000, рі/р2=2.5/2.5, ц.і/ц.2=0.3/0.3

Рис. 10. Линии уровня ф / p при Е1/Е2=50000/100000, р1/р2=3/2, ц.1/ц.2=0.3/0.3.

Рис. 11. Линии уровня Оz / p при Е1/Е2=50000/100000, р1/р2=2/3, ц.1/ц.2=0.3/0.3.

риалов: Е1/Е2 = 50000/100000, рі/р2 = =(2,5/2,5, 3/2, 2/3), Ц1/Ц2 = 0,3/0,3. Индекс 1 относится к вмещающей толще, индекс 2 - к пласту.

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

Из рис. 9 видно, что справа и слева под пластом симметрично расположены две области относительного растяжения. При соотношении рі/р2=3/2 (рис. 10) обе области относительного растяжения смещаются ближе к куполу антиклинали, чем при одинаковой плотности пласта и вмещающей толщи. Если соотношение рі/р2=2/3 (рис. 11), то области относительного растяжения смещаются ниже и ближе к краям и в глубину. Во всех случаях ближе к краям наблюдаются области относительного сжатия с высоким градиентом относительных напряжений. На рис. 9-11 они выделены сплошным белым цветом у боковых границ области. Вероятно, быстрое изменение напряжений в этих зонах указывает на наличие особенности в решении задачи теории упругости в угловой точке границы пласта.

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

--------------- СПИСОК ЛИТЕРАТУРЫ

1. Ластовецкий В.П. Применение метода конечных элементов для оценки напряженного состояния и прогноза фильтрационно-емкостных свойств нефтяной залежи // Краевые задачи и математическое моделирование: Сб. трудов 3-й межотраслевой научной конференции/ Новокузнецк: НФИ КемГУ, 2000. - С. 37-41.

2. Ломакин В.А. Теория упругости неоднородных тел. - М.: Изд-во МГУ, 1976. -368 с.

3. Зенкевич О. Метод конечных элементов в технике. -М.: Мир, 1975. - 541 с.

4. Каледин В.О. О математическом моделировании физических полей в сложно построенных массивах горных пород. -Горный информ.-аналит. бюллетень. - М.: Изд-во МГГУ, №10, 2000. - С. 15-19.

5. Оден Дж. Конечные элементы в нелинейной механике сплошных сред. - М.: Мир, 1976. - 464 с.

6. Сегерлинд Л. Применение метода конечных элементов. - М.: Мир, 1979. - 392 с.

7. Каледин В.О., Цветков А.Б. Программа нормирования напряжений для пакета решения трехмерных задач геомеханики // Информационные технологии в экономике, промышленности и образовании: Вып. 2 / М.: «Электрика», 1999. -С. 13-19.

КОРОТКО ОБ АВТОРАХ -------------------------------------------------------------------------------------------

Каледин Валерий Олегович, Бурнышева Татьяна Витальевна, Цветков Андрей Борисович — Новокузнецкий филиал - институт Кемеровского государственного университета.