"""Расчёт напряжений в осях разрезов и отношения 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()