CC BY
56
ПРОБЛЕМЫ . ФИЗИЧЕСКИХ ПРОЦЕССОВ ГОРНОГО . ПРОИЗВОДСТВА :
...© В.О. Каледин, 2000
УДК 539.3:518.6
В.О. Каледин
О МАТЕМАТИЧЕСКОМ МОДЕЛИРОВАНИИ ФИЗИЧЕСКИХ ПОЛЕЙ В СЛОЖНО ПОСТРОЕННЫХ МАССИВАХ ГОРНЫХ ПОРОД*
Введение
Интерпретация данных полевых измерений, в первую очередь - в задачах разведочной геофизики, сопряжена с построением математических моделей физических полей различной природы. Такие модели традиционно строятся путем решения краевых задач для уравнений математической физики, что в случае сложно построенного массива горных пород с включениями неправильной формы обычно требует трудоемкого численного исследования. К настоящему времени накоплен хорошо разработанный аппарат решения краевых задач на основе сеточных методов (разностные схемы, конечные элементы) и методов типа граничных интегральных уравнений. Использование этих методов невозможно без использования вычислительной техники; этим определяется первостепенная важность программной реализации соответствующих математических методов с учетом ориентации на конкретную область их применения.
В настоящей работе рассматриваются математические модели полей, построенные в рамках эллиптических крае-
Рис. 1. Пример расчетной области в сложно построенном массиве
вых задач. К ним относятся задачи теории упругости и электропроводности, на базе решений которых решаются задачи расчета горного давления, прямые задачи электроразведки постоянным током и электроразведки низкочастотным переменным током, в том числе с измерением магнитного поля. В основу решения перечисленных задач положен сеточ-
ный метод конечных элементов [1,2].
1. Задачи о деформировании массива при гравитационных и тектонических воздействиях
Моделирование механических напряжений (горного давления) в массиве горных пород наталкивается на сложности в описании механических свойств. Общеизвестно, что свойства горной породы на глубине ее залегания невозможно установить с исчерпывающей точностью. Этим вызвано разнообразие моделей - линейно и нелинейно-упругих, учитывающих пластические и вязкостные свойства и других моделей деформируемого тела, для описания физико-механических свойств которых требуется знание экспериментально определенных констант, которые на практике известны лишь с очень невысокой точностью. Поэтому расчеты полей напряжений в сложно построенных массивах являются в значительной степени приближенными, и при построении модели требуется исследовать, насколько существенно влияет на результат учет того или иного фактора.
Учитывая, что решение нелинейных задач, в том числе - задач о пластическом и вязком деформировании, может быть построено в виде последовательности решений задачи о линейно-упругом деформировании, представляется целесообразным рассмотреть решение задачи о деформировании сложно построенного массива в рамках теории упругости кусочно-однородной среды. Наиболее существенными факторами, безусловно подлежащими учету, в этом случае являются: различие модулей упругости слагающих пород; форма внутренних неоднородностей; различие плотности пород; форма рельефа дневной поверхности; наличие воздействий тектонического происхождения (сил, сжимающих или растягивающих массив в горизонтальном направлении).
Выделим объем участка земной коры (рис. 1), ограниченный дневной поверхностью сверху, горизонтальной плоскостью снизу и четырьмя вертикальными плоскостями с боков, и отнесем его к декартовой системе координат XYZ, в которой ось Z направлена вниз от дневной поверхности, а две горизонтальные оси параллельны границам полученной расчетной области. На образовавшихся границах необходимо поставить граничные условия в напряжениях либо в перемещениях. Вид граничных условий определяется содержанием рассматриваемой задачи; целесообразно рассмотреть, по крайней мере, следующие условия: а) все искусственно введенные границы после деформации остаются плоскими, а размеры расчетной области не изменяются (перемещения граничных точек лежат в плоскости границы, перемещения вдоль нормали отсутствуют); б) границы остаются плоскими, но размеры области изменяются при деформации; в) имеет место депланация границ. Воздействиями на объект расчета будут: собственный вес массива (гравитационные воздействия) и давление прилежащих объемов (тектонические воздействия). Свойства пород, слагающих массив, будем считать известными в некотором приближении и характеризуемыми их плотностью и модулями упругости первого и второго рода, т.е. массив считаем
*Работа выполнена при финансовой поддержке Российского Фонда фундаментальных исследований (грант № 00-07-90141)
кусочно-однородным и кусочно-изотропным. При наличии соответствующих кадастровых данных нетрудно обобщить свойства массива и на случай анизотропных слагающих пород.
Как известно, поставленная таким образом задача теории упругости сводится к отысканию минимума функционала Лагранжа (потенциальной энергии) вида [2]:
П =11 ст(е) • гdV - | и ■pg • dV -1 и • q • dV , (1)
2 V V £
где и(ху^) - искомое поле перемещений, Е - деформации, п(е) - напряжения, выраженные через деформации из обобщенного закона Гука, р - плотность, g - ускорение силы тяжести, д - давления на границе; в произведениях тензорных величин предполагаются свертки. Функционал Лагранжа обладает свойством выпуклости, благодаря которому подстановка в (1) искомого поля перемещений в виде линейной комбинации некоторых базисных функций и варьирование по коэффициентам разложения приводит к хорошо обусловленной системе линейных алгебраических уравнений с симметричной положительно определенной матрицей. Таким образом, задача сводится к выбору системы базисных функций и решению системы уравнений высокого порядка.
Следуя обычной процедуре метода конечных элементов, представим искомое поле перемещений сплайном, а коэффициенты разложения в этом случае будут иметь смысл перемещений узлов. Выбор порядка и степени гладкости сплайна и определяемый этим выбор формы конечных элементов нетривиален. Следует ожидать, что гладкость решения нарушается на границах неоднородностей. В самом деле, в силу различия модулей упругости напряжения, действующие вдоль границы раздела, имеют разрыв, поскольку деформации границы непрерывны, и наоборот - деформации нормали в точках, примыкающих к границе раздела с разных сторон, различны. Таким образом, имеет смысл использовать сплайн нулевой степени гладкости, обеспечивающий только непрерывность самих искомых функций, а не их производных.
Порядок сплайна при учете гравитационных воздействий целесообразно выбирать не ниже второго; в противном случае, как показывают численные эксперименты, точность решения сильно зависит от геометрии сетки конечных элементов. Проведенные исследования позволили также утверждать, что изопараметрические конечные элементы [1] мало пригодны для решения рассматриваемого класса задач. Оба отмеченных обстоятельства связаны с тем, что в точном решении данной задачи напряжения растут с увеличением глубины, причем эта зависимость близка к линейной. В то же время симплексные элементы первого порядка аппроксимации дают кусочно-постоянное поле напряжений, а изопараметрические второго порядка - дробнорациональную зависимость напряжений от глубины. Поэтому как симплексные, так и изопараметрические элементы не дают решения с приемлемой точностью. Наиболее эффективным оказывается использование субпараметрических элементов второго порядка в форме тетраэдра.
Способ постановки граничных условий оказывает следующее влияние на погрешность. Постановка граничных условий в напряжениях приводит к появлению в решении пограничной зоны, в которой вследствие краевых эффектов
напряженное состояние быстро изменяется с удалением от границы и поэтому не соответствует реальному. При жестком закреплении границ пограничная зона отсутствует, и решение имеет одинаковую погрешность как в середине расчетной области, так и у ее границ; однако отсутствует возможность учета тектонических воздействий. Наконец, при запрете депланации границ каждая грань модели сохраняет одну степень свободы смещения по нормали, что позволяет учесть тектонические воздействия в виде интегральной силы, действующей на сечение; недостаток такого решения состоит в том, что резко возрастают затраты вычислительного времени вследствие увеличения ширины ленты (фронта) матрицы системы уравнений. По-видимому, наилучшим способом учета тектонических воздействий можно считать их задание в виде начальных напряжений сто, действующих в недеформированном массиве. Тогда в функционал (1) следует добавить работу этих напряжений:
П = — | [ст(є) +ст0 ] • єdV - | и ■pg • dV -1 и • q • dV (2)
2 V V 5
Такое решение краевой задачи позволяет находить поле напряжений в сложно построенных массивах горных пород. По крайней мере, вблизи дневной поверхности получаемые результаты хорошо согласуются с данными геоморфологических наблюдений [З].
2. Моделирование гравитационных аномалий
Конечно-элементная модель расчетной области, построенная для решения краевой задачи, всегда представлена в виде большого по объему набора данных о геометрии внутренних неоднородностей, физико-механических констант слагающих пород, рельефе и т.п. Это создает предпосылки для использования однажды подготовленных данных для решения задач различной физической природы. С другой стороны, такие задачи взаимосвязаны по данным: результаты решения одной из них могут представлять исходные данные для решения другой. С этой точки зрения представляется целесообразным объединение программ решения таких задач в единый программный комплекс. Тогда расчет гравитационных аномалий может быть проведен с использованием базы данных, созданной для проведения расчета деформирования.
Для расчета гравитационной аномалии используем следующую расчетную модель. Расчетная область, имеющая форму прямоугольного параллелепипеда и разбитая на конечные элементы, окружена однородными плитами большого размера, также имеющими форму прямоугольных параллелепипедов, причем их плотность постоянна и одинакова и равна средней плотности земной коры. Аномалия напряженности поля сил тяжести (ускорения свободного падения) в любой точке на малой высоте над дневной поверхности находится как сумма составляющих, обусловленных вкладом каждой из плит и каждого из конечных элементов расчетной области [4]:
* , (3)
V Я
где у - гравитационная постоянная; р - разность фактической и средней плотности; 7 и 20 -координаты текущей точки и точки наблюдения; Я - расстояние от текущей точки до точки наблюдения.
Составляющие, соответствующие плитам, находятся аналитическим интегрированием по известным формулам.
Вклад расчетной области вычисляется как сумма интегралов по объему каждого конечного элемента, причем для вычисления этих интегралов также используется аналитическое интегрирование. Получаемые при этом результаты иллюстрируются на рис. 2.
Значения ускорения свободного падения, найденные расчетом, принципиально могут быть в дальнейшем использованы для расчета деформирования при гравитационных воздействиях, но обусловленное этим уточнение значительно меньше погрешности решения краевой задачи упругости (и тем более - пластичности), поэтому такое решение едва ли целесообразно. Большее значение имеет их использование для интерпретации данных полевых измерений гравитационных аномалий.
3. Краевая задача электропроводности и приложения к прямой задаче электроразведки постоянным током
Задача электропроводности неоднородного массива математически формулируется как краевая задача для эллиптического уравнения второго порядка в частных производных [5]:
Lu=div(-аVu)=f, (4)
и\г = и *’ ЛІГ = Л’ (ки + Іп 1 = О,
11 1 1 1 2 1 Г3
где f - плотность источника тока, а - тензор электропроводности, и -потенциал; /*п - вектор плотности тока, п -вектор единичной нормали, к - коэффициент пропорциональности между потенциалом на границе и плотностью тока через границу области, причем и *, /п и к - известные величины. Во всей области предполагается непрерывность потенциала.
Как показано В.А. Шеметовым [6], решение и краевой задачи (2), удовлетворяющее главным граничным условиям на Г—, можно искать в пространстве L2(Q) как удовлетворяющее интегральному тождеству (слабое решение)
|(иЬф)йУ - |ф/П<^ + |фкЫЗ + ^ |и(и,Уф)п№ = |/<р<3¥
□ Г2 Г3 і дО, О
(5)
для всех р Є W2 (□) . В отличие от вариационной постановки, это позволяет моделировать точечные источники тока, в том числе — находящиеся на границах раздела анизотропных сред. Представив искомое поле потенциала в виде сплайна, из уравнения (5) получим систему линейных алгебраических уравнений высокого порядка относительно значений потенциала в узлах конечно-элементной модели. Эта система имеет симметричную положительно определенную матрицу, и ее решение не представляет непреодолимых трудностей. Особо отметим, что коэффициенты системы разрешающих уравнений вычисляются, исходя из тех же данных, что и при решении задачи о деформировании массива; дополнительно требуется задать лишь электрическую проводимость слагающих пород. Таким образом, и в этих задачах имеется значительный выигрыш в технологии исследования, если расчеты ведутся в рамках одного комплекса программ на одной базе данных.
Как и для задач механики, ограниченность размеров расчетной области является возможным источником погрешности и требует непременного рассмотрения при решении прикладных задач. В равенствах (4)-(5) граничные слагаемые соответствуют трем видам граничных
условий на различных участках границы: заданию граничного потенциала (гранич-ное условие первого рода), заданию плотности тока через границу (второго рода) и заданию коэффициента пропорциональности между потенциалом на границе и плотностью тока через границу (условие третьего рода, аналогичное конвективному граничному условию в задачах теплопроводности). Исходя из результатов численных экспериментов, можно утверждать, что на искусственно введенной границе расчетной области следует задавать граничное условие третьего рода, в то время как на дневной поверхности - второго (условие непротекания тока). Граничные условия первого рода являются излишними.
Возможность учета точечных источников тока позволяет использовать решение краевой задачи для моделирования электроразведки постоянным током и переменным током на низкой частоте, когда волновыми эффектами можно пренебречь. Построены методики и алгоритмы расчета для дипольного электропрофилирования методом сопротивлений, метода мелкомасштабного заряда и метода заряда с измерением магнитного поля. В последнем случае после отыскания поля потенциала рассчитывается поле плотностей электрического тока и затем численным интегрированием находится поле вектора напряженности магнитного поля. Заметим, что, несмотря на такой опосредованный способ вычисления магнитного поля, результаты расчета удовлетворительно согласуются с точными решениями модельных задач (рис. 3).
4. Эвристические приемы повышения точности
Отметим, что точность решения краевой задачи сеточными методами, в том числе - методом конечных элементов, не вполне отвечает требованиям практики. Поскольку расчетная область велика по размерам и содержит большое количество неоднородных включений, размеры конечных элементов оказываются недостаточно малыми, и решение чувствительно к влиянию границ, а также к изменению формы и числа элементов. Поэтому необходимы дополнительные меры по уменьшению влияния перечисленных факторов. С этой целью В.П.
Рис. 2. Аномалия силы тяжести, вызванная рельефом: а - на возвышенности, б - на впадине
Рис. 3. Профили напряженности магнитного поля: _____- теория,______- использование граничных условий первого рода, - - - - использование граничных условий третьего рода
Ластовецким [7] разработаны некоторые эвристические приемы, строгое обоснование которых пока не найдено, но выигрыш от использования наглядно подтверждается численными экспериментами.
Упомянутые приемы основаны на использовании аналитических решений модельных задач. Основная идея состоит в том, что вначале строится конечноэлементная модель неоднородной расчетной области, а затем - вспомогательная модель той же задачи на однородной области, для которой точное аналитическое решение, как правило, не вызывает заметных трудностей. Влияние конечности расчетной области и конечности размеров элементов проявляется в том, что найденное численное решение вспомогательной задачи отличается от точного. Затем в численное решение основной задачи вводится поправка на фактически полученную погрешность вспомогательной задачи.
Для задач электроразведки такая методика в общих чертах сводится к отысканию нормального поля потенциала (или кажущегося сопротивления) и умножению численного решения основной задачи в каждой расчетной точке на коэффициент, равный отношению точного значения в этой точке к численно найденной величине нормального поля.
Для задач об упругом деформировании неоднородного массива точное решение известно в случае однородного массива с плоским рельефом. Поэтому дополнительная трудность состоит здесь в учете влияния рельефа. Такой учет может быть проведен методом малого параметра, но решение получается в рядах, а не в замкнутом виде.
5. Модели с переменным параметром и задачи идентификации
Рассмотренные выше краевые задачи при использовании сеточного метода решения сводятся к системе линейных алгебраических уравнений. Поскольку эта система имеет высокий порядок (до сотен тысяч уравнений), а численный метод принципиально дает решение лишь при известных числовых значениях всех входящих
в модель параметров, такие математические модели связаны с проведением длительных и трудоемких расчетов. Но с точки зрения их прикладного использования важна возможность оперативного анализа влияния изменения параметров на решение задачи. Поэтому целесообразно использование комбинированных численно-аналитических моделей [8], которые позволяют получить решение в виде аналитически заданной функции переменного параметра с числовыми коэффициентами. Имея набор таких решений, можно достаточно экономично оценивать параметры задачи, которые известны лишь приближенно, по данным полевых измерений, т.е. в некотором приближении решать задачи идентификации свойств массива горных пород.
Рассмотрим математическую модель в виде, обобщающем соотношения (1) и (5):
|(uL(X)ф)dV - ||фк(X)ис!8 + ^ |и(стгУф)п№ = |fфdV
О Г Г3 г дО* О
(6)
где X - варьируемый параметр. Им может быть модуль упругости породы либо ее коэффициент электропроводности. Путем дискретизации (5) получаем систему уравнений с переменным параметром:
К(Х)и = R, (7)
где К(Х) - матрица системы, зависящая от варьируемого параметра.
Решение системы уравнений (6) может быть получено в виде разложения по степенным и дробнорациональным функциям [8], причем алгоритмы решения оказываются существенно более экономичными, чем последовательный расчет многих вариантов с перебором значений параметра. Существенно то, что такой подход позволяет оценить, насколько неточность исходных данных влияет на изменение решения краевой задачи, поэтому уменьшается возможность неверной интерпретации данных полевых измерений.
Заключение
Математическое моделирование полей различной физической природы в сложно построенных массивах горных пород на базе общей методологии и единого пакета прикладных программ позволяет многократно использовать однажды подготовленные исходные данные для решения задач разного физического содержания. Представляется, что такой подход повышает качество интерпретации данных полевых измерений и делает возможной оперативную интерпретацию наблюдений в задачах разведочной геофизики. Как показал имеющийся опыт эксплуатации разработанного пакета программ, обработка данных может проводиться непосредственно в геофизических экспедициях, независимо от сложности строения массива.
СПИСОК ЛИТЕРАТУРЫ
1. Зенкевич О. Метод конечных элемен-
тов в технике. - М.: Мир, 1975.-541 с.
2. Оден Дж. Конечные элементы в не-
линейной механике сплошных сред. - М.: Мир, 1976. - 464 с.
3. Каледин В.О., Ластовецкий В.П. Математическое моделирование напряженно-деформированного состояния горных пород применительно к нефтегазопоисковым задачам. - Геофизика, №3, 1999, с. 63-68.
4. Гравиразведка /ред.: Мудрецова Е.А., Веселова К.Е. - М.: Недра, 1990.-607 с.
5. Якубовский Ю.В., Ляхов Л.Л. Электроразведка. - М.: Недра, 1988.-395 с.
6. Каледин В. О., Ластовецкий В.П.,
Шеметов ВА. Теория и практика решения задач электроразведки постоянным током с использованием пер-сонального компьютера. -Техника и технология разработки месторождений полезных ископаемых /между нар. науч.-технич. сборник/ Вып. 4. Новокузнецк: Сибирский государственный индустриальный университет, 1998. - С.60-68.
7. Ластовецкий В.П. Учет физических особенностей полей при их численном моделировании для уточнения результатов расчета // Численно-аналитические методы
решения краевых задач. Сборник трудов 2й межвузовской конференции. - Новокузнецк: 1999. - С. 85-88.
8. Протасов В.Д., Каледин В.О., Каледин В.О., Суханов А.В. О численно - аналитическом исследовании термоупругого деформирования многосвязных анизотропных конструкций размеро-стабильных платформ. - Доклады Академии наук, 1996, т. 346, № 6,с.757-761.
Каледин Валерий Олегович - профессор, доктор технических наук, Новокузнецкий филиал-институт Кемеровского государственного университета.
