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

= Методика численного моделирования

== Метод конечно-дискретных элементов

Метод конечно-дискретных элементов (англ. FDEM, finite-discrete element method) разработан в конце 90-х -- начале 2000-х гг. группой британских учёных во главе с А. Муньизой. Подробное описание метода приведено в диссертации Ильясова Б. Т. @ilyasovIssledovanieKinetikiDeformatsiy2016.

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

В методе конечно-дискретных элементов материал модели представляется в виде треугольных конечных элементов, соединённых четырёхугольными трещинными элементами. Сетка элементов изображена на #vref(<FDEM_Mesh>, "предл"). В ненапряжённом состоянии трещинные элементы имеют нулевую толщину, поэтому на рисунке конечные элементы показаны уменьшенными, а трещинные -- расширенными.

#figure(
  image("/assets/images/FDEM_Mesh.jpg", width: 60%),
  caption:[Конечно-дискретно-элементная сетка]
)<FDEM_Mesh>

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

#formula(
  $x_c=x_p+v dot Delta t,$
) <eq_node_position>

где $x_c$ -- положение узла в пространстве в текущий момент времени; $x_p$ -- положение узла в предыдущий момент времени; $v$ -- скорость узла; $Delta t$ -- временной шаг.

Скорость узла в общем случае определяется в соответствии с выражением:

#formula(
 $Delta v=F dot Delta t/m,$
) <eq_node_velocity>

в котором $F$ -- сила, действующая на данный узел; $m$ -- масса узла (треть массы треугольного конечного элемента). Полная скорость получается прибавлением этого приращения к скорости предыдущего шага.

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

#formula(
 $F=m dot g + F_"int" + F_"ext" + F_c +C_v$
) <eq_node_force>

где $g$ -- ускорение свободного падения; $F_"int"$ -- внутренние силы; $F_"ext"$ -- внешние приложенные силы; $F_c$ -- сила контактного взаимодействия; $C_v$ -- постоянная вязкостного демпфирования.
Внутренние силы $F_"int"$ находятся сложением сил $F_e$, вызванных упругой деформацией конечных элементов, и сил связи $F_j$, возникающих в трещинных элементах.

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

Обобщённо алгоритм обнаружения контактов можно представить в виде четырёх этапов:

#numbered-list(scheme: "decimal")[
  + Разбивка области поиска прямоугольной сеткой. Для оптимизации расчётов размер ячейки подбирается несколько больше окружности, описанной вокруг наибольшего элемента.
  + Запись дискретных элементов в столбцы сетки.
  + Запись дискретных элементов в отдельные ячейки сетки.
  + Нахождение потенциально контактирующих элементов внутри текущей ячейки и в соседних четырёх ячейках.
]

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

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

Изменение напряжений с увеличением смещений при сдвиге до и после достижения предельного смещения описывается полной диаграммой деформирования, изображённой на 
#vref(<Strain_stress_curve>, "предл").

#figure(
  image("/assets/images/Strain_stress_curve.jpg", width: 60%),
  caption:[Полная диаграмма деформирования при сдвиговом разрушении])
<Strain_stress_curve>

$f_"s"$ -- прочность трещинного элемента на сдвиг, $tau$ -- касательное напряжение, $s_p$ -- предельное касательное смещение, при котором достигается пиковая прочность трещинного элемента. Оно задаётся формулой
#formula(
  $s_p= (2 dot h dot f_s) / p$
)

где $h$ -- размер трещинного элемента, $p$ -- параметр штрафа, $s_r$ -- касательное смещение в момент окончательного разрушения элемента, то есть при полной потере связности по принятой схеме разупрочнения. Через $s_p$ и $s_r$ задаётся степень разрушения:

#formula(
  $D=(s-s_p) / (s_r - s_p)$
)

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

Прочность на сдвиг рассчитывается по критерию Кулона:

#formula(
 $f_s=C+sigma_n dot tg phi$
) <eq_shear_strength>

Остаточная прочность рассчитывается по формуле:

#formula(
 $f_r=sigma_n dot tg phi_r$
) <eq_residual_strength>

Метод конечно-дискретных элементов позволяет использовать не только дискретизованные сетки, в которых каждый конечный элемент отделён от другого трещинным элементом, но и конечно-элементные сетки, в которых узлы в вершинах элементов являются общими. Эта особенность метода широко используется в ПО Prorock для обеспечения корректности задания граничных условий с сохранением вычислительной эффективности моделирования. Главное окно ПО Prorock представлено на #vref(<Prorock_main_window>, "предл").

== Особенности реализации метода в Prorock <features>

#figure(
  image("/assets/images/Prorock_main_window.png", width: 100%),
  caption:[Главное окно ПО Prorock])
<Prorock_main_window>

Реализация FDEM в Prorock отличается алгоритмами и улучшениями оптимизации, направленными на значительное повышение производительности и точности моделирования. Часть из них описана в диссертации Ильясова Б. Т. @ilyasovIssledovanieKinetikiDeformatsiy2016, наиболее важные особенности приведены ниже.

=== Выполнение расчётов на графических ускорителях общего назначения

Главным недостатком метода конечно-дискретных элементов являются высокие требования к вычислительным мощностям. Чтобы обеспечить ускорение расчётов при приемлемой дискретизации моделей, в Prorock применяются параллельные вычисления. Код решателя написан на языке C++ с поддержкой платформы CUDA. Применение CUDA вместе с современными графическими устройствами позволяет выполнять расчёты в разы быстрее, чем на обычных процессорах, поэтому количество элементов в моделях также может быть увеличено.

=== Оптимизация алгоритма контактов

Алгоритм поиска контактов в Prorock работает значительно более эффективно благодаря избирательному применению поиска. Данный алгоритм обеспечивает заметное ускорение расчётов (ориентировочно в 1,5--2 раза) ценой некоторого допустимого снижения стабильности вычислений, что сказывается на удобстве применения программы, но не на точности вычислений.

=== Алгоритм принудительной стабилизации модели

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

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

=== Учёт дилатансии и влияния типа разрушения на параметры сдвиговой прочности

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

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

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

=== Учёт порового давления <pore_pressure>

В Prorock для анализа гидростатического воздействия задаётся уровень воды. Давление имеет положительное значение в дискретизированной области ниже этого уровня. Окно настройки порового давления представлено на #vref(<Pore_pressure_window>, "предл").

#figure(
  image("/assets/images/Pore_pressure_window.png", width: 60%),
  caption:[Окно настройки порового давления])
<Pore_pressure_window>

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

#formula(
 $P=gamma dot h dot H_u$
) <eq_pore_pressure>


где $gamma = rho dot g$ -- удельный вес флюида, $h$ -- вертикальное расстояние от узла до уровня грунтовых вод, $H_u$ -- коэффициент гидростатического воздействия, $rho$ -- плотность флюида, $g$ -- ускорение свободного падения.

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

#formula(
  $F_3=(3 dot P_a + P_b) dot L_1/8$
)

#formula(
  $F_2=(P_a + 3 dot P_b) dot L_1/8$
)

#formula(
  $F_4=(3 dot P_a + P_b) dot L_2/8$
)

#formula(
  $F_5=(P_a + 3 dot P_b) dot L_2/8$
)

=== Метод снижения прочностных свойств <srf_method>

Метод снижения прочностных свойств (Strength Reduction Method) -- это численный метод, используемый для расчёта коэффициента запаса устойчивости (Factor of Safety, FoS) области модели, например целика или откоса. Его суть заключается в последовательном и одновременном уменьшении прочностных характеристик массива (сцепления -- C, угла внутреннего трения -- $phi$ и предела прочности на растяжение -- $"UTS"$) до тех пор, пока массив не достигнет состояния предельного равновесия (т. е. произойдёт разрушение). Окно настройки представлено на #vref(<Srf>, "предл").

#figure(
  image("/assets/images/Srf.png", width: 45%),
  caption:[Окно настройки метода снижения прочностных свойств])
<Srf>

Формулы @hammahShearStrengthReduction2005 для пересчёта свойств представлены ниже. $"SRF"$ -- безразмерный коэффициент снижения прочности, значение которого последовательно увеличивается по шагам процедуры снижения прочности.

#formula(
  $C _"current"= 1/"SRF" dot C_"init"$
)


#formula(
  $tan(phi _"current")= 1/"SRF" dot tan(phi_"init")$
)


#formula(
  $"UTS" _"current"= 1/"SRF" dot "UTS"_"init"$
)

=== Импорт данных из блочной модели <bm_import>

В Prorock доступен импорт параметров физико-механических свойств конечных элементов $rho$ (плотность), $E$ (модуль Юнга) и трещинных элементов $C$ (сцепление), $phi$ (угол трения), $"UTS"$ (предел прочности при растяжении). Окно импорта блочной модели представлено на #vref(<BM_import_window>, "предл"). Процесс импорта заключается в переносе значений параметров материалов из ячеек блочной модели (далее -- БМ) в конечные и трещинные элементы сетки модели в Prorock.

#figure(
  image("/assets/images/BM_import_window.png", width: 75%),
  caption:[Окно импорта блочной модели])
<BM_import_window>



=== Геометрическое соответствие 

Для корректного импорта необходимо установить соответствие между 2D разрезом модели Prorock и 3D сечением БМ вертикальной плоскостью. Данное соответствие устанавливается с помощью двух пар опорных точек: в каждой паре точка на 2D разрезе соответствует точке на 3D разрезе. При этом важно соблюдение следующих правил:

- Выбирать точки необходимо слева направо;
- Вертикальные координаты точек в каждой паре должны быть одинаковыми ($y_"2D" = z_"3D"$);
- Две точки на 2D разрезе и две точки на 3D разрезе не должны находиться на одной вертикали (для задания направления секущей плоскости);
- Расстояния между точками на 2D разрезе и на 3D разрезе должны отличаться не более чем на 1% (для предотвращения ошибок при масштабировании).

При выполнении всех условий в Prorock запускается процедура импорта.

=== Критерии прочности

Импорт возможен по двум критериям прочности. 

- Критерий Мора--Кулона -- прямой импорт значений параметров из ближайших ячеек БМ. 
- Критерий Хука--Брауна (Hoek--Brown) -- прямой импорт обязательных вспомогательных параметров $("UCS",m_i,"GSI",D,sigma_3 )$, по которым вычисляются $E,C,phi,"UTS"$ (при этом $rho$ можно импортировать так же напрямую, как в случае выбора критерия Мора--Кулона).

=== Определение параметров с помощью критерия Хука--Брауна <calculation>
В обобщённом критерии Хука--Брауна @hoekHoekBrownFailureCriterionGSI2019 параметры $m_b$, $s$ и $a$ представляют собой материальные константы скального массива, учитывающие снижение прочности массива относительно интактной породы под влиянием трещиноватости, структурной нарушенности и техногенного повреждения. Параметр $m_b$ является аналогом константы $m_i$ для интактной породы, приведённым к условиям массива, и характеризует чувствительность прочности массива к всестороннему сжатию, влияя на наклон и общую форму огибающей разрушения. Параметр $s$ отражает степень сохранности и структурной целостности массива и фактически задаёт уровень снижения его прочности по сравнению с интактной породой. Параметр $a$ является показателем степени в уравнении критерия и определяет нелинейность огибающей прочности, то есть характер её кривизны при изменении бокового давления. Значения $m_b$, $s$, $a$ определяются по эмпирическим зависимостям через индекс геологической прочности $"GSI"$, константу интактной породы $m_i$ и коэффициент нарушенности массива $D$.


Вспомогательные величины $(sigma_"ci"= "UCS")$ @hoekHoekBrownFailureCriterionGSI2019:

#formula($
m_b = m_i exp(("GSI" - 100) / (28 - 14 D)),
$)
#formula($
s = exp(("GSI" - 100) / (9 - 3 D)),
$)

#formula($
a = 1 / 2 + 1 / 6 (e^(-"GSI" / 15) - e^(-20 / 3)),
$)

#formula($
sigma_(3n) = sigma_3 / sigma_"ci".
$)

Расчёт $E,C, phi, "UTS" (sigma_t = "UTS")$ @hoekRockMassProperties:

Модуль упругости Юнга: 
#formula($
E = 10^5 (1 - D/2) / (1 + exp((75 + 25 dot D - "GSI") / 11))
$)

Угол внутреннего трения:
#formula($
phi = sin^(-1)(
  (6 dot a dot m_b (s + m_b sigma_(3n))^(a - 1)) /
  (2 (1 + a) (2 + a) + 6 dot a dot m_b dot (s + m_b dot sigma_(3n))^(a - 1))
)
$)
Сцепление:
#formula($
C = (
  sigma_"ci" [(1 + 2a) dot s + (1 - a) dot m_b dot sigma_(3n)] (s + m_b dot sigma_(3n))^(a - 1)
) / (
  (1 + a) (2 + a) sqrt(
    1 + [6 dot a dot m_b (s + m_b dot sigma_(3n))^(a - 1)] / [(1 + a) (2 + a)]
  )
)
$)