diff --git a/E03_StabSection/E03_StabSection_Engine_Describe.py b/E03_StabSection/E03_StabSection_Engine_Describe.py index 6e3498b63..eb2e9e3d9 100644 --- a/E03_StabSection/E03_StabSection_Engine_Describe.py +++ b/E03_StabSection/E03_StabSection_Engine_Describe.py @@ -7,6 +7,7 @@ import numpy as np from E03_StabSection import E03_StabSection_Engine_GLE as gle from E03_StabSection import E03_StabSection_Engine_LEM as lem +from E03_StabSection import E03_StabSection_Engine_Support as support from E03_StabSection.E03_StabSection_Engine_Surface import ( circle_slices, evaluate, @@ -107,6 +108,10 @@ def describe(prep, kind, geom, method, n, min_depth, seed_F=None) -> dict: out["f0"] = float(res["f0"][0]) if prep.crack_depth > 0: out["crack_water_force"] = float(sl["Pw"][0]) + if prep.supports is not None: + out["support_forces"] = ( + support.forces(sl, 0, F, prep.supports["active"], s) if np.isfinite(F) else [] + ) return out diff --git a/E03_StabSection/E03_StabSection_Engine_GLE.py b/E03_StabSection/E03_StabSection_Engine_GLE.py index 7cee25909..c6c86962e 100644 --- a/E03_StabSection/E03_StabSection_Engine_GLE.py +++ b/E03_StabSection/E03_StabSection_Engine_GLE.py @@ -39,8 +39,8 @@ def forces(sl: dict, F, lam, f): g = 1.0 / F[:, None] m = sl["ca"] + g * p["sat"] A = sl["sa"] - g * p["tca"] - B0 = p["khws"] - g * p["cuca"] - rhs = sl["W"] - g * p["cusa"] + B0 = p["hz"] - g * p["cuca_n"] + rhs = p["wv"] - g * p["cusa_n"] # 보강 칸 포함(LEM `_pre`) S, n = rhs.shape N = np.empty((S, n)) den = np.empty((S, n)) diff --git a/E03_StabSection/E03_StabSection_Engine_LEM.py b/E03_StabSection/E03_StabSection_Engine_LEM.py index 33b74db07..e5ea8e44c 100644 --- a/E03_StabSection/E03_StabSection_Engine_LEM.py +++ b/E03_StabSection/E03_StabSection_Engine_LEM.py @@ -7,11 +7,14 @@ m_α = cosα + sinα·tanφ/F · A = sinα − tanφ·cosα/F · B0 = kh·Ws − (c·l − u·l·tanφ)·cosα/F λ = 절편간 힘 기울기(X = λ·E) — Bishop · Janbu 는 λ = 0 · Spencer 는 모멘트 · 힘 평형을 함께 맞추는 λ. 모멘트 F_m = Σ a·R / D (점 O 둘레) · 힘 F_f = Σ R·cosα / (Σ(N·sinα + kh·Ws) + Pw) · R = c·l + (N − u·l)·tanφ +보강재(`Engine_Support`) = 절편 바닥 점하중 — 능동은 W · kh·Ws 칸에 더함 · 수동(힘 / F)은 c' 칸(1/F 곱)에서 빼고 +전체 평형의 수동 항은 분자로 옮김(F = R/(D − T/F) ⇔ F = (R + T)/D · 같은 근). 외부 코드 없음(NumPy 만) — 상용 · 오픈소스 결과는 대조 시험에서만 씀. """ import numpy as np +from E03_StabSection import E03_StabSection_Engine_Support as support from E03_StabSection.E03_StabSection_Types import GAMMA_W, groundwater, material, surcharges _TINY = 1e-9 @@ -76,6 +79,7 @@ class Prepared: crack = model.get("tension_crack") or {} self.crack_depth = float(crack.get("depth") or 0.0) self.crack_water = float(crack.get("water") or 0.0) + self.supports = support.prepare(model.get("supports"), self.direction) def _build_columns(self, polys: list, pad: int) -> None: xs = np.unique(np.concatenate([e[:, 0] for e, _ in polys])) @@ -221,7 +225,7 @@ def slices(prep: Prepared, xe: np.ndarray, ze: np.ndarray) -> dict: valid = np.all(found & (zb >= prep.bottom_z - 1e-6) & np.isfinite(zt), axis=1) valid &= np.all(b > 0, axis=1) valid &= ~np.any(prep.blocked[base_mat] & found, axis=1) - return { + out = { "xe": xe, "ze": ze, "b": b, @@ -246,6 +250,9 @@ def slices(prep: Prepared, xe: np.ndarray, ze: np.ndarray) -> dict: "Pw": np.zeros(len(xe)), "zp": ze[:, 0].copy(), } + if prep.supports is not None: + support.apply(out, prep.supports) + return out def add_crack_water(sl: dict, depth: np.ndarray, fill: float) -> None: @@ -264,7 +271,11 @@ def _pre(sl: dict) -> dict: t, sa, ca = sl["tanphi"], sl["sa"], sl["ca"] cu = sl["c"] * sl["l"] - sl["u"] * sl["l"] * t khws = sl["kh"] * sl["Ws"] - lw = khws.copy() + # 보강 — 능동은 연직 · 수평 하중에 더함 · 수동(× 1/F)은 c' 칸(g 곱)에서 뺌 · 없으면 같은 배열 + sup = "sup_va" in sl + wv = sl["W"] + sl["sup_va"] if sup else sl["W"] + hz = khws + sl["sup_ha"] if sup else khws + lw = hz.copy() lw[:, 0] += sl["Pw"] # 머리 경계 수압은 전단 없음 → X_0 = 0 p = { "cu": cu, @@ -274,8 +285,12 @@ def _pre(sl: dict) -> dict: "tca": t * ca, "khws": khws, "lw": lw, - "fd0": khws.sum(1) + sl["Pw"], - "rc0": (cu * ca).sum(1), + "fd0": hz.sum(1) + sl["Pw"], + "rc0": (cu * ca).sum(1) - (sl["sup_hp"].sum(1) if sup else 0.0), + "wv": wv, + "hz": hz, + "cusa_n": cu * sa - sl["sup_vp"] if sup else cu * sa, + "cuca_n": cu * ca - sl["sup_hp"] if sup else cu * ca, } sl["_pre"] = p return p @@ -293,25 +308,35 @@ def _mom(sl: dict, center: np.ndarray): gn = dz * sl["sa"] - dx * sl["ca"] d0 = ((sl["xw"] - xo) * sl["W"] + (sl["zg"] - zo) * p["khws"]).sum(1) d0 = d0 + (sl["zp"] - center[:, 1]) * sl["Pw"] - out = (a * sl["tanphi"], gn, d0, (a * p["cu"]).sum(1)) + ma = (a * p["cu"]).sum(1) + if "sup_va" in sl: # 능동 = 활동 모멘트 · 수동(/F) = 분자로 옮김 + d0 = d0 + _sup_moment(sl, xo, zo, "a") + ma = ma - _sup_moment(sl, xo, zo, "p") + out = (a * sl["tanphi"], gn, d0, ma) sl["_mom"] = (center, out) return out +def _sup_moment(sl, xo, zo, tag): + """보강 점하중의 활동 모멘트 Σ (x − xo)·연직(아래 +) + (z − zo)·수평(+x).""" + v, h = sl[f"sup_v{tag}"], sl[f"sup_h{tag}"] + return (sl[f"sup_xv{tag}"] - xo * v + sl[f"sup_zh{tag}"] - zo * h).sum(1) + + def _normal(sl, F, lam, grad=0): """절편 수직력 N · 분모(m_α − λA) · grad ≥ 1 dN/dg · grad 2 dN/dλ — F · λ 꼴 (S,).""" p = _pre(sl) g = 1.0 / np.asarray(F, float)[:, None] L = np.asarray(lam, float)[:, None] if L.any(): - P1 = p["cusa"] + L * p["cuca"] + P1 = p["cusa_n"] + L * p["cuca_n"] Q1 = p["sat"] + L * p["tca"] den = sl["ca"] - L * sl["sa"] + g * Q1 - top = sl["W"] + L * p["lw"] + top = p["wv"] + L * p["lw"] else: # λ = 0(Bishop · Janbu) — 같은 식에서 λ 칸만 뺌 - P1, Q1 = p["cusa"], p["sat"] + P1, Q1 = p["cusa_n"], p["sat"] den = sl["ca"] + g * Q1 - top = sl["W"] + top = p["wv"] safe = np.where(np.abs(den) < _TINY, _TINY, den) N = (top - g * P1) / safe if not grad: @@ -319,7 +344,7 @@ def _normal(sl, F, lam, grad=0): dg = -(P1 + N * Q1) / safe if grad == 1: return N, den, dg, None - dl = (p["lw"] - g * p["cuca"] + N * (sl["sa"] - g * p["tca"])) / safe + dl = (p["lw"] - g * p["cuca_n"] + N * (sl["sa"] - g * p["tca"])) / safe return N, den, dg, dl @@ -334,7 +359,7 @@ def _parts(sl: dict, F, lam): N, den = _normal(sl, F, lam) g = 1.0 / F[:, None] A = sl["sa"] - p["tca"] * g - B0 = p["khws"] - p["cuca"] * g + B0 = p["hz"] - p["cuca_n"] * g return N, p["cu"] + N * sl["tanphi"], den, A, B0 @@ -408,11 +433,27 @@ def _iterate(sl, center, lam, F0, moment=True): return F, iters, done & np.isfinite(F) +def _ordinary_normal(sl: dict, F=None): + """Fellenius 바닥 수직 분해 N = W·cosα − kh·Ws·sinα (+ 보강 · 수동은 / F).""" + p = _pre(sl) + N = p["wv"] * sl["ca"] - p["hz"] * sl["sa"] + if "sup_vp" in sl and F is not None: + N = N + (sl["sup_vp"] * sl["ca"] - sl["sup_hp"] * sl["sa"]) / np.asarray(F, float)[:, None] + return N + + def ordinary(sl: dict, center: np.ndarray) -> dict: - """Fellenius — 절편간 힘 무시 · 바닥 수직 분해 N = W·cosα − kh·Ws·sinα · 모멘트 평형.""" - N = sl["W"] * sl["ca"] - _pre(sl)["khws"] * sl["sa"] + """Fellenius — 절편간 힘 무시 · 바닥 수직 분해 · 모멘트 평형(수동 보강이 있으면 N 이 F 를 따라 고정점 반복).""" at, gn, d0, ma = _mom(sl, center) - F = _ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1)) + F = np.full(len(sl["W"]), 1.5) + for _ in range(MAX_ITER if "sup_vp" in sl and sl["sup_vp"].any() else 1): + N = _ordinary_normal(sl, F) + F_new = _ratio(ma + (at * N).sum(1), d0 + (gn * N).sum(1)) + done = np.abs(F_new - F) < TOL * np.maximum(1.0, np.abs(F)) + F = np.where(np.isfinite(F_new) & (F_new > 0), F_new, F) + if done.all(): + break + F = F_new ok = ( sl["valid"] & np.isfinite(F) @@ -604,7 +645,7 @@ def slice_table( """ one = _subset(sl, np.arange(len(sl["W"])) == k) if method == "ordinary": - N = one["W"] * one["ca"] - one["kh"] * one["Ws"] * one["sa"] + N = _ordinary_normal(one, np.array([F])) E = np.zeros((1, N.shape[1] + 1)) X = E elif method == "morgenstern_price": diff --git a/E03_StabSection/E03_StabSection_Engine_Search.py b/E03_StabSection/E03_StabSection_Engine_Search.py index 145d4e529..d67fb9de7 100644 --- a/E03_StabSection/E03_StabSection_Engine_Search.py +++ b/E03_StabSection/E03_StabSection_Engine_Search.py @@ -13,6 +13,7 @@ import time import numpy as np from E03_StabSection import E03_StabSection_Engine_LEM as lem +from E03_StabSection import E03_StabSection_Engine_Support as support from E03_StabSection.E03_StabSection_Engine_Describe import describe, safety_map from E03_StabSection.E03_StabSection_Engine_Surface import ( # noqa: F401 — 밖에서 Search 로 부르던 이름 METHODS, @@ -45,8 +46,13 @@ FAR_TOL = 0.02 # ---- 원호 --------------------------------------------------------------------- +def _exit_lo(limits) -> float: + """진출 아래 한계 — limits 셋째 칸(보강으로 진입만 넓힌 때 원래 한계) · 없으면 limits[0].""" + return limits[2] if len(limits) > 2 else limits[0] + + def entry_exit_circles(prep, limits, n_points=25, thetas_deg=(10, 20, 35, 50, 65, 80, 100, 120)): - xl, xr = limits + xl, xr = limits[:2] pts = np.linspace(xl, xr, n_points) i, j = np.triu_indices(n_points, k=1) xE, xX = pts[i], pts[j] @@ -54,11 +60,12 @@ def entry_exit_circles(prep, limits, n_points=25, thetas_deg=(10, 20, 35, 50, 65 xX = np.repeat(xX, len(thetas_deg)) th = np.tile(np.radians(thetas_deg), len(i)) xc, zc, R, ok = circle_from_chord(prep, xE, xX, th) + ok &= xX >= _exit_lo(limits) - 1e-6 return xc[ok], zc[ok], R[ok], xE[ok], xX[ok] def grid_centers(prep, limits, nx=15, nz=15): - xl, xr = limits + xl, xr = limits[:2] z_in = prep.surface(np.linspace(xl, xr, 50)) H = max(np.nanmax(z_in) - np.nanmin(z_in), 1.0) L = max(xr - xl, 1.0) @@ -80,7 +87,7 @@ def grid_radius_circles(prep, limits, gx, gz, n_radii=10): owner = np.repeat(np.arange(len(cx)), n_radii) xc, zc = cx[owner], cz[owner] xE, xX, ok = circle_cross_surface(prep, xc, zc, R) - ok &= (xE >= limits[0] - 1e-6) & (xX <= limits[1] + 1e-6) + ok &= (xE >= limits[0] - 1e-6) & (xX <= limits[1] + 1e-6) & (xX >= _exit_lo(limits) - 1e-6) return xc[ok], zc[ok], R[ok], xE[ok], xX[ok], owner[ok] @@ -99,7 +106,7 @@ def _circle_eval(prep, xc, zc, R, xE, xX, n, min_depth, methods, chunk=4000): def refine_circle(prep, limits, method, seeds, n, min_depth, seed=0, maxiter=None): """원호 다듬기 — (진입 x, 진출 x, 중심각) 차분진화 · 첫 무리 = 거친 탐색의 좋은 면.""" fun = objective(method) - lo = np.array([limits[0], limits[0], np.radians(3)]) + lo = np.array([limits[0], _exit_lo(limits), np.radians(3)]) hi = np.array([limits[1], limits[1], np.radians(170)]) def obj(p): @@ -129,7 +136,7 @@ def _cloud(prep, limits, method, p0, n, min_depth, seed, k=CLOUD_SIZE, keep=10): """원 (진입 x, 진출 x, 중심각) 둘레 구름 k 개를 method 로 한 번에 셈 → 풀리는 좋은 것 keep 개.""" rng = np.random.default_rng(seed) span = limits[1] - limits[0] - lo = np.array([limits[0], limits[0], np.radians(3)]) + lo = np.array([limits[0], _exit_lo(limits), np.radians(3)]) hi = np.array([limits[1], limits[1], np.radians(170)]) p = np.clip(p0 + rng.normal(size=(k, 3)) * [0.03 * span, 0.03 * span, np.radians(4)], lo, hi) xc, zc, R, ok = circle_from_chord(prep, p[:, 0], p[:, 1], p[:, 2]) @@ -301,9 +308,9 @@ def optimize_polyline( ): """비원호 꼭짓점 최적화 — 차분진화 뒤 가장 좋은 면을 Nelder-Mead 로 국부 다듬기(상용 SAO · Optimize 대응).""" fun = objective(method) - xl, xr = limits + xl, xr = limits[:2] depth_max = max(prep.surface(np.linspace(xl, xr, 50)).max() - prep.bottom_z, 1.0) - lo = np.array([xl, xl] + [0.02] * m + [0.0] * m) + lo = np.array([xl, _exit_lo(limits)] + [0.02] * m + [0.0] * m) hi = np.array([xr, xr] + [0.98] * m + [depth_max] * m) def obj(p): @@ -370,7 +377,7 @@ def _ray_to_surface(prep, x0, z0, angle_deg, sign): def block_seeds(prep, limits, m, n_points=12, angles=((55, 35), (65, 25))): """블록 씨앗(Slide Block Search 대응) — 지층마다 바닥선 위 두 점을 잇는 토막 + 머리 쪽 능동 쐐기 · 끝 쪽 수동 쐐기를 지표까지 올림. 안쪽 꼭짓점 = 토막 두 끝 + 사이 (m − 2) 점.""" - xl, xr = limits + xl, xr = limits[:2] xs = np.linspace(xl, xr, n_points) i, j = np.triu_indices(n_points, k=1) xa, xb = xs[i], xs[j] @@ -386,6 +393,7 @@ def block_seeds(prep, limits, m, n_points=12, angles=((55, 35), (65, 25))): xE = _ray_to_surface(prep, a, za, act, -1.0) xX = _ray_to_surface(prep, b, zb, pas, 1.0) good = np.isfinite(xE) & np.isfinite(xX) & (xE >= xl) & (xX <= xr) + good &= xX >= _exit_lo(limits) for e, x0, x1, xx in zip(xE[good], a[good], b[good], xX[good]): ix = np.linspace(x0, x1, m) d = layer_bottom_depth(prep, ix, mat) @@ -421,6 +429,13 @@ def analyze( lim = (float(prep.surface_x[0]), float(prep.surface_x[-1])) else: lim = tuple(sorted((s * float(limits[0]), s * float(limits[1])))) + if prep.supports is not None: + # 보강 — 임계면이 보강재 끝 뒤(더 깊이 · 더 멀리)로 옮겨 가는 것도 찾게 진입 한계만 끝 뒤까지 넓힘 + # (진출 = 원래 한계 셋째 칸 — 설계 비탈 밖 자연 비탈만 무너지는 면은 보강 전처럼 안 봄) + reach = float(np.min(support.tips(prep.supports))) - max( + 2.0, 0.25 * prep.supports["L"].max() + ) + lim = (max(min(lim[0], reach), float(prep.surface_x[0])), lim[1], lim[0]) results: dict = {m: {} for m in methods} counts = {} contour = smap = None @@ -546,7 +561,7 @@ def analyze( if v.get("F") is not None and not v.get("reference_only") ] results[m]["min"] = min(cands, key=lambda v: v["F"]) if cands else None - return { + out = { "direction": s, "limits": sorted([s * lim[0], s * lim[1]]), "methods": results, @@ -555,6 +570,16 @@ def analyze( "counts": counts, "elapsed_s": round(time.perf_counter() - t0, 3), } + if prep.supports is not None: + # 보강력 결과 = 판정 임계면(방법별 최소 중 가장 낮은 면) 기준 + mins = [(m, v["min"]) for m, v in results.items() if v.get("min")] + if mins: + m, gov = min(mins, key=lambda kv: kv[1]["F"]) + out["support_forces"] = gov.get("support_forces", []) + out["support_forces_method"] = m + else: + out["support_forces"], out["support_forces_method"] = [], None + return out def _stack(prep, kind, geoms, n, min_depth): diff --git a/E03_StabSection/E03_StabSection_Engine_Station.py b/E03_StabSection/E03_StabSection_Engine_Station.py index 79da300be..2d4a0c3b9 100644 --- a/E03_StabSection/E03_StabSection_Engine_Station.py +++ b/E03_StabSection/E03_StabSection_Engine_Station.py @@ -28,12 +28,15 @@ def analyze_station( min_depth: float = MIN_DEPTH_M, case_workers: int = 1, groundwater_for=None, + supports=None, **model_kw, ) -> dict: """측점 하나 — case_workers > 1 이면 비탈 × 조건 경우를 프로세스로 나눠 돎(측점 하나만 다시 셀 때). groundwater_for(condition, model) → 엔진 지하수 꼴 — 주면 조건마다 모형 지하수를 갈아 끼움 (E02 `station_groundwater` 꼴) · None 이면 모형 기본(우기 = 지표 수위). + supports = 네일 · 앵커(`Engine_Support` 꼴) — 보강재마다 "side"(없으면 머리 x 가 있는 쪽 = 노면 가운데 + 기준 좌 +)의 비탈 경우에만 넣음 · 결과 cases[i].result.support_forces · None 이면 지금 그대로. 서버 프로세스 안에서는 1(차례로) 그대로 둠 — 자식이 main.py 를 다시 읽음(`Worker` 머리글). """ @@ -54,7 +57,9 @@ def analyze_station( if groundwater_for is not None: model["groundwater"] = norm_groundwater(groundwater_for(cond, model)) out.setdefault("areas", model["areas"]) - jobs += [(cond, model, side) for side in model["sides"]] + for side in model["sides"]: + own = _side_supports(supports, side["side"], design) + jobs.append((cond, {**model, "supports": own} if own else model, side)) kw = {"methods": methods, "noncircular": noncircular, "min_depth": min_depth} if case_workers > 1 and len(jobs) > 1: from concurrent.futures import ProcessPoolExecutor @@ -71,11 +76,29 @@ def analyze_station( ] for (cond, _m, side), res in zip(jobs, results): mins = {m: (res["methods"][m]["min"] or {}).get("F") for m in methods} + if supports: + res.setdefault("support_forces", []) # 보강재 없는 쪽 비탈 = 빈 칸 out["cases"].append({"condition": cond, **side, "min_F": mins, "result": res}) out["elapsed_s"] = round(time.perf_counter() - t0, 3) return out +def _side_supports(supports, side: str, design: dict) -> list[dict]: + """그 비탈 보강재 — "side" 칸 · 없으면 머리 x 가 노면 가운데(좌 · 우 노면 끝 가운데)보다 + 면 좌.""" + if not supports: + return [] + edges = design.get("road_edges") or {} + mid = 0.5 * ( + float((edges.get("left") or {}).get("offset_m", 2.0)) + + float((edges.get("right") or {}).get("offset_m", -2.0)) + ) + return [ + s + for s in supports + if (s.get("side") or ("left" if float(s["x"]) >= mid else "right")) == side + ] + + def analyze_stations(items: list[dict], workers: int | None = None, **kw) -> list[dict]: """전 측점 — 측점마다 따로 도니 프로세스로 나눠 돎. items = [{"cross", "design", "ground"?}, …] · 순서 그대로.""" import os diff --git a/E03_StabSection/E03_StabSection_Engine_Support.py b/E03_StabSection/E03_StabSection_Engine_Support.py new file mode 100644 index 000000000..4da53b2fa --- /dev/null +++ b/E03_StabSection/E03_StabSection_Engine_Support.py @@ -0,0 +1,152 @@ +"""E03 보강재(네일 · 앵커) — 활동면을 지나는 보강력을 절편 바닥 점하중으로(Slide2 Support 와 같은 꼴). + +입력 꼴(E03 보강 검토 화면과 한 벌) supports = [{ + "type": "nail" | "anchor", "x", "z" (머리 · 모형 좌표 m), "angle_deg" (수평 아래로 +), "length_m", + "spacing_m" (종방향 간격 · 단위 폭 힘 = 힘 / 간격), "tensile_kN" (강재 인장 한도), + "bond_kPa" (그라우트 · 지반 부착 · 네일 = 전 길이 · 앵커 = 끝 정착장만), "bond_length_m" (앵커 정착장), + "hole_dia_mm", "prestress_kN" (앵커만 · 0 이면 수동), "plate_kN"? (머리 판 · 없으면 tensile_kN), + "side"? ("left" | "right" · 고르는 것은 `Engine_Station`)}] +힘 크기 = 지점마다 Slide2 Force Diagram 의 최소 — 인장(tensile) · 인발(pullout = 부착 × 활동면 밖 부착 길이) · + 벗겨짐(stripping = 판 + 부착 × 활동면 안 부착 길이) · 부착 = π·D·bond_kPa (kN/m) · 자유장은 부착 없음. + 인발 안전율은 입력 부착값에 넣어 둠(Slide2 와 같음 — 엔진은 나누지 않음). +방향 = 보강재와 나란히(Slide2 Parallel to Reinforcement) · 비탈 속으로(활동 방향 반대 = −x) 아래로 angle. +능동 · 수동 = Slide2 Active/Passive Force Application: + 능동(프리스트레스 앵커) = 힘 그대로 — 활동력을 줄임 F = R / (D − T) · 힘 = min(프리스트레스, 한도) + 수동(네일 · 프리스트레스 0 앵커) = 힘 / F — 저항을 늘림 F = (R + T) / D +머리가 활동 덩이 밖이면(활동면이 보강재를 안 지나거나 덩이 밖에서 드나듦) 0 = not_crossing. +절편 칸(꼴 (S, n) · 활동 방향 +x 틀): sup_va · sup_ha = 능동 연직(아래 +) · 수평(+x) · sup_vp · sup_hp = 수동(× 1/F 전) + · sup_xva · sup_zha · sup_xvp · sup_zhp = 모멘트용 x·연직 · z·수평. 보강재 칸(꼴 (S, K)) sup_T · sup_mode · sup_x · sup_z. +""" + +import numpy as np + +MODES = ("tensile", "pullout", "stripping", "not_crossing", "prestress") +_NOT, _PRE = 3, 4 +_FIELDS = ("va", "ha", "vp", "hp", "xva", "zha", "xvp", "zhp") + + +def prepare(raw: list[dict] | None, direction: int) -> dict | None: + """입력 꼴 → 배열(활동 방향 +x 틀 · x 에 direction 곱함) · 없으면 None.""" + if not raw: + return None + s = 1.0 if direction >= 0 else -1.0 + cols = {k: [] for k in ("hx", "hz", "dx", "dz", "L", "free", "q", "tens", "plate", "sp", "pre")} + active = [] + for r in raw: + kind = r.get("type", "nail") + if kind not in ("nail", "anchor"): + raise ValueError(f"보강재 종류 모름: {kind}") + L, sp = float(r["length_m"]), float(r["spacing_m"]) + tens, bond = float(r["tensile_kN"]), float(r["bond_kPa"]) + dia = float(r["hole_dia_mm"]) / 1000.0 + lb = float(r.get("bond_length_m") or 0.0) if kind == "anchor" else L + pre = float(r.get("prestress_kN") or 0.0) if kind == "anchor" else 0.0 + plate = float(r["plate_kN"]) if r.get("plate_kN") is not None else tens + if L <= 0 or sp <= 0 or dia <= 0 or tens < 0 or bond < 0 or plate < 0 or pre < 0: + raise ValueError(f"보강재 값 범위 밖: {r}") + if not 0 < lb <= L: + raise ValueError(f"앵커 정착장은 0 초과 · 길이 이하: {r}") + a = np.radians(float(r["angle_deg"])) + for k, v in ( + ("hx", s * float(r["x"])), + ("hz", float(r["z"])), + ("dx", -np.cos(a)), + ("dz", -np.sin(a)), + ("L", L), + ("free", L - lb), + ("q", np.pi * dia * bond), + ("tens", tens), + ("plate", plate), + ("sp", sp), + ("pre", pre), + ): + cols[k].append(v) + active.append(pre > 0) + out = {k: np.array(v, float) for k, v in cols.items()} + out["active"] = np.array(active, bool) + return out + + +def tips(sup: dict) -> np.ndarray: + """보강재 끝 x (활동 방향 +x 틀) — 탐색 한계를 끝 뒤까지 넓힘.""" + return sup["hx"] + sup["L"] * sup["dx"] + + +def apply(sl: dict, sup: dict) -> None: + """절편 묶음(sl)에 보강력 칸을 덧붙임 — 면마다 보강재마다 첫 교점 · 힘 · 파괴 꼴.""" + xe, ze = sl["xe"], sl["ze"] + S, n = sl["b"].shape + hx, hz, dx, dz = (sup[k][None, :, None] for k in ("hx", "hz", "dx", "dz")) + ax, az = xe[:, None, :-1], ze[:, None, :-1] + ex, ez = np.diff(xe, axis=1)[:, None, :], np.diff(ze, axis=1)[:, None, :] + den = dx * ez - dz * ex + safe = np.where(np.abs(den) < 1e-12, np.nan, den) + wx, wz = ax - hx, az - hz + t = (wx * ez - wz * ex) / safe # 보강재 위 머리에서 거리 + u = (wx * dz - wz * dx) / safe # 활동면 마디 위 비율 + hit = (u >= 0) & (u < 1) & (t > 0) & (t <= sup["L"][None, :, None]) + tt = np.where(hit, t, np.inf) + j = np.argmin(tt, axis=-1) # (S, K) 첫 교점 절편 + tc = np.take_along_axis(tt, j[..., None], -1)[..., 0] + # 머리가 활동 덩이 안(진입 ~ 진출 사이 · 활동면 위)일 때만 — 첫 교점 = 덩이에서 나가는 점 + k_head = np.clip((xe[:, None, :] <= sup["hx"][None, :, None]).sum(-1) - 1, 0, n - 1) + x0 = np.take_along_axis(xe, k_head, 1) + x1 = np.take_along_axis(xe, k_head + 1, 1) + z0 = np.take_along_axis(ze, k_head, 1) + z1 = np.take_along_axis(ze, k_head + 1, 1) + zslip = z0 + (sup["hx"][None, :] - x0) * (z1 - z0) / np.maximum(x1 - x0, 1e-12) + inside = (sup["hx"] >= xe[:, :1]) & (sup["hx"] <= xe[:, -1:]) & (sup["hz"] > zslip) + cross = inside & np.isfinite(tc) & sl["valid"][:, None] + tc = np.where(cross, tc, 0.0) + lb = sup["L"] - sup["free"] + bond_in = np.clip(tc - sup["free"], 0.0, lb) + caps = np.stack( + [ + np.broadcast_to(sup["tens"], tc.shape), + sup["q"] * (lb - bond_in), + sup["plate"] + sup["q"] * bond_in, + ], + axis=-1, + ) + mode = np.argmin(caps, axis=-1) + T = np.min(caps, axis=-1) + pre_rules = sup["active"] & (sup["pre"] < T) + T = np.where(pre_rules, sup["pre"], T) + mode = np.where(pre_rules, _PRE, mode) + T = np.where(cross, T, 0.0) / sup["sp"] + mode = np.where(cross, mode, _NOT) + xs = sup["hx"] + tc * sup["dx"] + zs = sup["hz"] + tc * sup["dz"] + acc = {f: np.zeros((S, n)) for f in _FIELDS} + rows = np.arange(S) + for k in range(len(sup["L"])): + v, h = -T[:, k] * sup["dz"][k], T[:, k] * sup["dx"][k] # 아래 + · +x(비탈 속 = −) + tag = "a" if sup["active"][k] else "p" + for f, val in ((f"v{tag}", v), (f"h{tag}", h), (f"xv{tag}", xs[:, k] * v)): + np.add.at(acc[f], (rows, j[:, k]), val) + np.add.at(acc[f"zh{tag}"], (rows, j[:, k]), zs[:, k] * h) + for f in _FIELDS: + sl[f"sup_{f}"] = acc[f] + sl["sup_T"], sl["sup_mode"], sl["sup_x"], sl["sup_z"] = T, mode, xs, zs + + +def forces(sl: dict, k: int, F: float, active: np.ndarray, direction: int) -> list[dict]: + """면 k 의 보강재별 결과 칸(support_forces) — 한도(단위 폭) · 실제 넣은 힘(수동 = / F) · 파괴 꼴 · 교점.""" + out = [] + for i, (T, m) in enumerate(zip(sl["sup_T"][k], sl["sup_mode"][k])): + crossing = int(m) != _NOT + out.append( + { + "index": i, + "mobilized_kN_per_m": float(T), + "applied_kN_per_m": float(T if active[i] else T / F), + "mode": MODES[int(m)], + "force_application": "active" if active[i] else "passive", + "point": ( + [float(direction * sl["sup_x"][k, i]), float(sl["sup_z"][k, i])] + if crossing + else None + ), + } + ) + return out diff --git a/resources/tester/test_536_3_stab_supports.py b/resources/tester/test_536_3_stab_supports.py new file mode 100644 index 000000000..a234dc796 --- /dev/null +++ b/resources/tester/test_536_3_stab_supports.py @@ -0,0 +1,271 @@ +# -*- coding: utf-8 -*- +"""536-3 뒷부분 보강재(네일 · 앵커) 2D 한계평형 — 힘 크기 · 능동/수동 · 탐색 넓힘 · 보강 없을 때 불변. + +- 손계산: 평면 쐐기 닫힌 꼴(능동 F = (cL + (W cosα + T sin(α+δ)) tanφ) / (W sinα − T cos(α+δ)) · 수동 T → T/F) + · 교점 · 인발 · 벗겨짐 · 정착장 · 프리스트레스. +- Slide2 Verification(Slide_SlopeStabilityVerification.pdf): #47 Amherst 네일 벽(수동 · 인발 · 평면 Janbu 0.890 + · Sheahan 0.887) · #85 Duncan & Wright 그림 6.34(앵커 1단 · 능동 Bishop 1.531 · GLE 1.575 · 수동 1.324 · + 1.378) · #86 그림 7.28(5단 · Bishop 1.629 · Spencer 1.620 · GLE 1.622). #85 Slide Spencer(1.884 · 1.872)는 + 같은 표 Bishop · GLE · 참값(1.51 · 1.32)과 25 % 넘게 어긋난 값이라 대조에서 뺌. + + ./venv/Scripts/python.exe -m pytest resources/tester/test_536_3_stab_supports.py -q +""" + +from __future__ import annotations + +import math +import sys +from pathlib import Path + +import pytest + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT)) +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from test_532_2_stab_section_model import DESIGN, GROUND # noqa: E402 + +from E03_StabSection.E03_StabSection_Engine_Search import analyze, evaluate_specified # noqa: E402 +from E03_StabSection.E03_StabSection_Engine_Station import analyze_station # noqa: E402 + +MP = "morgenstern_price" + + +def near(v, ref, tol): + assert v is not None and abs(v - ref) <= tol * abs(ref), (v, ref) + + +def nail(x, z, angle, L, sp, tens, bond_kn_per_m, **kw): + """부착 = 길이당 힘(kN/m) 으로 줌 — bond_kPa = 그 값 / (π·D) · D 100 mm.""" + return { + "type": "nail", "x": x, "z": z, "angle_deg": angle, "length_m": L, "spacing_m": sp, + "tensile_kN": tens, "bond_kPa": bond_kn_per_m / (math.pi * 0.1), "hole_dia_mm": 100, **kw, + } # fmt: skip + + +def anchor(x, z, angle, L, sp, tens, bond_kn_per_m, lb, pre, **kw): + return {**nail(x, z, angle, L, sp, tens, bond_kn_per_m, **kw), "type": "anchor", + "bond_length_m": lb, "prestress_kN": pre} # fmt: skip + + +# ---- 손계산: 평면 쐐기 --------------------------------------------------------------- +# 활동 +x · 비탈 (10, 10) → (20, 0) · 면 (4, 10) → 기슭 (20, 0) · 쐐기 넓이 30 m² +WEDGE = { + "materials": [{"id": "s", "gamma_t": 18.0, "c": 5.0, "phi": 30.0}], + "regions": [{"material_id": "s", "polygon": [ + [0, -5], [30, -5], [30, 0], [20, 0], [10, 10], [0, 10]]}], +} # fmt: skip +PLANE = [[4.0, 10.0], [20.0, 0.0]] +HEAD, DELTA = (15.0, 5.0), 15.0 + + +def _wedge_closed(T, passive): + a, d = math.atan2(10, 16), math.radians(DELTA) + W, L, c, t = 30 * 18.0, math.hypot(16, 10), 5.0, math.tan(math.radians(30)) + F = 1.5 + for _ in range(200): + Tf = T / F if passive else T + F = (c * L + (W * math.cos(a) + Tf * math.sin(a + d)) * t) / ( + W * math.sin(a) - Tf * math.cos(a + d) + ) + return F + + +@pytest.mark.parametrize("passive", [False, True]) +def test_wedge_matches_closed_form(passive): + """쐐기 하나 · 보강 한 개(인장 100 kN · 간격 2 m → 50 kN/m) — Janbu 간편 · Spencer = 닫힌 꼴(±0.1 %).""" + sup = ( + nail(*HEAD, DELTA, 12, 2, 100, 1e6) + if passive + else anchor(*HEAD, DELTA, 12, 2, 100, 1e6, 12, 100) + ) + r = evaluate_specified( + {**WEDGE, "supports": [sup]}, polyline=PLANE, methods=("janbu", "spencer") + ) + F = _wedge_closed(50.0, passive) + near(r["janbu"]["F0"], F, 1e-3) + near(r["spencer"]["F"], F, 1e-3) + f = r["janbu"]["support_forces"][0] + assert f["mode"] == "tensile" and f["mobilized_kN_per_m"] == pytest.approx(50.0) + assert f["force_application"] == ("passive" if passive else "active") + assert f["applied_kN_per_m"] == pytest.approx(50.0 / r["janbu"]["F0"] if passive else 50.0) + # 보강 효과 — 보강 없는 쐐기보다 높음 · 능동 > 수동(Slide2 와 같음) + bare = evaluate_specified(WEDGE, polyline=PLANE, methods=("janbu",))["janbu"]["F0"] + assert F > bare and _wedge_closed(50.0, False) > _wedge_closed(50.0, True) + + +def _crossing_t(): + """머리 (15, 5) 에서 −x · 아래 15° 로 간 거리 t 가 면 z = 10 − (x − 4)·10/16 과 만나는 곳.""" + c, s, k = math.cos(math.radians(DELTA)), math.sin(math.radians(DELTA)), 10 / 16 + # 5 − t·s = 10 − (15 − t·c − 4)·k + return (5 - 10 + 11 * k) / (s + k * c) + + +def test_capacity_modes_hand_calc(): + """교점 거리 t · 인발(부착 × 밖 길이) · 벗겨짐(판 + 부착 × 안 길이) · 앵커 정착장 · 프리스트레스 · 안 지남.""" + t = _crossing_t() + L, q = 8.0, 20.0 # 길이 · 길이당 부착 kN/m + sups = [ + nail(*HEAD, DELTA, L, 1, 500, q), # 인발 = q·(L − t) + nail(*HEAD, DELTA, L, 1, 500, q, plate_kN=10), # 벗겨짐 = 10 + q·t (< 인발) + anchor(*HEAD, DELTA, L, 1, 500, q, 3.0, 0), # 교점이 자유장 — 인발 = q·3(정착장 전부) + anchor(*HEAD, DELTA, L, 1, 500, q, 3.0, 40), # 능동 · 프리스트레스 40 < 한도 60 + nail(*HEAD, DELTA, t - 0.5, 1, 500, q), # 면에 못 미침 + nail(25.0, 0.0, 0.0, 6, 1, 500, q), # 머리가 덩이 밖 + ] + r = evaluate_specified({**WEDGE, "supports": sups}, polyline=PLANE, methods=("janbu",)) + out = r["janbu"]["support_forces"] + expect = [ + ("pullout", q * (L - t)), + ("stripping", 10 + q * t), + ("pullout", q * 3.0), + ("prestress", 40.0), + ("not_crossing", 0.0), + ("not_crossing", 0.0), + ] + assert [o["index"] for o in out] == list(range(len(sups))) + for o, (mode, T) in zip(out, expect): + assert o["mode"] == mode and o["mobilized_kN_per_m"] == pytest.approx(T, rel=1e-6) + x = HEAD[0] - t * math.cos(math.radians(DELTA)) + assert out[0]["point"] == pytest.approx([x, HEAD[1] - t * math.sin(math.radians(DELTA))]) + assert out[4]["point"] is None + + +# ---- Slide2 Verification ------------------------------------------------------------- +def _p47(): + return { + "materials": [{"id": "a", "gamma_t": 18.9, "c": 25, "phi": 0}], + "regions": [{"material_id": "a", "polygon": [ + [-6, -2], [11, -2], [11, 6.1], [0, 6.1], [0, 0], [-6, 0]]}], + "surcharge": [{"x0": 0, "x1": 0.02, "q": 14.6 / 0.02}], # 숏크리트 무게 14.6 kN/m + "supports": [nail(0, z, 20, 4.9, 1.5, 118, 15, plate_kN=86) for z in (4.9, 3.4)], + } # fmt: skip + + +def test_slide2_47_soil_nail_wall_planar(): + """수동 네일 2단 · 인발이 지배 · 평면 면 Janbu 0.890(Sheahan 0.887) · 탐색도 같은 평면.""" + r = evaluate_specified( + _p47(), polyline=[[0, 0], [6.279, 6.1]], direction=-1, methods=("janbu",) + ) + near(r["janbu"]["F"], 0.890, 0.03) + near(r["janbu"]["F"], 0.887, 0.01) + k = 6.1 / 6.279 + c, s = math.cos(math.radians(20)), math.sin(math.radians(20)) + for f, z in zip(r["janbu"]["support_forces"], (4.9, 3.4)): + t = z / (s + k * c) + assert f["mode"] == "pullout" + assert f["mobilized_kN_per_m"] == pytest.approx(15 * (4.9 - t) / 1.5, rel=1e-6) + a = analyze(_p47(), direction=-1, methods=("bishop", "janbu"), limits=[-6, 11]) + near(a["methods"]["janbu"]["min"]["F"], 0.890, 0.03) + + +def _p85(active): + return { + "materials": [{"id": "c", "gamma_t": 98, "c": 350, "phi": 0}], # ft · psf · pcf + "regions": [{"material_id": "c", "polygon": [[15, 10], [57, 10], [57, 30], [25, 30]]}], + "supports": [anchor(20, 20, 0, 37, 1, 9000, 1e7, 37, 9000 if active else 0, plate_kN=1e9)], + } + + +@pytest.mark.parametrize( + "active, bishop, gle", + [(True, 1.531, 1.575), (False, 1.324, 1.378)], +) +def test_slide2_85_active_passive(active, bishop, gle): + a = analyze(_p85(active), direction=-1, methods=("bishop", "spencer", MP)) + ms = a["methods"] + near(ms["bishop"]["min"]["F"], bishop, 0.03) + near(ms[MP]["min"]["F"], gle, 0.03) + near(ms["bishop"]["min"]["F"], 1.51 if active else 1.32, 0.03) # Duncan & Wright 참값 + f = a["support_forces"][0] + assert f["mode"] == "tensile" and f["mobilized_kN_per_m"] == pytest.approx(9000) + assert f["force_application"] == ("active" if active else "passive") + + +def _p86(capacity=800): + return { + "materials": [{"id": "f", "gamma_t": 130, "c": 0, "phi": 37}], + "regions": [{"material_id": "f", "polygon": [[0, 0], [75, 0], [75, 24], [30, 24]]}], + "supports": [ + anchor(5 * k, 4 * k, 0, 20, 1, capacity, 1e7, 20, capacity, plate_kN=1e9) + for k in range(1, 6) + ], + } + + +def test_slide2_86_five_rows(): + """Slide 임계 원(중심 (−3.383, 40.466) · R 40.568) 값 · 탐색(전체 붕괴만 — 얕은 면 3 m 아래 뺌).""" + r = evaluate_specified(_p86(), circle=[-3.383, 40.466, 40.568], direction=-1) + near(r["bishop"]["F"], 1.629, 0.03) + near(r["spencer"]["F"], 1.620, 0.03) + near(r[MP]["F"], 1.622, 0.03) + assert all(f["mode"] == "tensile" for f in r["bishop"]["support_forces"]) + a = analyze(_p86(), direction=-1, methods=("bishop", "spencer"), min_depth=3.0) + near(a["methods"]["bishop"]["min"]["F"], 1.629, 0.03) + near(a["methods"]["spencer"]["min"]["F"], 1.620, 0.03) + + +def test_search_moves_behind_strong_supports(): + """아주 센 보강 — 진입 쪽 한계가 보강재 끝 뒤까지 넓어지고 임계면이 보강재를 안 지나는 깊은 면으로 감.""" + bare = {k: v for k, v in _p86().items() if k != "supports"} + lim = [0, 45] + a0 = analyze(bare, direction=-1, methods=("bishop",), min_depth=3.0, limits=lim) + a = analyze(_p86(1e5), direction=-1, methods=("bishop",), min_depth=3.0, limits=lim) + assert a0["limits"] == [0, 45] and a["limits"][1] > 45 # 끝 x 45 + 여유 5 + gov = a["methods"]["bishop"]["min"] + assert gov["entry"][0] > 45 and a["support_forces_method"] == "bishop" + assert all(f["mode"] == "not_crossing" for f in a["support_forces"]) + assert gov["F"] > a0["methods"]["bishop"]["min"]["F"] + + +# ---- 보강 없을 때 불변 · 측점 꼴 ------------------------------------------------------ +def test_no_supports_unchanged(): + """supports 없음 · 빈 목록 · 덩이 밖 보강재 = 같은 값 · 결과 칸 support_forces 는 보강 줄 때만.""" + base = evaluate_specified(WEDGE, polyline=PLANE) + empty = evaluate_specified({**WEDGE, "supports": []}, polyline=PLANE) + far = evaluate_specified({**WEDGE, "supports": [nail(25, 0, 0, 3, 1, 100, 50)]}, polyline=PLANE) + for m, v in base.items(): + assert empty[m]["F"] == v["F"] and far[m]["F"] == pytest.approx(v["F"], abs=1e-12) + assert "support_forces" not in v and far[m]["support_forces"][0]["mode"] == "not_crossing" + + +def _station(supports): + cross = {"chainage_m": 0.0, "samples": [{"offset_m": x, "elevation_m": z} for x, z in GROUND]} + design = { + "section_mode": "right_cut", + "design_line": [{"offset_m": x, "elevation_m": z} for x, z in DESIGN], + "road_edges": {"left": {"offset_m": 2.0}, "right": {"offset_m": -2.0}}, + } + return analyze_station( + cross, + design, + conditions=("dry",), + methods=("bishop",), + noncircular=False, + supports=supports, + ) + + +def test_station_side_and_result_slot(): + """우 깎기 면 네일 — 우 비탈 경우에만 힘 · 좌 비탈은 보강 없는 값 그대로 · cases[i].result.support_forces.""" + plain = _station(None) + sups = [nail(-4.0, 11.48, 15, 8, 1.5, 200, 30)] + res = _station(sups) + by = {c["side"]: c for c in res["cases"]} + old = {c["side"]: c for c in plain["cases"]} + assert "support_forces" not in old["right"]["result"] + assert ( + by["left"]["min_F"] == old["left"]["min_F"] and by["left"]["result"]["support_forces"] == [] + ) + f = by["right"]["result"]["support_forces"] + assert len(f) == 1 and set(f[0]) >= {"index", "mobilized_kN_per_m", "mode"} + assert by["right"]["min_F"]["bishop"] > old["right"]["min_F"]["bishop"] + # "side" 칸이 머리 x 보다 먼저 — 좌로 돌리면 좌 비탈에만(좌 성토 면 밖이라 not_crossing) + left = _station([{**sups[0], "side": "left"}]) + by_l = {c["side"]: c for c in left["cases"]} + assert by_l["right"]["result"]["support_forces"] == [] + assert by_l["right"]["min_F"] == old["right"]["min_F"] + + +if __name__ == "__main__": + sys.exit(pytest.main([__file__, "-q"]))