CC BY
30
© Л.А. Назаров, Л.А. Назарова, 2013
Л.А. Назаров, Л.А. Назарова
ОБРАТНЫЕ ЗАДАЧИ ГЕОМЕХАНИКИ: ОПРЕДЕЛЕНИЕ СВОЙСТВ И НАПРЯЖЕНИЙ В МАССИВЕ ГОРНЫХ ПОРОД*
Для антропогенных и природных объектов рассмотрены обратные задачи, решение которых позволяет оценить естественные горизонтальные напряжения, найти константы уравнений состояния, описывающих процесс вязкоупругого деформирования породных массивов, а также определить фокальные параметры предстоящего сейсмического события.
Ключевые слова: массив горных пород, обратные задачи, вязко-упругое деформирование, природные напряжения.
Знание напряжений и свойств массивов горных пород необходимо для геомеханического обоснования технологий разработки месторождений полезных ископаемых, интерпретации данных мониторинга состояния приконтурных зон, прогнозных расчетов длительной прочности конструктивных элементов подземных сооружений, а также моделирования техногенных геодинамических процессов. Практически универсальным подходом к решению этой проблемы является формулировка и решение обратных задач. Построенная соответствующим образом целевая функция позволяет не только включить в анализ комплекс разнородной информации (например, объем выделившегося газа, температуру призабойной зоны и вызванное подработкой оседание дневной поверхности), но и варьированием весовых коэффициентов устанавливать приоритеты типов этой информации в зависимости от ее значимости и качества.
В отличие от геофизики, объем входных данных для обратных задач геомеханики сравнительно невелик. Это обуслов-
* Работа выполнена при частичной поддержке Российского Фонда Фундаментальных Исследований (грант № 12—05—00843), Интеграционных проектов СО РАН № 76 и № 100.
лено высокой стоимостью проведения натурных экспериментов, существенной «дискретностью» расположения пунктов измерения и, как правило, квазистатическим характером процессов деформирования и разрушения горных пород. Данные обстоятельства, по-видимому, объясняют сравнительно редкое использование обратных задач в классических постановках [1, 2, 3] для решения практических задач механики горных пород. Только в последнее десятилетие благодаря широкому распространению GPS технологий мониторинга земной поверхности [4], а также разработке новых методов космической геодезии для детального определения поверхностных деформаций (радарная интерферометрия [5]), появилась возможность получать большие массивы данных, интерпретация которых может дать информацию о геомеханических и геодинамических [6] процессах в верхней части земной коры.
Рассмотрим несколько примеров решения граничных и коэффициентных обратных задач для геомеханических объектов различного масштабного уровня. Все расчеты проводились методом конечных элементов по оригинальным кодам, реализующим вязкоупругие модели деформирования сред в различных системах координат в 2D и 3D постановках [7, 8]. Поиск минимумов целевых функций осуществлялся методом сопряженных градиентов [9] в модификации [10].
1. Обратные граничные задачи.
1.1. Реконструкция краевых условий по информации о геомеханических полях
во внутренних точках области.
Компоненты упругих полей sij (напряжений, деформаций) в некоторой области А можно найти [9] из соотношения
s9 (a) = {у(Ь^ (p -b)dB , (1)
в
где a e A , b e B = dA ; Sij - соответствующая функции Грина; функция у описывает распределение смещений и/или напряжений на B . Введем в А сетку с узлами am , тогда дискретный аналог (1) имеет вид s.. (a ) = V у S.. ,
. v m / / , Т n ljmn ' neNB
где V = У(Ь„), = Бц (ат - Ьп )АВп, АВп — элемент поверхности В , МВ -множество номеров узлов, принадлежащих В . Пусть в узловых точках области С с А известны функции (ат) (1 = 1,..., L), зависящие от параметров геомеханических
полей, которые можно измерить в натурных условиях. Введем целевую функцию
ь 2
ф(¥п ) =!Т1 I [1 (V п , ат ) - ^т )] , (2)
1 =1 теМС
где у1 +у2 +... + уь = 1, у1 е [0,1] — весовые коэффициенты, а функции соответствуют ', но с произвольными значениями уп. Точки минимума Ф доставляют решение обратной граничной задачи определения неизвестных краевых условий по информации во внутренних точках области.
Если некоторые из функций 11 линейные, то процедуру
поиска минимума ф можно существенно ускорить, поскольку часть производных Фу вычисляются аналитически. В частности, при 11 = ви она сводится к решению системы линейных уравнений
3 3
I Уп I I Б
Итп Б Итк II БИтк^И (ат ) , к е ^В .
пеМВ 1 =1 теМС 1 =1 теМС
Проиллюстрируем описанный подход на примере крупномасштабного геологического объекта — Центральной Азии и ее обрамления. В [12] разработана объемная геомеханическая модель (рис. 1, дневная поверхность), учитывающая основные структурные элементы (плиты и микроплиты), тектонические разломы и соответствующие типы геодинамических режимов, вертикальную стратификацию физических свойств в каждом блоке.
Анализ сейсмического режима, при котором учитываются местоположение гипоцентров землетрясений, амплитуда и тип подвижки берегов разрыва, магнитуда М позволяет установить параметры поля сейсмотектонических деформаций (ориентацию главных деформаций, а также их знак и относительную
80
40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 ВД
_ тектонические \ область задания входных данных
разломы I------- при решении обратной задачи
Рис. 1. Геомеханическая модель Центральной Азии: блоковая структура, сетка конечных элементов и граничные условия
величину). Такое исследование с использованием методики [13] проведено (О.А.Кучай, ИНГГ СО РАН) для центральной части рассматриваемого объекта (штриховой прямоугольник на рис. 1). На рис. 2 кругами различного диаметра показаны распределенные по конечным элементам сейсмодеформации, а также величина и ориентация главных горизонтальных деформаций в1 и в2 там, где суммарная энергия всех событий (за
рассмотренный период времени около 40 лет) превышает таковую для землетрясений с М > 6.5. Эти данные послужили входной информацией для оценки смешений Ц3 на границах
периферийных плит (рис. 1), вызвавших такое распределение деформаций.
В качестве функций 1 (L = 1) выберем удвоенный угол между первой главной деформацией и направлением на восток
11 (ат) = агф 1 т ее т ,
Ч (ат )
60 66 70 75 80 85 90 95 100 105 110 115 ВД
Рис. 2. Сейсмотектонические деформации — входные данные для решения обратной задачи
где ат = (9т, фт) — географические координаты (широта и
долгота) центра элемента (поскольку глубина гипоцентров определяется, как правило, весьма приближенно, то все данные осреднялись по верхнему слою земной коры мощностью 30 км, в котором находится подавляющее большинство очагов землетрясений).
Неизвестными аргументами уп целевой функции Ф являлись модули ип векторов граничных смещений ип (рис. 1): их ориентация считалась известной и совпадающей с направлением (задавалось азимутом ап) движения литосферных плит [14].
В результате минимизации Ф были найдены значения и п
(табл. 1), при которых величина и ориентация модельных деформаций соответствовала реальным (рис. 3, на врезке показан укрупненный фрагмент для сравнения расчетных и натурных данных).
Верифицированная таким образом геомеханическая модель Центральной Азии позволяет при сохранении тенденции
Таблица 1
Значения граничных смешений
№ платформы п Модуль ип мм Азимут ап ,град
диапазон итог
Индо-Австралийская, запад 1 55—75 63 45
Тихоокеанская 2 80—100 95 315
Индо-Австралийская, восток 3 - 50 5
Европа, север 4 20—40 25 90
Европа, юг 5 - 5 90
Рис. 3. Распределение деформаций, полученное в результате решения обратной задачи
движения тектонических плит выполнить прогнозную оценку параметров процесса эволюции геодинамических полей и сделать выводы о предстоящем изменении сейсмического режима территории.
1.2. Определение параметров предстоящего сейсмического события по деформациям дневной поверхности.
Любому динамическому событию в массиве горных пород предшествует стадия квазистатического деформирования [15],
причем в окрестности будушего очага (например, участка тектонического разлома Я с высокой жесткостью, рис. 4) возникает область концентрации напряжений и деформаций, «след» которой проявляется на дневной поверхности в виде аномальной зоны повышенных деформаций. Последние, особенно при небольших глубинах очагов, могут быть обнаружены и количественно оценены методам космической геодезии: на больших базах абсолютная точность измерения деформаций достигает 10-7 [4, 16].
Определить детальную структуру и свойства Я практически нереально, поэтому воспользуемся широко распространенным в сейсмологии приемом [17]: построим эквивалентный точечный источник, генерируюший на поверхности поле деформаций г.., близкое к таковому от аномальной зоны в*.. Не вдаваясь в подробности, которое можно найти в [6], конкретизируем структуру целевой функции (2). Здесь в качестве искомых параметров уп выступают декартовы компоненты Рх , ¥у и ¥х
двойной силы с моментом (типичной модели очага землетрясения [17]), а также координаты точки ее приложения х0, у0и
ок и Он- минимальное и максимальное горизонтальные напряжения в нетронутом массиве
Рис. 4. Схема расчетной области и модель эквивалентного источника 50
Рис. 5. Изолинии касательных деформаций, индуцированные на дневной поверхности аномальной зоной и эквивалентным источником
z0 (рис. 4). Весовые коэффициенты принимали значения
у1 = 1, у 2 = у 3 = 0, т.е. для синтеза источника использовалась
только компонента деформаций sxx. На рис. 5 в эпицентраль-
ной зоне (штриховой прямоугольник на рис. 4) показаны изолинии касательных деформаций, индуцированных аномальной зоной (exy ) и эквивалентным источником (exy ): несмотря на то,
что эта компонента не фигурировала в функции Ф, относительная ошибка не превышает 10 %.
Найденные таким образом параметры эквивалентного источника позволяют на основе эмпирических соотношений [18] оценить магнитуду готовящегося сейсмического события. 2. Обратные коэффициентные задачи. 2.1. Определение упругих характеристик конструктивных элементов технологий выемки запасов по смещениям контура выработанного пространства
Многие рудные месторождения с пологим залеганием пластов разрабатываются сплошной слоевой системой с закладкой выработанного пространства или оставлением целиков различного назначения [19, 20]. Деформационные и реологические характеристики несущих элементов влияют на технологию отработки, поэтому их достоверное определение in situ позволяет повысить эффективность горных работ.
W = O O^ = о (J¡. (y) =pgy - литостатическое напряжение Ц - коэффициент бокового отпора Рис. 6. Структура расчетной области, схема измерения и граничные условия
Продемонстрируем способ оценки модуля Юнга Е и коэффициент Пуассона v закладочного массива по данным измерения относительных смещений контура выработок. На рис. 6 изображен фрагмент типичной конфигурации подземного пространства при отработке рудного месторождения с закладкой. Очистные работы моделировались увеличением протяженности камеры Ln: на каждом этапе n забой продвигается на Ах = 0.125h (h -мощность отрабатываемого пласта) в горизонтальном направлении. Поведение массива описывается уравнениями линейной теории упругости, выполнены условия плоской деформации.
Пусть в камере C регистрируется изменение расстояния Aun и Awn между точками контура C , вызванное продвижением забоя (рис. 6). Введем целевую функцию
N 2 N 2
Ф(Е, v) = Yi ^[l -Aun (E, v)/Aun ] +Y2 Z[l "AWn (E, v)/Awn ]
n=1 n=1
где y1 + y2 = 1, а число N характеризует максимальную протяженность C , при которой закладочный массив еще находится в упругом состоянии. Точка (Е., v.) глобального минимума
Ф и является решением обратной коэффициентной задачи определения упругих характеристик закладки.
В качестве входных использовались синтетические данные Áu'n = (1 + Áun (E,,v,), = (1 + q2)Awn (E,,v,),
где Aun и Awn - решение прямой задачи при искомых значениях модуля Юнга E, и коэффициента Пуассона v,, а равномерно распределенные на отрезке [-Е, Е] случайные величины и q2 моделировали ошибку во входных данных, Е - амплитуда помехи.
На рис. 7 показаны изолинии Ф для N = 10 , у1 = 0.25 и у2 = 0.75 при различных значениях Е . Как и следовало ожидать, с увеличением Е у целевой функции появляются несколько минимумов, среди которых необходимо выбирать самый «глубокий», что замедляет процедуру поиска решения. Тем не менее, даже при заметной помехе коэффициент Пуассона v, находится практически точно, а модуль Юнга E, с ошибкой 10-15 %.
2.2. Оценка реологических свойств.
Пусть породы отрабатываемого пласта или закладка проявляют реологические свойства. Для описания вязкоупругого деформирования таких материалов используем следующие уравнения состояния [21] t
Ke(t) = c(t) + о[(t -т)-ус(x)dх ,
0
t
цё. (t) = с.. (t) + pj(t - т)-УС.. (x)d т , (3)
0
где K и ц — модули объемного сжатия и сдвига, с — среднее напряжение, e — объемная деформация, t — время, а и Р - реологические параметры (размерность с1-7), у — показатель степени в ядре Абеля, черта означает девиатор соответствующего тензора. В качестве входной информации для определения а , Р и у будем использовать изменение во времени конвергенции кровли и почвы Aw'N(t), зарегистрированное при неподвижном забое. При этом по описанной в п. 2.1 методике деформационные характеристики E и v (по которым можно рассчитать K и ц) уже установлены.
Рис. 7. Изолинии целевой функции Ф при различных значениях параметров
Рис. 8. Изолинии целевой функции и схема поиска глобального минимума
Введем целевую функцию
т 2
Ф(а, р, у) = |[1 -AwN (а, р, у,()/ Aw"N( т) ] 6 т ,
где и Т — моменты начала и окончания измерений (значение Т во многом зависит от чувствительности датчика смеще-54
ний), а Aw"n(t) = (1 + g2)Awn (а,,р., у,^). Численные эксперименты показали, что величина параметра а практически не влияет на структуру Ф, поэтому на рис. 8 показаны изолинии целевой функции в сечении а = 0.0003 (среднее значение для соляных пород [21]) при Е, = 0.1.
Можно видеть, что функция Ф имеет несколько локальных минимумов, расположенных на одной прямой L . Поэтому для решения обратной задачи достаточно отыскать два любых минимума М1 и М 2, а затем найти наименьшее значение Ф на
прямой М1М 2.
Описанный метод интерпретации результатов измерений смещений точек контура выработанного пространства позволяет, например, в процессе выемки запасов осуществлять мониторинг изменения реологических свойств закладочного массива и на этой основе проводить прогнозные расчеты устойчивости последнего.
- СПИСОК ЛИТЕРАТУРЫ
1. Тихонов А.Н., Арсенин В.Я. Методы решения обратных задач. - М.: Наука, 1974. - 224 с.
2. Tarantola A. Inverse Problems. Theory and Methods. Siam, 2005. 358 p.
3. Романов В. Г. Обратные задачи математической физики. - М.: Наука, 1984. 264 с.
4. Teunissen P.J.G., Kleusberg A. GPS for Geodesy. Springer, 1998. 650
p.
5. Akcin H., Degucci T., Kutoglu H.S. Monitoring Mining Induced Subsidence Using GPS and InSAR. Proc. XXIII FIG Congress, Munich, Germany, 2006, October 8—13, pp. 1—12.
6. Назарова Ё.А., Назаров Ё.А. Метод определения параметров очага готовящегося землетрясения на основе данных о смещениях дневной поверхности // Доклады Академии Наук. 2009. Т. 427. № 4. C.534—538.
7. Назарова Ё. А. Напряженное состояние наклонно-слоистого массива горных пород вокруг выработки // ФТПРПИ. 1985. № 2. C.33—37.
8. Назарова Ё.А. Моделирование объемных полей напряжений в раз-ломных зонах земной коры // Доклады Академии Наук. 1995. Т.342. № 6. C.804—808.
9. Васильев Ф.П. Численные методы решения экстремальных задач. М.: Наука, 1988, 550 с.
10. Карчевский А.Ё. Численное решение одномерной обратной задачи для системы упругости // Доклады Академии Наук. 2000. Т. 375. № 2. С. 235—238.
11. Новацкий В. Теория упругости. - М.: Мир, 1975. - 872 с.
12. Ёеви К.Г., Назаров Ё.А., Назарова Ё.А. и др. Актуальные вопросы современной геодинамики Центральной Азии. - Новосибирск: Издательство СО РАН, 2005. -445 с.
13. Юнга С.Ё. Методы и результаты изучения сейсмотектонических деформаций. - М.: Наука, 1990. - 191 с.
14. http ://www.unavco.org/community_science/science-support/crustal_motion
15. Добровольский И.П. Теория подготовки тектонического землетрясения. - М.: Наука, 1991. - 219 с.
16. Hofmann-Wellenhof B., Lichtenegger H., and Collins J. GPS Theory and practice. Springer, New York, 2001. 383 p.
17. Аки К., Ричардс П. Количественная сейсмология. Т. 1. - М.: Мир, 1983. 519 с.
18. Ризниченко Ю.В. Проблемы сейсмологии. - М.: Наука, 1985. 408 с.
19. Бондаренко В.И., Кузьменко А.М., Грядущий Ю.Б. и др. Технология подземной разработки пластовых месторождений полезных ископаемых. Днепропетровск, 2002. 206 с.
20. Бронников Д.М., Замесов Н.Ф., Богданов Г.И. Разработка руд на больших глубинах. - М.: Недра, 1982. - 292 с.
21. Барях А.А., Константинова С.А., Асанов В.А. Деформирование соляных пород. - Екатеринбург: УрО РАН, 1996. - 204 с. ii^J-i
КОРОТКО ОБ АВТОРАХ
Назаров Леонид Анатольевич, доктор физико-математических наук, заведующий лабораторией, naz@misd.nsc.ru,
Назарова Лариса Алексеевна — доктор физико-математических наук, заведующий отделом, larisa@misd.nsc.ru,
Институт горного дела им. Н.А. Чинакала Сибирского отделения РАН.
