Building/chapters/IFS.typ

#import "/.template/lib/index.typ": corp-table, formula, vref, vrefs, eqref, numbered-list

#let img-width = 14cm
#import "@preview/rexllent:0.4.1": xlsx-parser

= Определение природного поля напряжений на участке Кварцитовый

== Введение

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

Расчёт выполнен в три этапа:
#numbered-list(scheme: "decimal")[
+ По механизмам очагов землетрясений определены направление максимального горизонтального напряжения и коэффициент формы тензора напряжений $R$.
+ По величине вертикального напряжения, поровому давлению, исходным свойствам тектонических разломов, полученным в ходе инженерно-геологических изысканий @geomechKvarcitovyi2026 и моделирования, рассчитаны величины максимального и минимального горизонтального напряжения.
+ Полученные значения природного поля напряжений переведены в систему координат моделей разрезов и заданы в ПК Prorock.
]

В главе используются следующие основные термины.

- *Событие* -- зафиксированное землетрясение, рассматриваемое как элемент выборки данных (термин сейсмологических каталогов).
- *Механизм очага* -- описание разрыва в очаге землетрясения как смещения двух блоков пород по плоскости. Механизм задаётся ориентацией двух узловых плоскостей и направлением смещения по каждой из них.
- *Узловые плоскости* -- две взаимно перпендикулярные плоскости, проходящие через очаг и разделяющие области сжатия и растяжения первых колебаний продольных волн. Одна из них совпадает с плоскостью реального разрыва, вторая называется вспомогательной.
- *Невязка* (misfit angle, $beta$) -- угол между наблюдаемым вектором смещения в очаге и рассчитанным направлением максимального касательного напряжения.
- *Обратный расчёт напряжений* -- подбор по совокупности механизмов очагов единого тензора напряжений, наилучшим образом объясняющего наблюдаемые направления смещений в очагах.
- *Коэффициент формы тензора* $R$ -- коэффициент, показывающий относительное положение промежуточного главного напряжения между наибольшим и наименьшим.
- *Сдвиговый режим напряжений* -- режим, при котором максимальное и минимальное главные напряжения горизонтальны, а вертикальное занимает промежуточное положение.

== Сбор исходных данных

При помощи сервиса Google Earth определены точные координаты искомых объектов. Маломырское месторождение находится в координатах 53°03′34,92″ с. ш. и 131°42′28,85″ в. д. Участок Кварцитовый расположен на правом борту ручья Маломыр -- левого притока реки Нижняя Стойба.

Для определения ориентации и величины природного поля напряжений использована база данных World Stress Map (WSM, 2025 г.) @heidbachWorldStressMap2025. Ближайшая к участку точка наблюдения имеет код wsm144823 и расположена в координатах 51,53° с. ш., 131,95° в. д., на расстоянии $approx 172$ км от участка Кварцитовый. 

Расположение участка Кварцитовый и точки наблюдения показано на #vref(<pointsmap>, "п").

#figure(
  image("/assets/images/points_on_map.png", width: 70%),
  caption: [Расположение точки wsm144823 относительно участка Кварцитовый на WSM.],
  gap: 10pt,
)<pointsmap>

Точка wsm144823 соответствует землетрясению 22.07.2013 г. на юго-востоке Сибири
с моментной магнитудой $M_w=4","8$ на глубине 24 км. Исходный механизм очага
опубликован в каталоге #link("https://www.globalcmt.org/")[Global CMT]
@ekstromGlobalCmtProject2012 с идентификатором `C201307221508A` со следующими
данными (#vref(<wsm144823>, "и").)


#figure(
  image("/assets/images/wsm144823.png", width: img-width),
  caption: [Данные из точки wsm144823 с ресурса WSM],
  gap: 10pt,
)<wsm144823>

Обе плоскости разлома, приведённые в записи, характеризуются сдвиговым типом смещения, что соответствует сдвиговому режиму напряжений, представленному на #vref(<strike-slip_faulting>, "п").

#figure(
  image("/assets/images/strike-slip_faulting.png", width: 30%),
  caption: [Представление сдвигового режима.],
  gap: 10pt,
)<strike-slip_faulting>

Для выполнения обратного расчёта напряжений база WSM предъявляет ко входным данным численные критерии качества: для данных класса B необходимо не менее 8 независимых механизмов очагов и средняя невязка $beta$ не более $20degree$ @heidbachWorldStressMap2025. Схема определения угла невязки показана на #vref(<misfit-angle>, "п").
Рейтинговая система оценки данных WSM приведена на #vref(<ranking-scheme>, "п").

#figure(
  image("/assets/images/misfit-angle.png", width: 40%),
  caption: [ Схема определения угла невязки (misfit angle, $beta$) между наблюдаемым направлением скольжения и направлением максимального касательного напряжения на плоскости разлома.],
  gap: 10pt,
)<misfit-angle>

== Методика обратного расчёта природного поля напряжений

Механизмы очагов описывают разрыв в источнике землетрясения как смещение двух блоков по плоскости и задаются ориентацией и направлением смещения. Механизмы восстанавливают по записям сейсмических волн на станциях, расположенных в разных направлениях от очага. Классический подход опирается на знак значения колебания продольной волны: в одних направлениях от очага первый импульс регистрируется как движение от очага (сжатие), в других -- к очагу (растяжение). Нанесённые на воображаемую сферу колебания разделяют её на четыре чередующихся сектора -- два сжатия и два растяжения. Границами секторов служат две взаимно перпендикулярные плоскости, проходящие через центр сферы, -- *узловые плоскости* (#vref(<mechanism>, "и")). 

#figure(
  image("/assets/images/mechanism.png", width: 95%),
  caption: [Первые колебания механизма очага землетрясения],
  gap: 10pt,
)<mechanism>

Современные каталоги, включая использованный здесь Global CMT, определяют те же параметры подбором полной волновой записи: ориентации двух узловых плоскостей и направление смещения по каждой из них. Такой результат принято изображать чёрно-белой диаграммой очага.

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


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

В обратном расчёте принимается гипотеза Уоллеса--Ботта @wallaceGeometryShearingStress1951, @bottMechanicsObliqueSlip1959. Согласно гипотезе, направление смещения по действительной плоскости разрыва совпадает с направлением максимального касательного напряжения, рассчитанного на этой плоскости для подбираемого тензора напряжений. Каждый механизм очага задаёт две возможные плоскости разрыва. Реализуется только одна из них, но по сейсмическим данным заранее неизвестно, какая именно. Для каждой из них вычисляют невязку β, в расчёт принимают плоскость с меньшей невязкой. 

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

Обратный расчёт выполнен по линейной схеме Майкла
@michaelDeterminationStressSlip1984. Метод опирается на два допущения: поле
напряжений однородно в пределах рассматриваемой области, и скольжение в очаге
происходит в направлении максимального касательного напряжения на плоскости
разрыва (гипотеза Уоллеса--Ботта @wallaceGeometryShearingStress1951,
@bottMechanicsObliqueSlip1959). <assumption> Касательное напряжение, действующее на
плоскость разрыва с единичной нормалью $arrow(n)$, равно:

#formula[$tau = sigma dot arrow(n) - (arrow(n) dot sigma dot arrow(n)) dot arrow(n)$,]

где $tau$ -- вектор касательного напряжения, действующего в плоскости разрыва; $sigma$ -- тензор напряжений; $arrow(n)$ -- единичный вектор нормали к плоскости разрыва. Первый член $sigma dot arrow(n)$ задаёт полный вектор напряжения на площадке с нормалью $arrow(n)$. Величина $arrow(n) dot sigma dot arrow(n)$ представляет собой нормальное напряжение на этой площадке, а второй член -- его векторную составляющую, направленную по нормали. Их разность $tau$ лежит в плоскости разрыва и определяет направление касательного напряжения (#vref(<stress-vector>, "и")).

#figure(
  image("/assets/images/stress-vector.png", width: 50%),
  caption: [Разложение вектора касательного напряжения на плоскости разрыва],
  gap: 10pt,
)<stress-vector>

Условие совпадения вектора $tau$ по направлению и величине с единичным вектором наблюдаемого смещения $arrow(s)$ линейно относительно компонент тензора $sigma$, поэтому искомый девиатор напряжений можно определить решением системы уравнений методом наименьших квадратов. Поскольку механизмы очагов сообщают направление смещения, но не абсолютную величину касательного напряжения, в линейной схеме предполагается, что на плоскостях разрывов, по которым произошли землетрясения выборки, эти величины приблизительно одинаковы. Поэтому расчёт определяет относительные, а не абсолютные напряжения. Для удобства модуль касательного напряжения принимают равным единице @barthNewConstraintsIntraplate2010.

Форма тензора характеризуется коэффициентом $R$
@riveraSpatialHeterogeneityTectonic2002:

#formula[$R = frac(sigma_2 - sigma_3, sigma_1 - sigma_3)$]
#formula[$quad 0 <= R <= 1$]
где $R$ -- коэффициент, показывающий относительное положение
промежуточного главного напряжения между наибольшим и наименьшим.

Обратный расчёт реализован специально разработанной программой на языке Python, её исходный код приведён в #ref(<R-calc>, supplement: [приложении]). 

Входными данными служат параметры механизмов очагов: простирание, угол падения и угол скольжения. Вычисления ведутся в географической системе координат, в которой каждый механизм по стандартным формулам преобразуется в единичные векторы нормали плоскости разрыва $arrow(n)$ и наблюдаемого смещения $arrow(s)$. 

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

Достоверность реализации проверена воспроизведением опубликованных обратных расчётов @barthNewConstraintsIntraplate2010 по данным таблиц 1 и 2 указанной публикации. Для Буреинского блока получен азимут $S_("Hmax") = 70","1degree$ (в публикации $70","2degree$) при $R = 0","40$ и средней невязке $12","3degree$ (в публикации $12","2degree$), контрольные невязки отдельных событий воспроизводятся с расхождением не более $0","6degree$. Для Станового района совпадают коэффициент формы ($R = 0","47$), средняя невязка ($15","1degree$, в публикации $15","0degree$) и наибольшая невязка очищенной выборки ($30","1degree$ при заявленном пределе $32degree$), однако азимут $S_("Hmax")$ составляет $63","6degree$ вместо опубликованного $59","1degree$. Расхождение $4","5degree$ лежит в пределах собственной неопределённости станового решения (доверительные области ориентации осей по бутстреп-анализу составляют около $plus.minus 20degree$) и на принятую в настоящем отчёте выборку не влияет: она основана на данных Буреинского блока, для которого обратный расчёт воспроизводится с расхождением $0","1degree$. Полный протокол сверки приведён в #ref(<R-calc>, supplement: [приложении]).

== Обоснование выборки данных о природном поле напряжений

Опубликованные данные обратных расчётов для территории Амурской плиты получены
в работе Barth & Wenzel @barthNewConstraintsIntraplate2010 по трём районам:
Хинган, Становой складчатый пояс на северной границе плиты и Буреинский блок
юго-восточнее него. Из них ближайшим к участку Кварцитовый является Становой.

Опубликованный для него обратный расчёт имеет рейтинг качества B по критериям WSM:
азимут $S_("Hmax") = 59","1degree$, $R = 0","47$, $N = 15$ механизмов при средней
невязке $bar(beta) = 15","0degree$. 

В ходе анализа данных из выборки исключено событие №47 (2 ноября 1973 г.) с невязкой $96degree$. Расхождения
такого масштаба авторы объясняют либо действием локальных сил, не связанных
с региональным полем, либо неопределённостью ориентации узловых плоскостей
механизма: по данным синтетических тестов @michaelSpatialVariationsStress1991,
стандартное отклонение ориентации плоскостей $15degree$ способно давать
среднюю невязку $bar(beta)$ около $43degree$. После исключения аномальных значений остальные не превышают $32degree$.

С обратным расчётом Станового района данные wsm144823 не согласуются. Невязка
этого события относительно тензора, восстановленного описанным выше расчётом
по опубликованной очищенной выборке ($N = 15$), составляет $61","3degree$. Если
включить событие в выборку и пересчитать тензор, невязка снижается до
$49","4degree$: такой вариант соответствует протоколу публикации, в котором
невязки потенциальных выбросов оцениваются по тензору, найденному до их
исключения. 

Оба значения превышают наибольшую невязку очищенной выборки
($32degree$). Для сравнения: из буреинской выборки исключено событие
с невязкой $42degree$, из становой -- событие с невязкой $96degree$.

По результатам @barthNewConstraintsIntraplate2010, доверительные области
ориентаций осей $sigma_1$ и $sigma_3$ при уровне доверия 95 %, построенные
бутстреп-анализом, составляют около $plus.minus 20degree$ по азимутам. О
решении для Станового района авторы пишут: «Хотя оптимальное решение даёт
ось $sigma_3$ с крутым погружением, отвечающим преимущественно взбросовому
режиму, широкие доверительные области в большей своей части накрывают
область сдвигового режима. Этот неоднозначный результат анализа
указывает на сильное влияние отдельных решений за механизм очага, которые
переключают режим либо на взбросовый, либо на сдвиговый. Дальнейшее разбиение
этого района невозможно из-за малого числа механизмов очагов»
@barthNewConstraintsIntraplate2010 (здесь и далее перевод авторов отчёта). В
обсуждении публикации указано: «Столь широкие доверительные области позволяют
интерпретировать поле напряжений как сдвиговое или преимущественно
взбросовое» @barthNewConstraintsIntraplate2010. Поэтому использовать
опубликованную становую инверсию в качестве базовой для участка со сдвиговым
режимом нежелательно.

Ближайшее событие согласуется с данными обратного расчёта Буреинского блока. В опубликованной
выборке также исключён один выброс -- событие 3 сентября 1999 г.
с невязкой $42degree$, поэтому $N = 11$. Невязка wsm144823 относительно
тензора Буреинского блока, восстановленного по опубликованной выборке,
составляет $19","6degree$ и не выходит за пределы невязок этой выборки. 

После добавления wsm144823 в обратный расчёт по $N = 12$ получаем:
$S_("Hmax") = 70","3degree$ (в публикации $70","2degree$), $R = 0","42$
(в публикации $0","40$) при $bar(beta) = 12","2degree$ и наибольшей невязке
$25","6degree$. Критерии WSM для данных класса B выполняются.


Поэтому в качестве выборки использованы данные обратного расчёта Буреинского блока и
механизм ближайшего к участку землетрясения wsm144823.

== Результаты обратного расчёта 

Результат обратного расчёта подтверждает принятый сдвиговый режим напряжений:
#numbered-list(scheme: "decimal")[
+ Максимальное главное напряжение $sigma_1$ соответствует максимальному горизонтальному напряжению $S_("Hmax")$.
+ Промежуточное напряжение $sigma_2$ соответствует вертикальному напряжению $S_v$.
+ Минимальное напряжение $sigma_3$ -- минимальному горизонтальному напряжению $S_("hmin")$.
]
Полные результаты обратного расчёта приведены в #vrefs((<inversion_results>, <inversion_orientation>), "п").
#block(breakable: false)[
#let inversion-results = csv("/assets/tables/inversion_results.csv", delimiter: ";")

#figure(
  table(
    columns: (1.5fr, 1fr),
    align: right + horizon,
    fill: (x, y) => if y == 0 { rgb("#FFD35F") },
    table.header(..inversion-results.first().map(strong)),
    // Ячейки CSV содержат Typst-разметку для математических обозначений.
    ..inversion-results.slice(1).flatten().map(s => eval(s, mode: "markup")),
  ),
  caption: [Результаты обратного расчёта],
) <inversion_results>
]

#block(breakable: false)[
#let inversion-orientation = csv("/assets/tables/inversion_orientation.csv", delimiter: ";")

#figure(
  table(
    columns: (1.5fr, 1fr, 1fr),
    align: right + horizon,
    fill: (x, y) => if y == 0 { rgb("#FFD35F") },
    table.header(..inversion-orientation.first().map(strong)),
    // Ячейки CSV содержат Typst-разметку для математических обозначений.
    ..inversion-orientation.slice(1).flatten().map(s => eval(s, mode: "markup")),
  ),
  caption: [Ориентация главных осей тензора напряжений],
) <inversion_orientation>
]

Расчёт включает 12 независимых механизмов очагов: 11 механизмов землетрясений
Буреинского блока по данным @barthNewConstraintsIntraplate2010 и механизм
ближайшего к участку землетрясения wsm144823. Полный список приведён
в #vref(<inversion_input_data>, "п").
#block(breakable: false)[
#let inversion-data = csv("/assets/tables/inversion_data.csv", delimiter: ";")

#figure(
  table(
    columns: (1.2fr, 1fr, 1.1fr, 1.4fr, 1.2fr),
    align: right + horizon,
    fill: (x, y) => if y == 0 { rgb("#FFD35F") },
    table.header(..inversion-data.first().map(strong)),
    ..inversion-data.slice(1).flatten(),
  ),
  caption: [Сводные данные механизмов землетрясений, использованные для обратного расчёта напряжений],
) <inversion_input_data>
]

Примечание: Strike (простирание) -- азимут горизонтальной линии на плоскости разлома, отсчитываемый от севера.
Dip (угол падения) -- наклон плоскости разлома относительно горизонтали.
Rake (угол скольжения) -- направление смещения пород по плоскости разлома, параметр характеризует преобладающий тип смещения: сдвиговый, сбросовый или взбросовый.

== Методика расчёта абсолютных значений природного поля напряжений

Вертикальное напряжение принято равным весу вышележащей толщи:

#formula[$S_v = rho dot g dot H$]

где $S_v$ -- вертикальная компонента природного поля напряжений, Па, $rho$ -- принятая средняя
плотность массива, кг/м³, $g$ -- ускорение свободного падения, м/с²,
$H$ -- вертикальная мощность перекрывающей толщи, м. 
Режим напряжений принят сдвиговым по данным точки наблюдения WSM. При сдвиговом
режиме главные напряжения соотносятся как

#formula[$S_(italic("Hmax")) > S_v > S_(italic("hmin")),$]

то есть максимальное горизонтальное напряжение $S_(italic("Hmax"))$ является
наибольшим из главных напряжений ($sigma_1$), вертикальное $S_v$ занимает
промежуточное положение ($sigma_2$), а минимальное горизонтальное
$S_(italic("hmin"))$ -- наименьшее ($sigma_3$).

Соотношение между горизонтальными напряжениями задаётся коэффициентом формы
тензора. Выражая $S_(italic("Hmax"))$ из определения $R$, получаем:

#formula[$R dot (S_(italic("Hmax")) - S_(italic("hmin"))) = S_v - S_(italic("hmin"))$]

#formula[$R dot S_(italic("Hmax")) - R dot S_(italic("hmin")) = S_v - S_(italic("hmin"))$]

#formula[$S_(italic("Hmax")) = frac(S_v - (1 - R) dot S_(italic("hmin")), R)$]

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

#formula[$|tau| <= mu dot sigma_n'$]

где $tau$ -- касательное напряжение на плоскости разлома, $sigma_n'$ --
эффективное нормальное напряжение на той же плоскости, $mu$ -- коэффициент
трения, принятый по данным отчёта о свойствах тектонических разломов
($mu = tan phi.alt$). Скольжение начинается раньше всего по плоскости,
составляющей с направлением $sigma_1$ угол $45degree - frac(phi.alt, 2)$:
это критически ориентированная плоскость, которой на круге Мора отвечает
точка касания с прямой Кулона. Условие устойчивости для неё, записанное
в главных эффективных напряжениях, принимает вид ограничения на отношение
крайних главных напряжений. При сдвиговом режиме $sigma_1 = S_(italic("Hmax"))$
и $sigma_3 = S_(italic("hmin"))$, поэтому

#formula[$frac(S_(italic("Hmax")) - P_p, S_(italic("hmin")) - P_p) <= A$]

#formula[$A = (sqrt(1 + mu^2) + mu)^2$]

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

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

#formula[$S_(italic("Hmax")) - P_p = A dot (S_(italic("hmin")) - P_p)$]

или, после выражения $S_(italic("Hmax"))$,

#formula[$S_(italic("Hmax")) = A dot S_(italic("hmin")) - (A - 1) dot P_p$]

Таким образом, искомые напряжения связываются системой двух уравнений:

#formula[
  $ 
    cases(
      S_(italic("Hmax")) = frac(S_v - (1 - R) dot S_(italic("hmin")), R),
      S_(italic("Hmax")) = A dot S_(italic("hmin")) - (A - 1) dot P_p,
    )
   $
]

== Результаты расчёта абсолютных значений природного поля напряжений

Вертикальная компонента природного поля напряжений составляет:

#formula[$S_v = rho dot g dot H = 2750 dot 9","81 dot 400 = 10791000 "Па" = 10","791 "МПа"$]

Коэффициент трения принят по углу трения 33,4° из результатов инженерно-геологических изысканий @geomechKvarcitovyi2026:

#formula[$mu = tan(33","4°) = 0","659$]

#formula[$A = (sqrt(1 + 0","659^2) + 0","659)^2 = 3","45$]

Поровое давление при глубине 400 м принято равным 3,250 МПа по результатам расчётных гидрогеологических моделей в ПО Rocscience Slide2. Заметного влияния на результаты точность принятого значения не оказывает: анализ чувствительности (#vref(<sensitivity_analysis>, "п")) показывает, что при изменении глубины от 200 до 600 м и порового давления от 1,750 до 4,875 МПа задаваемое в расчётных моделях отношение горизонтального напряжения к вертикальному изменяется лишь от 0,750 до 0,741.

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

#formula[
  $ cases(
    S_(italic("Hmax")) = frac(
      10791000 - (1 - 0","42) dot S_(italic("hmin")),
      0","42
    ),
    S_(italic("Hmax")) = 3","45 dot S_(italic("hmin"))
    - (3","45 - 1) dot 3250 dot 10^3,
  ) $
]

#formula[$S_(italic("Hmax")) = frac(10791000 - (1 - 0","42) dot S_(italic("hmin")), 0","42)$]

#formula[$S_(italic("Hmax")) = 3","45 dot S_(italic("hmin")) - (3","45 - 1) dot 3250 dot 10^3$]

#formula[$S_(italic("Hmax")) = 1","607 dot 10^7 "Па", quad S_(italic("hmin")) = 6966609","167 "Па"$]

#formula[$S_(italic("Hmax")) = 1","607 dot 10^7 "Па" = 16","07 "МПа"$]

#formula[$S_(italic("hmin")) = frac(6966609","167, 10^6) = 6","97 "МПа"$]

Принятое природное поле напряжений на глубине $H = 400$ м представляется
тензором в главных осях:

#formula[$bold(sigma) = mat(delim: "(", 16","07, 0, 0; 0, 10","79, 0; 0, 0, 6","97), quad "МПа"$]

Компоненты отнесены к главным осям тензора: $sigma_1 = S_(italic("Hmax"))$
действует горизонтально в направлении азимута $70","3degree$, $sigma_2 = S_v$ --
вертикально, $sigma_3 = S_(italic("hmin"))$ -- горизонтально в перпендикулярном
направлении ($159","9degree$).

== Перевод напряжений в систему координат расчётной модели

Компоненты $S_("Hmax")$ и $S_("hmin")$ заданы в географических осях, тогда как
расчётная модель строится в плоскости вертикального разреза. Поэтому горизонтальные
напряжения развёрнуты к осям разреза по формулам поворота:

#formula[$Delta = A_(italic("sec")) - A_H$]

#formula[$sigma_s = sigma_H dot cos^2 Delta + sigma_h dot sin^2 Delta$]

#formula[$sigma_o = sigma_H dot sin^2 Delta + sigma_h dot cos^2 Delta$]

где $A_(italic("sec"))$ -- азимут разреза, $A_H$ -- азимут действия максимального
горизонтального напряжения, $sigma_s$ -- горизонтальное напряжение, действующее
вдоль плоскости разреза, $sigma_o$ -- в перпендикулярном горизонтальном направлении.

Расчёт выполнен по девяти профильным линиям. Исходные данные приведены в #vref(<profile_stresses>, "п").
Для профильной линии №1, азимут которой
$A_(italic("sec")) = 0degree$, угол поворота и напряжения составляют:

#formula[$Delta = 0 - 70","3 = -70","3°$]

#formula[$sigma_s = 16","07 dot cos^2(-70","3) + 6","97 dot sin^2(-70","3) = 8","00 "МПа"$]

#formula[$sigma_o = 16","07 dot sin^2(-70","3) + 6","97 dot cos^2(-70","3) = 15","04 "МПа"$]

Вдоль плоскости разреза действует напряжение $sigma_s$. В качестве входного
параметра модели используется отношение этого напряжения к вертикальному:

#formula[$sigma_x slash sigma_y = sigma_s slash S_v$]

Для профильной линии №1 оно составляет:

#formula[$sigma_x slash sigma_y = 8","00 slash 10","791 = 0","741$]

Значения $Delta$, $sigma_s$ и $sigma_x slash sigma_y$ для остальных профильных
линий, рассчитанные по тем же формулам, приведены в #vref(<profile_stresses>, "п").
Расчёт для удобства выполнен на языке Python с использованием библиотеки
SymPy, исходный код программы приведён в #ref(<profile-stress-calc>, supplement: [приложении]).

#block(breakable: false)[
#let profile-stresses = csv("/assets/tables/profile_stresses.csv", delimiter: ";").slice(1)

#figure(
  kind: table,
  caption: [Отношение природного поля напряжений $sigma_x slash sigma_y$ для расчётных моделей],
  table(
    columns: (1.3fr, 0.8fr, 0.7fr, 0.9fr, 1fr),
    align: right,
    fill: (_, y) => if y == 0 { rgb("#ffd35f") } else { rgb("#ffff") },
    [*№ профильной линии *],
    [*Аз. разреза, °*],
    [*$Delta$, °*],
    [*$sigma_s$, МПа*],
    [*$sigma_x slash sigma_y$*],
    ..profile-stresses.flatten(),
  )
) <profile_stresses>
]

В ПК Prorock это отношение задаётся в окне Initial field stress для каждой
модели разреза, как показано на #vref(<field_stress>, "п").

#figure(
  image("/assets/images/ifs_window.png", width: 75%),
  caption: [Окно ввода параметров природного поля напряжений],
  gap: 10pt,
)<field_stress>

Примечание: Для материала закладки подземных выработок природное поле напряжений не
учитывалось.

=== Анализ чувствительности <sensitivity_analysis>

В ПК Prorock напряжённое состояние в модели задаётся одним из двух способов
@prorockInitialFieldStress. В расчёте применён способ с использованием
уровня земной поверхности. Природное поле при этом неоднородно,
вертикальная компонента в каждой точке вычисляется как $gamma H$ от уровня
земной поверхности, а задаётся соотношение компонент $sigma_x$ и $sigma_y$. 

Для оценки чувствительности задания начального поля напряжений выполнен
анализ, показывающий, как при изменении мощности перекрывающей
толщи $H$ меняется коэффициент бокового распора $lambda = sigma_x slash sigma_y$, совпадающий со входным параметром модели. Рассмотрены мощности перекрывающей толщи $H = 200$, $400$ и $600$ м.

В качестве примера принята профильная линия № 1
с азимутом $A_(italic("sec")) = 0degree$. Для неё $sigma_x = sigma_s$ --
горизонтальное напряжение вдоль плоскости разреза, а $sigma_y = S_v$ --
вертикальное напряжение. 

Расчёт для H = 400 м. Вертикальное напряжение при
$rho = 2750$ кг/м³ и $g = 9,81$ м/с² составляет:

#formula[$S_v = 2750 dot 9","81 dot 400 / 10^6 = 10","791 "МПа"$]

При $R = 0","42$, $A = 3","45$ и $P_p = 3","250$ МПа. Минимальное и максимальное
горизонтальные напряжения равны:

#formula[
  $S_(italic("hmin")) = frac(
    10","791 + 0","42 dot (3","45 - 1) dot 3","250,
    1 - 0","42 + 0","42 dot 3","45
  ) = 6","967 "МПа"$
]

#formula[
  $S_(italic("Hmax")) = 3","45 dot 6","967 - (3","45 - 1) dot 3","250
  = 16","072 "МПа"$
]

Для профильной линии № 1 угол поворота составляет
$Delta = 0 - 70","3 = -70","3degree$, поэтому:

#formula[
  $sigma_x = sigma_s = 16","072 dot cos^2(-70","3degree) +
  6","967 dot sin^2(-70","3degree) = 8","001 "МПа"$
]

#formula[
  $sigma_x slash sigma_y = 8","001 slash 10","791 = 0","741$
]

Расчёт для H = 200 м при $P_p = 1","750$ МПа:

#formula[
  $S_v = 2750 dot 9","81 dot 200 / 10^6 = 5","396 "МПа"$
]

#formula[
  $S_(italic("hmin")) = frac(
    5","396 + 0","42 dot (3","45 - 1) dot 1","750,
    1 - 0","42 + 0","42 dot 3","45
  ) = 3","547 "МПа"$
]

#formula[
  $S_(italic("Hmax")) = 3","45 dot 3","547 - (3","45 - 1) dot 1","750
  = 7","949 "МПа"$
]

#formula[
  $sigma_x = 7","949 dot cos^2(-70","3degree) +
  3","547 dot sin^2(-70","3degree) = 4","047 "МПа"$
]

#formula[
  $sigma_x slash sigma_y = 4","047 slash 5","396 = 0","750$
]

Расчёт для H = 600 м при $P_p = 4","875$ МПа:

#formula[
  $S_v = 2750 dot 9","81 dot 600 / 10^6 = 16","187 "МПа"$
]

#formula[
  $S_(italic("hmin")) = frac(
    16","187 + 0","42 dot (3","45 - 1) dot 4","875,
    1 - 0","42 + 0","42 dot 3","45
  ) = 10","450 "МПа"$
]

#formula[
  $S_(italic("Hmax")) = 3","45 dot 10","450 - (3","45 - 1) dot 4","875
  = 24","108 "МПа"$
]

#formula[
  $sigma_x = 24","108 dot cos^2(-70","3degree) +
  10","450 dot sin^2(-70","3degree) = 12","002 "МПа"$
]

#formula[
  $sigma_x slash sigma_y = 12","002 slash 16","187 = 0","741$
]

Сводные результаты анализа чувствительности показаны на #vref(<stress_sensitivity>, "п").

// График анализа чувствительности, построенный средствами Typst (place и
// line): по оси абсцисс -- мощность перекрывающей толщи H, по оси ординат --
// коэффициент бокового распора lambda = sigma_x/sigma_y.
// Точки читаются из stress_sensitivity.csv; десятичная запятая меняется
// на точку только для пересчёта в координаты, подписи точек выводятся
// из таблицы как есть. Помимо опорных расчётов (H = 200, 400, 600 м)
// таблица содержит промежуточные точки через 50 м, полученные по тем же
// формулам при поровом давлении, интерполированном между принятыми
// значениями, -- они сглаживают кривую; маркеры и подписи наносятся
// только на опорные точки.
#let sens-chart() = {
  let rows = csv("/assets/tables/stress_sensitivity.csv", delimiter: ";").slice(1)

  // Геометрия области построения: начало координат (левый нижний угол
  // сетки) и размеры, pt от левого верхнего угла блока.
  let (x0, y0) = (48pt, 155pt)
  let (w, h) = (180pt, 130pt)
  // Диапазоны осей: H, м и lambda.
  let (xmin, xmax) = (0.0, 650.0)
  let (ymin, ymax) = (0.735, 0.755)
  let sx(v) = x0 + (v - xmin) / (xmax - xmin) * w
  let sy(v) = y0 - (v - ymin) / (ymax - ymin) * h
  let pts = rows.map(row => (float(row.at(0)), float(row.at(2).replace(",", ".")), row.at(2)))
  // Опорные точки -- выполненные расчёты для H = 200, 400 и 600 м.
  let anchor-hs = (200, 400, 600)

  block(width: 253pt, height: 197pt, clip: true)[
    // Сетка, деления и подписи оси абсцисс (H, м)
    #for t in (0, 200, 400, 600) {
      place(dx: sx(t), dy: y0 - h, line(end: (0pt, h), stroke: 0.4pt + luma(225)))
      place(dx: sx(t), dy: y0, line(end: (0pt, 3pt), stroke: 0.6pt))
      place(dx: sx(t) - 15pt, dy: y0 + 6pt,
        box(width: 30pt, align(center)[#text(9pt)[#str(t)]]))
    }
    // Сетка, деления и подписи оси ординат (lambda)
    #for (t, lab) in ((0.735, "0,735"), (0.74, "0,740"), (0.745, "0,745"), (0.75, "0,750"), (0.755, "0,755")) {
      place(dx: x0, dy: sy(t), line(end: (w, 0pt), stroke: 0.4pt + luma(225)))
      place(dx: x0 - 3pt, dy: sy(t), line(end: (3pt, 0pt), stroke: 0.6pt))
      place(dx: x0 - 37pt, dy: sy(t) - 4.5pt,
        box(width: 33pt, align(right)[#text(9pt)[#lab]]))
    }
    // Оси со стрелками
    #place(dx: x0, dy: y0, line(end: (w + 14pt, 0pt), stroke: 0.8pt))
    #place(dx: x0, dy: y0, line(end: (0pt, -(h + 14pt)), stroke: 0.8pt))
    #place(dx: x0 + w + 20pt, dy: y0,
      polygon((0pt, 0pt), (-7pt, -2.6pt), (-7pt, 2.6pt), fill: black))
    #place(dx: x0, dy: y0 - h - 20pt,
      polygon((0pt, 0pt), (-2.6pt, 7pt), (2.6pt, 7pt), fill: black))
    // Названия осей
    #place(dx: x0 + w / 2 - 40pt, dy: y0 + 24pt,
      box(width: 80pt, align(center)[#text(12pt)[$H$, м]]))
    #place(dx: 4pt, dy: y0 - h / 2 - 25pt,
      rotate(-90deg, reflow: true, text(12pt)[$lambda = sigma_x slash sigma_y$]))
    // Линия зависимости и точки со значениями lambda
    #for (p, q) in pts.zip(pts.slice(1)) {
      place(dx: sx(p.at(0)), dy: sy(p.at(1)),
        line(end: (sx(q.at(0)) - sx(p.at(0)), sy(q.at(1)) - sy(p.at(1))), stroke: 1pt))
    }
    #for (hx, ly, lab) in pts {
      if hx in anchor-hs {
        place(dx: sx(hx) - 2.2pt, dy: sy(ly) - 2.2pt,
          circle(radius: 2.2pt, fill: black, stroke: none))
        place(dx: sx(hx) - 4pt, dy: sy(ly) - 14pt, text(9pt)[#lab])
      }
    }
  ]
}

#figure(
  kind: image,
  sens-chart(),
  caption: [Зависимость коэффициента бокового распора $lambda = sigma_x slash sigma_y$ от мощности перекрывающей толщи $H$],
  gap: 10pt,
) <stress_sensitivity>