CC BY
52
УДК 539.3
Конечно-элементное моделирование напряженно-деформированного состояния горных пород с учетом ползучести
© Ю.И. Димитриенко, Ю.В. Юрин МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Предложена модель для расчета напряженно-деформированного состояния (НДС) осадочных горных пород с учетом их ползучести. Представлен алгоритм конечно-элементного решения трехмерной задачи ползучести, использующий конечно-разностные схемы метода Эйлера по времени. Разработано специализированное программное обеспечение, позволяющее строить компьютерные 3Б-модели областей горных пород по исходным сериям 2Б-изображений, полученных с помощью данных сейсморазведки, а также проводить конечно-элементный расчет изменения НДС горных пород во времени. Проведено численное моделирование напряженно-деформированного состояния горных пород на примере зоны из Астраханского нефтегазового месторождения. Установлено, что в одних точках происходит поднятие горной породы, в других — ее опускание. Скорость ползучести различных слоев различна — наибольшие значения скорости ползучести реализуются в глинистых слоях и в песчаных, заполненных жидкостью, которые обладают наиболее заметными свойствами ползучести. Разработанный алгоритм и программное обеспечение для численного моделирования показали себя достаточно эффективными и могут быть применены для исследования НДС горных пород.
Ключевые слова: горные породы, напряженно-деформированное состояние, ползучесть горных пород, метод конечных элементов, численное моделирование
Введение. Расчет напряженно-деформированного состояния (НДС) горных пород играет важную роль при решении ряда прикладных задач геофизики: прогнозирование безопасности эксплуатации подземных выработок в процессе добычи полезных ископаемых, более качественный прогноз коллекторских свойств осадочных горных пород, а также прогнозирование геодинамики горных массивов и возможных землетрясений [1-6]. Сложность проблемы моделирования НДС горных пород заключается в отсутствии достоверных данных о физико-механических свойствах горных пород, находящихся на больших глубинах, а также о массовых и поверхностных внешних механических нагрузках, действующих на изучаемую область горной породы. Только частично данные о физико-механических характерстиках отдельных слоев горных пород могут быть получены по данным сейсморазведки. Как правило, методы обработки сейсмоданных позволяют определить одну или две упругие константы каждого слоя горного массива — скорости продольных и (или) поперечных волн, по которым пересчитываются модуль упругости и коэффициент Пуассона.
Существуют методики экспериментально-расчетного определения упругих констант слоев горных пород с учетом их возможной анизотропии. Однако для корректного определения НДС необходима информация о реологических свойствах горных пород, главным образом о деформациях ползучести, которая развивается даже при относительно низких уровнях напряжений и является одной из причин геодинамического движения блоков горных пород [7]. Информацию о характеристиках ползучести слоев горных пород в настоящее время обычно получают только косвенно: по результатам испытаний аналогов — приповерхностных слоев. Число публикаций, в которых приведены расчеты НДС горных пород с учетом ползучести, невелико [8, 9].
Корректное задание граничных условий для задачи расчета НДС также составляет большую проблему. Для отдельного блока горного массива граничные условия зависят от общей тектонической обстановки в литосферной плите, определение которой, как правило, составляет не менее сложную задачу. Поэтому в настоящее время в основном используют модельное задание граничных условий, которое формулируется на основе обобщения экспериментальных данных — в виде задания тектонических горизонтальных напряжений, линейно изменяющихся по глубине горного массива [5].
В настоящей статье изложены результаты работ по разработке алгоритма численного решения 3D-задач расчета НДС в горных породах с учетом деформаций ползучести, а также по разработке собственного программного обеспечения, предназначенного для численного моделирования НДС в горных породах.
Постановка задачи расчета НДС горных пород с учетом ползучести. Рассмотрим краевую трехмерную задачу термоползучести на основе теории течения [10]:
Уо + рГ = 0;
о = С ••(£ - £ - £);
£с = Р(£с, о);
£ = def(u) = 2(У® и + (V® и)т); (1)
£с
,=с = 0
о • п у = ъе, и у = ие.
Здесь V — набла-оператор [11]; о — тензор напряжений; р — переменная плотность горной породы, зависящая от конкретного типа горной породы (известняк, песчаник, глина и т. п.); { = - ge 2 — вектор плотности силы тяжести; g — ускорение свободного падения;
4 С — тензор модулей упругости; £ — тензор малых деформаций; £С — тензор деформаций ползучести; £ = а(0-00) — тензор тепловых деформаций; Ж(£с, о) — дважды непрерывно дифференцируемая
в некоторой области О с К12 тензорная функция, описывающая модель скоростей деформаций термоползучести; и — вектор перемещений; ие — заданный вектор перемещений; п — вектор внешней нормали; 8е — заданный вектор напряжений на части поверхности тела , задающий тектонические напряжения; е2 — вектор базиса, ориентированный по направлению действия силы тяжести; а — тензор температурного расширения; ® — знак тензорного произведения; (■) — скалярное произведение.
Вектор напряжений 8е на границе блока горной породы, задающий тектонические напряжения, согласно [4] выберем линейно изменяющимся с глубиной горной породы:
8е = - а2п, (2)
где а = 5 • 104 Па/ м — экспериментальная константа, полученная по данным работы [4]; 2 = х3 — вертикальная координата (текущая глубина) горной породы; п — вектор нормали.
Слои горной породы будем считать изотропными, поэтому для них тензоры модулей упругости имеют следующий вид [10]:
4 С = ^Е ®Е + 2рА, где р — константы Ламе; Е — метрический
тензор; А — единичный тензор 4-го ранга.
Деформации ползучести большинства горных пород, как правило, обнаруживают нелинейную зависимость от напряжений и являются необратимыми, для их расчета применим теорию пластического течения [10—12], в которой дополнительно учтено влияние первого инварианта тензора напряжений на скорость деформации ползучести:
( Гт(-\\Ч \
Е(о) = -Л
1 +
( (о )
7 J
(3)
где 11 (о) = Е--о — первый инвариант тензора напряжений; аи — второй инвариант (интенсивность) тензора напряжений [11]; а^, л , г1, г2 — константы ползучести.
Вводя тензор упругих напряжений ое = 4С••£, тензор напряже-
с 4^1 с "6 4^1 6
ний ползучести о =- С ••£ и тензор термонапряжений о =- С••£,
тензор напряжений может быть записан в виде следующего разложения: о = ое + ос + о6 .
Температурное поле 6 будем считать известным, изменяющимся линейно с глубиной горной породы [6]:
г
0 = 00 +н (61 "бе), (4)
где 60 = 293 К; 01 = 305 К; Н = 4-103м.
Численный алгоритм решения задачи расчета НДС. Конечно-элементный расчет НДС конструкций с учетом деформаций термоползучести в настоящее время реализован в основных коммерческих программных пакетах, в том числе в АКБУБ [13, 14]. Однако в большинстве случаев для решения нелинейных уравнений ползучести применяют явные или неявные методы, основанные на методе Эйлера аппроксимации производных по времени. Метод Эйлера достаточно эффективен с точки зрения экономии оперативной памяти при проведении вычислений, но относительно затратен по времени вычислений и не всегда обеспечивает требуемую точность расчетов деформаций ползучести.
Для применения метода Эйлера тензорное соотношение ползучести (третье уравнение в (1)) рассмотрим в виде системы дифференциальных уравнений относительно тензора деформаций ползучести £с как функции от временного параметра т. Введем сетку по временному параметру: Лм = (т0 = 0,..., т м = Т) с шагом Бгр( Лм) =
= (Лт0 = т1 "тo,., = "т#-1).
Для начального шага (т = 0) на основе начального условия = 0 имеем следующую систему:
с
£
т=0
У-ое(0) +рГ = -V- о6; ое(0) = 4С ••£(0); £(0) = ёеГ(и(0));
(5)
ое(0) • п
= 8(й0) - о0
• п
и
(0)
= и (0) — и и ■
Из формулы для разложения тензора напряжений и указанного начального условия имеем следующие соотношения для тензоров на начальном шаге:
и
£с(0) = 0; о(0) = о£(0) + о6
(0) _ е(0ь е
(6)
Далее на т -м шаге (т еКN+1) имеем следующую вычислительную процедуру для метода Эйлера.
1. Вычисление тензоров деформаций и напряжений ползучести т -го шага:
[ес(т) = £с(т-1) +АттЖ(о(т-1));
I о<т) = "4 с ,.£с(т)
(7)
2. Решение краевой задачи и определение соответствующих тензоров упругих напряжений:
^(т) +рГ = -V • (ос( т) +о6); ое(т) = 4 С--£(т);
£(т) = ае^Ь;
ое(т) • п
= 8(йт) - (ос(т) + о6 ) • п
и
(т)
= и( т) — Ид
(8)
Вычисление тензора напряжений
(т) е(т) . с(т) , 6
о - о ^ ; + о ^ ; + о .
(9)
Метод конечных элементов для численного решения краевых задач расчета НДС. Рассмотрим задачу (8) для т -го шага метода Эйлера, где и^от) еЬ2(!и),8^) еЬ2(Еа). Будем предполагать, что
тензор модулей упругости С удовлетворяет условию положительной определенности [10], т. е. для всякого симметричного тензора
второго ранга £ справедливо неравенство( у > 0) £ ••4 С ••£ >у£ ••£
Умножив первое соотношение рассматриваемой системы (9) ска-лярно на векторное поле V е У0°(П), получим
| ёеГ( V) • о {т)с1У = | V • (V • (ое(т) + ос{т) + о6 ))ё¥ -
| ёеГ( V) ••(о с<-т) + о 6).
(10)
□
Применяя для первого интеграла в правой части формулу Остроградского и учитывая силовое граничное условие в системе, получим
и
| ёе^у) • •ое(т)йУ =| у • 8(ьт)йI -1 ёе^у) • • (ос(т) + о6) йУ.
^ V ) • • I -^гг
□ 1а а
Далее, подставляя определяющие соотношения, будем иметь
(11)
|ёеГ(у)••4С-ёе^и))йУ = { у^^й1 +{ёеГ(у)••4С••[ £«т) +1 IйУ. (12)
□ а V )
Тогда под слабым решением краевой задачи на « -й стадии метода Рунге — Кутты в пространстве Соболева Ж2(1)(а) будем понимать такую вектор-функцию и[т) е^2(1)(П), что если т) ё ^2(1)(П) — такая вектор-функция, что 7>£и (^(т)) = иЬт), то и[т) - ™(т) ё Z(□) и и«т) удовлетворяет последнему интегральному соотношению Уу ё 2(□) . Решение такой задачи существует и единственно [15].
Перейдем в последнем интегральном соотношении к матричной форме. Для этого запишем компоненты соответствующих тензоров и векторов в декартовых координатах [10—18]:
/ = (/1 /2 /1 )т
> та 3
/ = {иг\и{ът\$>{ът), у);
оГ" = о
/
е(т) _/ е(т)
'«11
е( т) 22
е (т) >«33
е (т) >«12
е( т)
Мз
е( т) ^ 23
1 = (/ /22 /11 /12 /13 /23)
се( т) = /„( т) с( т) с( т) ^с1 22 33
(т) «12
2с
(т)
«13
кб, /=|о е(т), о с(т), 1 с(т), 11
2в«т3 ))> к6,
(13)
где > Кп — пространство столбцов размерности п. Также введем матрицу дифференциальных операторов Б и матрицу модулей упругости С ё ЦК,6,6) [10]:
Б
0 0 I ' С1111 С1122 С1133 0 0 0
0 5 2 0 С2222 С 2233 0 0 0
0 0 53 ; С = С3333 0 0 0
5 2 51 0 С2323 0 0
а 3 0 51 сим. С1313 0
V0 5 3 5 2 ) V С1212
. (14)
Тогда интегральное соотношение можно записать в следующем виде: {/Оу)т СВи(т)йУ = | |(Бу)т С[с^ +е]йУ. (15)
а 2о а V )
Слабое решение краевой задачи будем искать на основе метода конечных элементов (МКЭ). В качестве конечного элемента будет исполь-
зован 10-узловой тетраэдр с квадратичной аппроксимацией. Введем далее конечно-элементное пространство (Он) в следующем виде:
Fh
(О h ) = {fh : О h ^ > M 3 Â ( * ) = ф * ( * )Ik , x e K e T, fK e > M30}, (16)
где ФK : K ^ L(M,3,30) — матрица базисных аппроксимирующих функций (функций формы); f* =( fKi fh fh • ••fK? fK\ fK3 )т — значения вектор-функции fh в узлах тетраэдра K e Th (называемые степенями свободы КЭ). Матрица ФK (.) имеет блочный вид:
ФK (x) = (Ф11 Ф22 Фзз Ф44 Ф12 Ф13 Ф14 Ф23 Ф24 Ф34 ); Ф, (x) = N, (x) E ;
Г L (x) ( 2L, (x) -1), i = j;
n,(x)=-L()L Л • • (17)
[4L, (x)L, (x), i * ,,
где E e L(M,3,3) — единичная матрица; Lt(x) — барицентрические координаты точки в тетраэдре K e Th, построенные по его вершинам,
i, j e N 5.
Пусть далее Euh, Ech — аппроксимация частей Eu и Ec границы Г = ЗО, полученная при триангуляции Th. В пространстве Fh (Оh ) выделим подпространство Zh (О h ) = { fh e Fh (О h )| fh ( x ) = 0, x e Euh сГ h}.
Тогда схемой МКЭ-поиска приближенного решения ослабленной задачи на s-й стадии метода Рунге — Кутты на регулярной триангуляции Th будем называть задачу поиска такой вектор-функции uh e Fh (он), что u%\x) = и{ьf(x) при x e Nd(Th ) oEuh (т. е. интерполирует вектор-функцию u^f на границе Euh ) и для любой вектор-функции vh e Zh (Оh ) справедливо соотношение
j (Dvh)т CDuSm)dV = j VhTSbm)dE + j (Dvh)т C^sj(+ l)dV. (18)
Он Ech ОН V y
Раскладывая интегралы по всей области Он на интегралы по конечным элементам и подставляя определение элементов пространства Fh (Оh ), получим
x (vk )т [Aku($k - fK ] = 0;
KeTh
Ak = j Bk tcb*dv ;
K
/к - | (Фк )т S(b'¡:)dЪ+\Бк ТС |е*(:) +11 ¿V, Вк (х) = ЭФК (х). (19)
дК п2пн К V )
Матрица Ак е ДМ,30,30) и вектор /к ->М30 представляет локальную матрицу жесткости и локальный вектор нагрузок соответственно. Пусть введена единая нумерация и : Тн ^{и>1,...,и30} степеней свободы триангуляции (инъективная функция, ставящая в соответствие каждому конечному элементу (КЭ) к е Тн кортеж номеров соответствующих ему степеней свободы), число которых в Тн обозначим N - 3еагё(Ш(Тн)). Тогда введем матрицу А? е £(М, Ы, Ы) и векторы /К, VI е > МЫ :
Ак -
Фк -
(А?)
I (Ак )) , р -о, (к), д -о, (к)
О, р, д еи(к)
)ы
(фк)
р -.)(фк ) , р -0(к) 0, р £и(к)
ЧЫ
, Ф-{/, V, }.
)
(20)
Построим далее на их основе следующую матрицу и векторы:
Ак - X Ак, Фк - X Фк . (21)
кеТн к еТн
Тогда схема МКЭ может быть записана в виде
V )т [А
К ~(:)К
Ы,
(*)
- 0.
(22)
Далее, поскольку vн е ), то (VК ) - 0, р е I, где I -- {р е о(к): У к е : :(дк п ЕЫЙ) Ф 0}. Кроме того, компоненты (и*:)К ) зафиксированы значениями иЬ: в узлах ^(Тн) пЕ„н. Но тогда указанные компоненты могут быть исключены из векторов VК и
; (: )К
(в результате получим векторы V , и]
К ,,(:)Ъ >Ш>Ы-д
, д - еагё(1)), а
известные компоненты в щ:)К могут быть перенесены в правую часть:
/'К -[(/'К ) -(/К)' -Е(АК)' (:)К)
V ре1 р
(23)
Тогда, исключив из матрицы АК строки и столбцы, а из вектора /К строки с номерами р е I (получив матрицу АК е £(М, N - д, N - д)
и вектор /а ё > К^ 4 ) и учитывая произвольность вектора vе ё > км д, приходим к разрешающей системе линейных алгебраических уравнений (СЛАУ) схемы МКЭ:
Ааи(т)а = /а. (24)
Семейство схем МКЭ при к ^ 0 сходится к слабому решению краевой задачи [13]. Решение разрешающей СЛАУ в данной работе осуществляется на основе метода сопряженного градиента.
Результаты численного моделирования. В качестве исходных данных для расчета НДС использованы данные сейсморазведки, представленные в форматах БЕСУ. Куб сейсмических данных разделен на последовательность двумерных изображений исследуемого горного массива. Типичное двумерное изображение сейсмоданных приведено на рис. 1. В качестве примера реализации разработанных
Рис. 1. Двумерное изображение сечения куба сейсмоданных Астраханского месторождения (данные ОАО ЦГЭ)
методик выбрано Астраханское нефтегазовое месторождение. Экспериментальные данные предоставлены ОАО «Центральная геофизическая экспедиция» (ОАО ЦГЭ). Для преобразования сейсмоданных в формат геометрических данных, позволяющих строить компьютерную 3Б-модель исследуемой области, разработано специализированное программное обеспечение, позволяющее проводить эти операции в полуавтоматизированном режиме, задавая поверхности раздела областей с помощью двумерных сплайнов. Компьютерная 3Б-модель исследуемого горного массива, представляющего часть (10*10 км в плане и 4 км по глубине) для Астраханского месторождения, построенная с помощью данного программного обеспечения, предсталена на рис. 2.
Использование собственного программного обеспечения для компьютерной генерации исследуемых 3Б-областей горного массива позволило сгенерировать высококачественные конечно-элементные сетки, применяя для этого в полуавтоматизированном режиме генератор Те1§еи, с последующим улучшением сетки в областях стыка областей.
Рис. 2. Компьютерная 3Б-модель исследуемого горного массива, представляющего часть Астраханского месторождения
В качестве КЭ выбран четырехузловой тетраэдр. Пример КЭ-сетки для исследуемой области горного массива Астраханского месторождения приведена на рис. 3. Число КЭ в сетке составило 9,456 млн элементов и 1,516 млн узлов.
Рис. 3. КЭ-сетка для исследуемой области горного массива Астраханского месторождения
Исследуемый горный массив в рамках компьютерной 3Б-модели был разделен на пять типов областей, которые на основании информации об их продольных скоростях звука были классифицированы как: 1) известняк; 2) глина; 3) песок; 4) песчаник; насыщенный жидкостью; 5) песчаник. Модуль упругости отдельных областей определен на основании данных сейсморазведки, плотность и коэффициент Пуассона — на основании характерных данных для представленных сред. Константы ползучести областей подобраны с использованием данных из [7]. Значения констант упругости и ползучести отдельных областей исследованного горного массива Астраханского месторождения приведены в таблице.
Значения констант упругости и ползучести отдельных областей исследованного горного массива
Материал Плотность р, г/см3 Модуль Юнга е, ГПа Коэффициент Пуассона v Коэффициент вязкости "Л, Гпа • с Показатель степени Показатель степени г2 Характерные напряжения стх, ГПа
Известняк 2,6 46,2 0,219 2 • 1010 1,75 0,2 105
Глина 2,35 3,4 0,405 105 1,1 0,2 105
Песок 1,6 0,179 0,405 2,5 • 106 1,3 0,2 105
Песчаник,
насыщен- 2,25 13,283 0,153 4 • 105 1,3 0,2 105
ный жидко-
стью
Песчаник 2,25 6,321 0,405 2 • 1011 1,3 0,2 105
Расчеты НДС горных массивов на основе разработанного алгоритма проведены с помощью собственного программного обеспечения, разработанного в НОЦ «СИМПЛЕКС» МГТУ им. Н.Э. Баумана, на базе программной платформы БМСМ. На рис. 4 показана ЗБ-картина вертикальных перемещений, м, ир = пг (/) - пг (0), куба горного массива за
счет деформаций ползучести. Перемещения распределены неравномерно по слоям, существуют зоны локальных максимумов вертикальных перемещений, связанных с сильной искривленностью слоев горного массива.
На рис. 5 приведены двумерные картины изменения вертикальных перемещений ир горного массива во времени. Несмотря на то
что характер распределения перемещений ир в целом на протяжении
исследуемого периода времени (1 год) остается примерно одинаковым, абсолютные значения перемещений заметно изменяются.
Рис. 4. Картина вертикальных перемещений ыгр куба горного массива (верхний слой массива не показан)
Рис. 5. Моделирование кинетики изменения поля вертикальных перемещений п2р исследуемого блока горной породы:
а — 1 месяц; б — 3 месяца; в — 7 месяцев; г — 1 год
Графики изменения вертикальных перемещений ы^ во времени в
разных точках горного массива представлены на рис. 6. Установлено, что в некоторых точках происходит поднятие горной породы, а в некоторых — опускание. На рис. 7 приведены графики изменения во
времени значений интенсивности деформаций ползучести гсы С течением времени деформации ползучести монотонно накапливаются во всех рассмотренных точках, но скорость ползучести их различна — наибольшие значения скорости ползучести реализуются в глинистых и в песчаных слоях, заполненных жидкостью, которые обладают наиболее заметными свойствами ползучести.
Ъиъ, мм 25
20
15
10
Точка 1
3
4 2
0 876 1752 2628 3504 4380 52566132 7008 7884 *,ч
Ъиз, мм 0
-0,13 -0,25 -0,38 -0,50 -0,63 -0,75 -0,88 -1,00 -1,13 -1,25
V
V
\ ч \ 5 /
ъ 1
4
0 876 1752 2628 3504 4380 52566132 7008 7884 *,ч
Рис. 6. Изменение во времени вертикальных перемещений ыгр исследуемого блока горной породы в шести различных точках
XX и а уу :
Расчетные ЗБ-картины горизонтальных напряжений а
Гпа, в рассматриваемом блоке горной породы представлены на рис. 8 Проведенные вычисления показывают, что напряжения а
XX и а уу
с увеличением глубины горной породы возрастают по модулю, являясь отрицательными по знаку (сжимающими). В приповерхностном слое напряжения ахх и ауу в некоторых зонах меняют знак и являются положительными, хотя абсолютные значения напряжений в этих зонах относительно невелики.
е£,% 0,25
0,20
0,15
0,10
0,05
Точка 1
4 р> 3
2 -
0 876 1752 2628 3504 4380 5256 6132 7008 7884 t, ч Рис. 7. Прогнозирование изменения интенсивности деформаций ползучести г°и исследуемого блока горной породы в его различных точках
Рис. 8. ЗБ-картины распределения полей нормальных напряжений а ^ и а уу в исследуемом блоке горной породы
На рис. 9 показана картина распределения интенсивности а и напряжений в рассматриваемом блоке горной породы. В целом интенсивность напряжений в горной породе возрастает с увеличением глубины, однако максимумы интенсивности напряжений достигают-
ся в областях с максимальным наклоном их границ (анитиклинали), поскольку в этих областях становятся положительными нормальные напряжения а . Касательные напряжения а Х2 достигают максимума на границах областей неоднородностей.
Рис. 9. ЗБ-картина распределения поля интенсивности напряжений аи в исследуемом блоке горной породы
Заключение. Предложена модель для расчета напряженно-деформированного состояния осадочных горных пород с учетом их ползучести. Предложен алгоритм КЭ-решения трехмерной задачи ползучести, использующий конечно-разностные схемы Рунге — Кутты различного порядка по времени. Разработано специализированное программное обеспечение, позволяющее строить компьютерные ЗБ-модели областей горных пород по исходным сериям 2Б-изображений, полученных с помощью данных сейсморазведки, а также проводить расчет изменения НДС горных пород во времени. Проведено численное моделирование НДС горных пород на примере зоны из Астраханского нефтегазового месторождения. Установлено, что в некоторых точках происходит поднятие горной породы, а в некоторых — ее опускание. Скорость ползучести разных слоев различна — наибольшие значения скорости ползучести реализуются в глинистых слоях и в песчаных, заполненных жидкостью, которые обладают наиболее заметными свойствами ползучести. Разработанный алгоритм и программное обеспечение для численного моделирования показали себя достаточно эффективными и могут быть применены для исследования НДС горных пород.
Исследование выполнено за счет гранта Российского научного фонда (проект №14-19-00847).
ЛИТЕРАТУРА
[1] Гущенко О.И. Сейсмотектонический стресс-мониторинг литосферы (структурно-кинематический принцип и основные элементы алгоритма). Докл. РАН, 1996, т. 346, № 3, с. 399-402.
[2] Жалковский Н.Д., Кучай О.А., Мучная В.И. Сейсмичность и некоторые характеристики напряженного состояния земной коры Алтай-Саянской области. Геология и геофизика, 1995, т. 36 (10), с. 20-30.
[3] Леонов Ю.Г. Напряжения в литосфере и внутриплитная тектоника. Геотектоника, 1995, № 6, с. 3-21.
[4] Ребецкий Ю.Л. Механизм генерации остаточных напряжений и больших горизонтальных сжимающих напряжений в земной коре внутриплитовых орогенов. Проблемы тектонофизики. К 40-летию создания М.В. Гзовским лаборатории тектонофизики в ИФЗ РАН. Москва, ИФЗ РАН, 2008, с. 431-466.
[5] Ребецкий Ю.Л. Механизм генерации тектонических напряжений в областях больших вертикальных движений землетрясений. Физическая мезо-механика, 2008, т. 11, № 1, с. 66-73.
[6] Гзовский М.В. Основы тектонофизики. Москва, Наука, 1975, 533 с.
[7] Van der Pluum В. А. Marble myionites in the Bancroft shear zone, Ontario, Canada: microstructures and deformation mechanisms. J. of Structural Geology, 1991, vol. 13, no. 10, pp. 1125-1135.
[8] Каюмов Р.А., Шакирзянов Ф.Р. Моделирование поведения и оценка несущей способности системы тонкостенная конструкция — грунт с учетом ползучести и деградации грунта. Ученые записки Казанского ун-та. Сер. Физико-математические науки, 2011, т. 153, № 4, с. 67-75.
[9] Стефанов Ю.П. Некоторые особенности численного моделирования поведения упруго-хрупкопластичных материалов. Физическая мезомеханика, 2005, т. 8, № 3, с. 129-142.
[10] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 4: Основы механики твердого тела. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2013, 624 с.
[11] Димитриенко Ю.И., Губарева Е.А., Юрин Ю.В. Асимптотическая теория термоползучести многослойных тонких пластин. Математическое моделирование и численные методы, 2014, № 4, с. 18-36.
[12] Димитриенко Ю. И., Губарева Е. А., Юрин Ю. В. Конечно-элементное моделирование процессов термоползучести на основе методов Рунге — Кут-ты. Наука и образование, 2015, № 3. doi: 10.7463/0315.0759406 http://technomag.bmstu.ru/doc/759406.html
[13] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 1: Тензорный анализ. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 367 с.
[14] Implicit Creep. URL: http://ansys.net/ansys/ papers/nonlinear/ con-flong creep.pdf
[15] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Москва, Бином, 2001, с. 363-375.
[16] Фалейчик Б.В. Одношаговые методы численного решения задачи Коши. Минск, БГУ, 2010, 42 с.
[17] Даутов Р.З., Карчевский М.М. Введение в теорию метода конечных элементов. Казань, КГУ, 2004, 239 с.
[18] Димитриенко Ю.И., Губарева Е.А., Сборщиков С.В. Конечно-элементное моделирование эффективных вязкоупругих свойств однонаправленных композиционных материалов. Математическое моделирование и численные методы, 2014, № 2, с. 28-49.
Статья поступила в редакцию 30.05.2015
Ссылку на эту статью просим оформлять следующим образом: Димитриенко Ю.И., Юрин Ю.В. Конечно-элементное моделирование напряженно-деформированного состояния горных пород с учетом ползучести. Математическое моделирование и численные методы, 2015, № 3, с. 101-118.
Димитриенко Юрий Иванович родился в 1962 г., окончил МГУ им. М.В. Ломоносова в 1984 г. Д-р физ.-мат. наук, профессор, директор Научно-образовательного центра «Суперкомпьютерное инженерное моделирование и разработка программных комплексов» МГТУ им. Н.Э. Баумана, заведующий кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор более 300 научных работ в области механики сплошных сред, вычислительной механики, механики и термомеханики композитов, математического моделирования в науке о материалах, вычислительной газодинамики. e-mail: dimit.bmstu@gmail.com
Юрин Юрий Викторович родился в 1988 г., окончил МГТУ им. Н.Э. Баумана в 2012 г. Аспирант кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор 10 научных работ в области вычислительной математики и механики. e-mail: yvyurin@yandex.ru
Finite element simulation of the rock stress-strain state under creep
© Yu.I. Dimitrienko, Yu.V. Yurin Bauman Moscow State Technical University, Moscow, 105005, Russia
A model for calculation of a rock stress-strain state considering creep is suggested. The algorithm for finite element solving the three-dimensional creep problem using finite-difference scheme of Euler's method with respect to time is presented. The specialized software is developed allowing the computer to build 3D-models of rock areas based on the initial series of 2D images, obtained with the seismic data, and to perform finite element calculation of variations in rock strain-stress state with time. Numerical simulation of rock stress-strain state was conducted on the example of a zone of the Astrakhan oil and gas field. It was found that there occurs rock mass rising in some points, and in the other points it can slope down with time. The creep rate of different layers is not the same — the highest values of the creep rate are realized in the layers of clay and sand, filled with fluid, which have the most notable creep properties. The developed algorithm and software for numerical simulation proved to be quite effective and can be applied to the study of rock stress-strain state.
Keywords: rock, stress-strain state, creep, finite element method, numerical simulation.
REFERENCES
[1] Guschenko O.I. Doklady RAN - Reports ofRAS, 1996, vol. 346, no. 3, pp. 399-402.
[2] Zhalkovskiy N.D., Kuchay O.A., Muchnaya V.I. Geologiya i geofizika - Geology and Geophysics, 1995, vol. 36 (10), pp. 20-30.
[3] Leonov Yu.G. Geotektonika - Geotectonics, 1995, no. 6, pp. 3-21.
[4] Rebetskiy Yu.L. Mehanizm generactii ostatochikh napryazheniy i bolshikh gori-zontalnikh szimayushikh napryazheniy v zemnoy kore vnutriplitovikh orogenov [The Mechanism of the Generation of Residual Stresses and Large Horizontal Compressive Stresses in the Earth Crust Intraplate Orogenes]. Problemy tec-tonofiziki. K 40-letiyu sozdaniya M.V. Gzovskim laboratorii tectonofiziki v IFE
RAS [On the 40th anniversary of creation the tectonophysics laboratory in the IPE RAS by M.V. Gzovskiy]. Moscow, IFE RAS Publ., 2008, pp. 431-466.
[5] Rebetskiy Yu.L. Fizicheskaya mezomekhanika - Physical mesomehanics, 2008, vol. 11, no. 1, pp. 66-73.
[6] Gzovskiy M.V. Osnovi tectonofiziki [Fundamentals of Tectonophysics]. Moscow, Nauka Publ., 1975, 533 p.
[7] Van der Pluum B.A. Marble myionites in the Bancroft shear zone, Ontario, Canada: Microstructures and Deformation Mechanisms. Journal of Structural Geology, 1991, vol. 13, no. 10, pp. 1125-1135.
[8] Kayumov R.A., Shakirzyanov F.R. Ucheniye zapiski Kazanskogo universi-teta.Ser. Physiko-Matematicheskie Nauki - Scientific notes of the Kazan University. Ser. Phys. and math, 2011, vol. 153, no. 4, pp. 67-75.
[9] Stephanov Yu.P. Fizicheskaya mezomekhanika - Physical mesomehanics, 2005, vol. 8, no. 3, pp. 129-142.
[10] Dimitrienko Yu.I. Mekhanika sploshnoi sredy [Continuum Mechanics]. In 4 vols. Vol. 4. Osnovy mekhaniki tverdogo tela [Fundamentals of Solid Mechanics]. Moscow, BMSTU Publ., 2013, 624 p.
[11] Dimitrienko Y.I., Gubareva E.A., Yurin Yu.V. Matematicheskoe modelirovanie i chislennye menody - Mathematical Modeling and Numerical Methods, 2014, no. 4, pp. 18-36.
[12] Dimitrienko Y.I., Gubareva E.A., Yurin Yu.V. Nauka i obrazovanie; elec-tronnyy nauchno-tekhnicheskiy zhurnal - Science and Education: Electronic Scientific and technical Journal, No. 03, March 2015. Available at: http://technomag.bmstu.ru/doc/759406.html, DOI: 10.7463/0315.0759406
[13] Dimitrienko Yu.I. Mekhanika sploshnoi sredy [Continuum Mechanics]. In 4 vols. Vol. 1. Tenzornyi analiz [Tensor Analysis]. Moscow, BMSTU Publ., 367 p.
[14] Implicit Creep. Available at: http://ansys.net/ansys/papers/nonlinear/ con-flong_creep.pdf
[15] Bakhvalov N.S., Zhidkov N.P., Kobelkov G.M. Chislennye metody [Numerical Methods]. Moscow, Binom Publ., 2001, pp. 363-375.
[16] Faleychik B.V. Odnoshagovye metody chislennogo resheniya zadach Koshi [One-Step Methods for the Numerical Solution to the Cauchy Problem]. Minsk. BGU Publ., 2010, 42 p.
[17] Dautov P.Z., Karachevskiy M.M. Vvedeniye v teoriyu metoda konechnykh ele-mentov [Introduction to the Finite-Element Theory]. Kazan. Kazan University Publ., 2004, 239 p.
[18] Dimitrienko Y.I., Gubareva E.A., Sborschikov S.V. Matematicheskoe modeliro-vanie i chislennye menody - Mathematical Modeling and Numerical Methods, 2014, no. 2, pp. 28-48.
Dimitrienko Yu.I. (b. 1962) graduated from Lomonosov Moscow State University in 1984. Dr. Sci. (Phys. & Math.), professor, head of the Department of Computational Mathematics and Mathematical Physics, director of the Scientific-Educational Center of Supercomputer Engineering Modeling and Program Software Development at Bauman Moscow State Technical University. Member of the Russian Academy of Engineering Science. Author of over 300 research publications in the fields of computational mechanics, gasdynamics, thermomechanics of composite materials, mathematical simulations in material science. e-mail: dimit.bmstu@gmail.com
Yurin Yu.V. (b. 1989) graduated from Bauman Moscow State Technical University in 2012. Post-graduate student at the Department of Computational Mathematics and Mathematical Physics at Bauman Moscow State Technical University. Author of 10 research publications in the field of computational mathematics and mechanics. e-mail: yvyurin@yandex.ru
