"""Средний коэффициент Hu порового давления по моделям Slide2. В ПК Prorock гидростатическое воздействие задаётся формулой P = gamma_w h H_u (раздел «Учёт порового давления»), где h -- вертикальное расстояние от узла до заданной поверхности грунтовых вод, H_u -- коэффициент вклада поровой воды. Поровое давление, рассчитанное в Slide2 фильтрационным анализом (метод Groundwater), заменяется одним коэффициентом H_u, наилучшим образом воспроизводящим рассчитанное поле. Поверхность грунтовых вод в Prorock при этом задаётся по рельефу дневной поверхности, поэтому h -- глубина узла ниже рельефа. Файл .slim -- zip-архив модели Slide2, содержащий три текстовых файла: .sli -- модель: удельный вес воды gamma_w (ключ gammaw), координаты вершин геометрии (секция vertices) и внешний контур (exterior), верхняя огибающая которого задаёт рельеф дневной поверхности; .slw -- сетка фильтрационного расчёта: координаты узлов (секция nodes); .w01 -- результаты расчёта: строки «номер узла, полный напор, пьезометри- ческий напор» в метрах; поровое давление u = gamma_w h_p, кПа. Для каждого узла сетки вычисляются глубина h ниже рельефа и поровое давление u, после чего строятся три оценки коэффициента: 1. H_u по МНК через начало координат u = H_u gamma_w h (узлы с h >= 1 м): H_u = sum(u h) / (gamma_w sum(h^2)); минимизирует среднеквадратичную невязку поля давлений -- основное значение для переноса в Prorock; 2. среднее и медиана поточечных отношений H_u,i = u / (gamma_w h) по узлам с h >= 5 м (порог отсекает неустойчивые отношения у поверхности); 3. среднее H_u,i в полосе изобары 1 МПа (900--1100 кПа) -- аналог прежнего ручного способа, в котором коэффициент определялся по глубине залегания изобары 1000 кПа. Запуск (из корня репозитория): python assets/code/hu_from_slide.py "E:\\...\\Исходники_Slide" [ещё пути...] Аргументы командной строки -- файлы .slim или каталоги; для каталога обрабатываются все файлы .slim в его подкаталогах. Результаты печатаются в консоль и записываются в assets/tables/hu_from_slide.csv. Имя модели «35 1.slim» разбирается как контур 2035 года и профильная линия 1. Используется только стандартная библиотека. """ import csv import math import re import statistics import sys import zipfile from pathlib import Path # Пороги глубины, м: для МНК и для поточечных отношений H_u,i. H_MIN_FIT, H_MIN_POINT = 1.0, 5.0 # Полоса изобары 1 МПа для сопоставления с прежним ручным способом, кПа. ISOBAR_LO, ISOBAR_HI = 900.0, 1100.0 # Префикс номера контура в имени файла -- год контура. CONTOUR_YEARS = {"26": 2026, "31": 2031, "33": 2033, "35": 2035} _NUM = r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?" def section(text, name): """Строки секции name. Заголовок -- слово с двоеточием в начале строки, данные -- всё до следующего заголовка (в .sli с отступом, в .slw нет).""" lines, inside = [], False for line in text.splitlines(): if line[:1].isspace() or line[:1].isdigit(): if inside: lines.append(line) else: inside = line.split(":")[0].strip() == name return lines def read_vertices(sli): """Координаты вершин геометрии: {номер: (x, y)}.""" rx = re.compile(rf"(\d+)\s+x:\s*({_NUM})\s+y:\s*({_NUM})") return {int(m[1]): (float(m[2]), float(m[3])) for m in map(rx.search, section(sli, "vertices")) if m} def read_exterior(sli): """Номера вершин внешнего контура (замкнутого многоугольника).""" for line in section(sli, "exterior"): m = re.search(r"\[([0-9,\s]+)\]", line) if m: return [int(v) for v in m[1].split(",")] return [] def surface_function(ids, verts): """Рельеф дневной поверхности: верхняя огибающая внешнего контура как функция от x.""" ring = [verts[i] for i in ids if i in verts] ring += ring[:1] # замыкание контура segments = [(a, b) for a, b in zip(ring, ring[1:]) if a[0] != b[0]] def y_at(x): top = None for (x1, y1), (x2, y2) in segments: lo, hi = (x1, x2) if x1 < x2 else (x2, x1) if lo <= x <= hi: y = y1 + (y2 - y1)*(x - x1)/(x2 - x1) top = y if top is None or y > top else top return top return y_at def read_nodes(slw): """Координаты узлов сетки фильтрационного расчёта: {номер: (x, y)}.""" rx = re.compile(rf"(\d+)\s+x:\s*({_NUM})\s+y:\s*({_NUM})") return {int(m[1]): (float(m[2]), float(m[3])) for m in map(rx.search, section(slw, "nodes")) if m} def read_heads(w01): """Напоры из результатов: {номер узла: (полный, пьезометрический)}, м.""" heads = {} rx = re.compile(rf"(\d+)\s+({_NUM})\s+({_NUM})$") for line in w01.split("# Fluxes")[0].splitlines(): # блок напоров m = rx.match(line.strip()) if m: heads[int(m[1])] = (float(m[2]), float(m[3])) return heads def model_stats(path): """Оценки H_u по одной модели .slim.""" with zipfile.ZipFile(path) as z: def pick(suffix): names = [n for n in z.namelist() if n.lower().endswith(suffix)] return z.read(names[0]).decode("cp1251", "replace") if names else "" sli, slw, w01 = pick(".sli"), pick(".slw"), pick(".w01") if not (sli and slw and w01): raise ValueError("нет .sli/.slw/.w01 -- модель без результатов расчёта") m = re.search(rf"gammaw:\s*({_NUM})", sli) gamma_w = float(m[1]) if m else 9.81 verts = read_vertices(sli) y_surface = surface_function(read_exterior(sli), verts) nodes, heads = read_nodes(slw), read_heads(w01) # Глубина ниже рельефа и поровое давление каждого узла, м и кПа. # Контроль сопоставления узлов: полный напор минус высота узла равен # пьезометрическому напору (допускается невязка до 0,5 м). pts = [] for n, (x, y) in nodes.items(): if n in heads: total, ph = heads[n] if abs(total - y - ph) > 0.5: raise ValueError("напоры не сходятся с координатами узлов") y_top = y_surface(x) if y_top is not None and y_top - y >= 0: pts.append((y_top - y, gamma_w*ph)) fit = [(h, u) for h, u in pts if h >= H_MIN_FIT] hu_lsm = sum(u*h for h, u in fit)/(gamma_w*sum(h*h for h, _ in fit)) rms = math.sqrt(sum((u - hu_lsm*gamma_w*h)**2 for h, u in fit)/len(fit)) point = [u/(gamma_w*h) for h, u in pts if h >= H_MIN_POINT and u > 0] isobar = [u/(gamma_w*h) for h, u in pts if h >= H_MIN_POINT and ISOBAR_LO <= u <= ISOBAR_HI] return { "n": len(fit), "h_max": max(h for h, _ in fit), "u_max": max(u for _, u in fit), "hu_lsm": hu_lsm, "rms": rms, "hu_mean": statistics.mean(point) if point else None, "hu_median": statistics.median(point) if point else None, "hu_isobar": statistics.mean(isobar) if isobar else None, } def model_name(stem): """«35 1» -> (год контура, номер линии) по имени файла модели.""" parts = stem.split() year = str(CONTOUR_YEARS.get(parts[0], parts[0])) return year, parts[1] if len(parts) > 1 else "--" def fmt(value, digits=3): """Число с десятичной запятой и типографским знаком минус.""" if value is None: return "--" return f"{value:.{digits}f}".replace(".", ",").replace("-", "\u2212") def main(): if len(sys.argv) < 2: sys.exit("укажите файлы .slim или каталоги с моделями Slide2") paths = [] for arg in sys.argv[1:]: p = Path(arg) paths.extend(sorted(p.rglob("*.slim")) if p.is_dir() else [p]) root = Path(__file__).resolve().parents[2] header = ("Год контура", "№ профильной линии", "Hu (МНК)", "Hu среднее", "Hu медиана", "Hu у изобары 1 МПа", "Узлов", "h_max, м", "u_max, кПа", "СКО поля, кПа") rows, values = [header], [] print(f"{'год':>4} {'лин':>4} {'Hu МНК':>8} {'среднее':>8} " f"{'медиана':>8} {'1 МПа':>8} {'узлов':>6} " f"{'h_max':>6} {'u_max':>7} {'СКО':>7}") for path in paths: try: st = model_stats(path) except (ValueError, KeyError, zipfile.BadZipFile) as e: print(f"{path.stem}: пропущена ({e})") continue year, line = model_name(path.stem) print(f"{year:>4} {line:>4} {fmt(st['hu_lsm']):>8} " f"{fmt(st['hu_mean']):>8} {fmt(st['hu_median']):>8} " f"{fmt(st['hu_isobar']):>8} {st['n']:>6} " f"{fmt(st['h_max'], 0):>6} {fmt(st['u_max'], 0):>7} " f"{fmt(st['rms'], 0):>7}") rows.append((year, line, fmt(st["hu_lsm"]), fmt(st["hu_mean"]), fmt(st["hu_median"]), fmt(st["hu_isobar"]), st["n"], fmt(st["h_max"], 1), fmt(st["u_max"], 0), fmt(st["rms"], 0))) values.append(st["hu_lsm"]) if values: print(f"среднее Hu (МНК) по всем моделям: {fmt(statistics.mean(values))}") out = root/"assets/tables/hu_from_slide.csv" with open(out, "w", encoding="utf-8", newline="") as f: csv.writer(f, delimiter=";").writerows(rows) print(f"таблица записана: {out}") if __name__ == "__main__": main()