Building/assets/code/profile_stress_ratio.py

"""Расчёт напряжений в осях разрезов и отношения sigma_x / sigma_y.

Природное поле напряжений задано в главных осях: максимальное горизонтальное
напряжение S_Hmax действует по азимуту A_H, минимальное S_hmin -- в
перпендикулярном горизонтальном направлении, вертикальное S_v занимает
промежуточное положение. Расчётная модель строится в плоскости вертикального
разреза, поэтому горизонтальные напряжения поворачиваются к осям разреза по
формулам преобразования компонент тензора при повороте вокруг вертикали:

    Delta = A_sec - A_H,
    sigma_s = sigma_H cos^2(Delta) + sigma_h sin^2(Delta),
    sigma_o = sigma_H sin^2(Delta) + sigma_h cos^2(Delta),

где A_sec -- азимут разреза; sigma_s действует вдоль плоскости разреза,
sigma_o -- в перпендикулярном горизонтальном направлении. Входным параметром
модели служит отношение sigma_x / sigma_y = sigma_s / S_v, задаваемое в окне
Initial field stress ПК Prorock.

Скрипт формирует три результата.

1. Сводная таблица по линиям (profile_stresses.csv) на опорной глубине
   H = 400 м: азимуты разрезов читаются из самой этой таблицы (столбец
   «Азимут разреза»), по ним пересчитываются Δ, σ_s и σ_x/σ_y; столбец
   азимутов при записи сохраняется. Напряжения S_Hmax = 16,07 МПа,
   S_hmin = 6,97 МПа и S_v = 10,791 МПа -- из раздела «Результаты расчёта
   абсолютных значений природного поля напряжений»

2. Расчётная таблица по разрезам (section_stresses.csv): для каждой
   профильной линии на каждый год контура (2026, 2031, 2033, 2035 гг.)
   природное поле оценивается на глубине карьера H, а горизонтальные
   напряжения определяются из системы уравнений того же раздела:

       S_Hmax = (S_v - (1 - R) S_hmin) / R,     -- коэффициент формы R,
       S_Hmax = A S_hmin - (A - 1) P_p,         -- фрикционное ограничение,

   где S_hmin = (S_v + R (A - 1) P_p) / (1 - R + R A); S_v = rho g H,
   поровое давление принято пропорциональным глубине: P_p = 3,25 МПа · H / 400
   (на опорной глубине 400 м -- принятое в отчёте значение 3,25 МПа).
   Азимут линии берётся из profile_stresses.csv, глубина карьера -- из
   profile_azimuths.csv.

3. Таблица анализа чувствительности (stress_sensitivity.csv) для профильной
   линии № 1. В ней сравнивается отношение sigma_x/sigma_y на глубинах от
   200 до 600 м с шагом 50 м. Поровое давление задано равным 1,750, 3,250
   и 4,875 МПа на опорных глубинах 200, 400 и 600 м, на промежуточных
   глубинах интерполируется квадратично по трём опорным значениям.

Исходные данные: assets/tables/profile_stresses.csv -- азимуты разрезов
по профильным линиям; assets/tables/profile_azimuths.csv -- глубина карьера
по разрезам (год контура, номер линии, глубина). Вычисления ведутся в точной
рациональной арифметике SymPy.

"""
import csv
from fractions import Fraction
from pathlib import Path

import sympy as sp

# Азимут действия максимального горизонтального напряжения, °.
A_H = sp.Rational(703, 10)
# Коэффициент формы тензора и фрикционное ограничение A = (sqrt(1+mu^2)+mu)^2.
R = sp.Rational(42, 100)
A = sp.Rational(345, 100)
# Средняя плотность толщи, кг/м3, и ускорение свободного падения, м/с2.
RHO, G = 2750, sp.Rational(981, 100)

# Опорное поле на глубине H = 400 м, МПа (раздел «Результаты расчёта
# абсолютных значений природного поля напряжений»). 
P_REF = sp.Rational(325, 100)
S_V_REF = sp.Rational(10791, 1000)
S_HMIN_REF = (S_V_REF + R*(A - 1)*P_REF) / (1 - R + R*A)
S_HMAX_REF = A*S_HMIN_REF - (A - 1)*P_REF

S_v, P_p, Delta = sp.symbols("S_v P_p Delta", real=True)

# Система двух уравнений решена относительно S_hmin в символьном виде.
S_hmin = (S_v + R*(A - 1)*P_p) / (1 - R + R*A)
S_Hmax = A*S_hmin - (A - 1)*P_p

# Напряжения в осях разреза как функции угла поворота.
sigma_s = S_Hmax*sp.cos(Delta)**2 + S_hmin*sp.sin(Delta)**2
sigma_o = S_Hmax*sp.sin(Delta)**2 + S_hmin*sp.cos(Delta)**2


def rotated(sigma_H, sigma_h, azimuth):
    """Напряжения в осях разреза при заданных главных напряжениях."""
    d = sp.rad(sp.Rational(azimuth) - A_H)
    s = sigma_H*sp.cos(d)**2 + sigma_h*sp.sin(d)**2
    o = sigma_H*sp.sin(d)**2 + sigma_h*sp.cos(d)**2
    return sp.Rational(azimuth) - A_H, s, o


def section(azimuth, depth, pp):
    """Поле на глубине карьера для данного разреза.

    Возвращает (S_v, S_Hmax, S_hmin, sigma_s, sigma_o, sigma_x/sigma_y).
    """
    subs = {S_v: RHO*G*depth/10**6, P_p: pp}
    d = sp.rad(sp.Rational(azimuth) - A_H)
    sv = S_v.subs(subs)           # вертикальное напряжение на глубине карьера
    sh = S_hmin.subs(subs)        # минимальное горизонтальное
    sH = S_Hmax.subs(subs)        # максимальное горизонтальное
    s = sigma_s.subs({**subs, Delta: d})
    o = sigma_o.subs({**subs, Delta: d})
    return sv, sH, sh, s, o, s/sv


def fmt(value, digits):
    """Число с десятичной запятой и типографским знаком минус."""
    text = f"{float(sp.N(value, 30)):.{digits}f}"
    return text.replace(".", ",").replace("-", "\u2212")


def num(cell):
    """Число из ячейки CSV с десятичной запятой."""
    return Fraction(cell.replace(",", "."))


def read_azimuths(root):
    """Азимуты линий: {номер линии: азимут} из profile_stresses.csv."""
    path = root / "assets/tables/profile_stresses.csv"
    azimuths = {}
    with open(path, encoding="utf-8", newline="") as f:
        for row in csv.reader(f, delimiter=";"):
            if len(row) == 5 and row[0].isdigit():
                azimuths[int(row[0])] = int(row[1])
    return azimuths


def read_sections(root):
    """Строки таблицы глубин: [(год, линия, H), ...]."""
    path = root / "assets/tables/profile_azimuths.csv"
    sections = []
    with open(path, encoding="utf-8", newline="") as f:
        for row in csv.reader(f, delimiter=";"):
            if len(row) == 3 and row[0].isdigit():
                sections.append((row[0], int(row[1]), num(row[2])))
    return sections


def write_csv(path, rows):
    with open(path, "w", encoding="utf-8", newline="") as f:
        csv.writer(f, delimiter=";").writerows(rows)


def reference_table(root, azimuths):
    """Сводная таблица по линиям на опорной глубине 400 м."""
    print("опорная глубина 400 м (поле из отчёта)")
    rows = [("№ профильной линии", "Азимут разреза, °", "Δ, °",
             "σ_s, МПа", "σ_x/σ_y")]
    for line in sorted(azimuths):
        delta, s, o = rotated(S_HMAX_REF, S_HMIN_REF, azimuths[line])
        print(f"линия {line}  азимут {azimuths[line]}  Δ={fmt(delta,1)}"
              f"  σ_s={fmt(s,2)}  σ_o={fmt(o,2)}")
        rows.append((str(line), str(azimuths[line]), fmt(delta, 1),
                     fmt(s, 2), fmt(s/S_V_REF, 3)))
    write_csv(root / "assets/tables/profile_stresses.csv", rows)


def section_table(root, sections, azimuths):
    """Расчётная таблица по разрезам: контур, линия, H, P_p."""
    print("\nрасчёт по разрезам (глубина карьера из исходной таблицы, P_p = 3,25·H/400 МПа)")
    print(f"{'год':>5} {'лин':>4} {'азимут':>7} {'H,м':>5} {'Pp,МПа':>7} "
          f"{'S_Hmax':>8} {'S_hmin':>8} {'σ_s,МПа':>8} {'σ_x/σ_y':>8}")
    rows = [("Год контура", "№ профильной линии", "Азимут разреза, °",
             "Глубина карьера, м", "P_p, МПа", "Δ, °", "σ_s, МПа", "σ_x/σ_y")]
    for year, line, depth in sections:
        if line not in azimuths:
            raise SystemExit(f"линия {line} отсутствует в profile_stresses.csv")
        azimuth = azimuths[line]
        # Давление привязано к опорному значению 3,25 МПа на глубине 400 м
        # и принято пропорциональным глубине карьера.
        pp = P_REF*depth/400
        sv, sH, sh, s, o, k = section(azimuth, depth, pp)
        print(f"{year:>5} {line:>4} {azimuth:>7} {fmt(depth,0):>5} "
              f"{fmt(pp,2):>7} {fmt(sH,2):>8} {fmt(sh,2):>8} "
              f"{fmt(s,2):>8} {fmt(k,3):>8}")
        rows.append((year, str(line), str(azimuth), fmt(depth, 0), fmt(pp, 2),
                     fmt(sp.Rational(azimuth) - A_H, 1), fmt(s, 2), fmt(k, 3)))
    write_csv(root / "assets/tables/section_stresses.csv", rows)


def stability_table(root):
    """Анализ чувствительности sigma_x/sigma_y относительно H"""
    # Линия № 1 имеет азимут разреза 0°. Опорное давление при H = 200 м
    # задано как 1750 кПа (в расчёте используются МПа),
    #Точки выводятся через 50 м; на промежуточных глубинах поровое давление интерполируется
    # квадратично по трём опорным значениям.
    anchors = ((200, sp.Rational(1750, 1000)),
               (400, sp.Rational(3250, 1000)),
               (600, sp.Rational(4875, 1000)))

    def pp_at(h):
        (h1, p1), (h2, p2), (h3, p3) = anchors
        return (p1*(h - h2)*(h - h3)/((h1 - h2)*(h1 - h3))
                + p2*(h - h1)*(h - h3)/((h2 - h1)*(h2 - h3))
                + p3*(h - h1)*(h - h2)/((h3 - h1)*(h3 - h2)))

    print("\nанализ чувствительности: профильная линия 1 (азимут 0°)")
    print(f"{'H,м':>5} {'Pp,МПа':>7} {'S_v':>8} {'S_Hmax':>8} "
          f"{'S_hmin':>8} {'σ_s,МПа':>8} {'σ_x/σ_y':>8}")
    rows = [("H, м", "P_p, МПа", "σ_x/σ_y")]
    for depth in range(200, 601, 50):
        pp = pp_at(sp.Rational(depth))
        sv, sH, sh, s, _, k = section(0, depth, pp)
        print(f"{depth:>5} {fmt(pp,3):>7} {fmt(sv,3):>8} {fmt(sH,2):>8} "
              f"{fmt(sh,2):>8} {fmt(s,2):>8} {fmt(k,3):>8}")
        rows.append((str(depth), fmt(pp, 3), fmt(k, 3)))
    write_csv(root / "assets/tables/stress_sensitivity.csv", rows)


def main():
    root = Path(__file__).resolve().parents[2]
    azimuths = read_azimuths(root)
    sections = read_sections(root)
    reference_table(root, azimuths)
    section_table(root, sections, azimuths)
    stability_table(root)


if __name__ == "__main__":
    main()