eomsangdonandClaude Opus 5.5
f2e8bd6d5f
- 엔진(자체 · NumPy · SciPy): 평면 · 쐐기(교선) · 전도 Markland(측방 한계 20° 기본) · 평면 Hoek-Bray(인장균열 없음 · 임계 · 깊이 · 수압 · kh · 수동 볼트) · 쐐기 SWedge 꼴(미끄럼 방향 8 경우 · 두 면 · 한 면 · 수압 중력식 · kh · 볼트) - 스테레오넷 좌표(등각 · 등면적 하반구 · 극점 · 대원 · 교선 · 마찰원 · 일광 엔벨로프 · 위험 영역 = 판정 부등식) - 길: 절리 작업본 E06_StabRock/joints.json(초기값 · 저장 · 초기화) · 해석 = 암 측점 깎기 쪽만 · 파생 저장본 stab_rock · 결과 · 측점 스테레오넷 - 화면: 좌 절리군 표 · 설정 / 가운데 위험 띠 · 스테레오넷 / 우 측점 결과 · 안전율 표 · E 단계 막대에 E06(끝에 덧붙임) - 시험 test_536_1(Hoek Sau Mau Ping 0.8184 · 1.0655 · Sharma 윗면 무게 · Hoek-Bray 쐐기 계수 · Kumsar · Priest 두 예제 · Markland 경계 · 투영 왕복 · 위험 영역 = 판정) · test_528 E 단계 수 4 → 5 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
445 lines
18 KiB
Python
445 lines
18 KiB
Python
"""E06 암반 비탈면 — Markland 판정 · Hoek-Bray 평면 · 쐐기 안전율 (요구 345 · 351 · PLAN 536-1).
|
||
|
||
근거 = 지식DB `original/안정성평가/03`(Markland 부등식 · Hoek-Bray 평면 · 쐐기) · 06(상용 SW) ·
|
||
Rocscience RocPlane · SWedge 이론서(힘 벡터 · 미끄럼 방향 8 경우 · 수직력 · 지지력은 수동력).
|
||
코드는 직접 짬(Dips · RocPlane · SWedge 코드 안 씀) · NumPy · SciPy 만.
|
||
|
||
좌표 = 동(+x) · 북(+y) · 위(+z). 면 = (경사방향 α · 경사 ψ) · 선 = (추향 · 선주각).
|
||
단위 = 길이 m · 단위중량 kN/m³ · 점착력 kPa · 각 ° · 힘 kN(평면 = 단위 폭당 kN/m).
|
||
|
||
⚠ 지식DB 03 의 평면 볼트 항(수직 T cos(θ+ψp) · 활동 −T sin(θ+ψp))은 θ = 볼트와 면 법선 사이 각일
|
||
때 맞음 — 여기는 θ = 볼트 경사(수평 아래)로 받아 힘 벡터로 풂(수직 T sin(θ+ψp) · 활동 −T cos(θ+ψp)).
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
from typing import Any
|
||
|
||
import numpy as np
|
||
from scipy.optimize import minimize_scalar
|
||
|
||
GAMMA_W = 9.81 # 물 단위중량 kN/m³ (E03 Types 와 같음)
|
||
LATERAL_PLANAR_DEG = 20.0 # 평면 측방 한계(Dips 기본 · Wyllie & Mah)
|
||
LATERAL_TOPPLING_DEG = 20.0 # 전도 측방 한계(Dips 기본 · 지식DB 03 은 30 — 화면에서 바꿈)
|
||
|
||
|
||
# ---- 방향 ----------------------------------------------------------------------
|
||
def ang_diff(a: float, b: float) -> float:
|
||
"""두 방위각 사이 최소 각거리 0~180."""
|
||
d = abs((a - b) % 360.0)
|
||
return min(d, 360.0 - d)
|
||
|
||
|
||
def plane_normal(dip_dir: float, dip: float) -> np.ndarray:
|
||
"""면의 위쪽 단위 법선."""
|
||
a, p = math.radians(dip_dir), math.radians(dip)
|
||
return np.array([math.sin(a) * math.sin(p), math.cos(a) * math.sin(p), math.cos(p)])
|
||
|
||
|
||
def line_vector(trend: float, plunge: float) -> np.ndarray:
|
||
"""선(추향 · 선주각) → 아래로 향한 단위 벡터."""
|
||
t, p = math.radians(trend), math.radians(plunge)
|
||
return np.array([math.sin(t) * math.cos(p), math.cos(t) * math.cos(p), -math.sin(p)])
|
||
|
||
|
||
def vector_line(v: np.ndarray) -> tuple[float, float]:
|
||
"""벡터 → (추향 · 선주각) 하반구(위로 향하면 뒤집음)."""
|
||
v = np.asarray(v, float)
|
||
if v[2] > 0:
|
||
v = -v
|
||
v = v / np.linalg.norm(v)
|
||
trend = math.degrees(math.atan2(v[0], v[1])) % 360.0
|
||
return trend, math.degrees(math.asin(max(-1.0, min(1.0, -v[2]))))
|
||
|
||
|
||
def pole(dip_dir: float, dip: float) -> tuple[float, float]:
|
||
"""면의 극점(하반구) = (α + 180 · 90 − ψ)."""
|
||
return (dip_dir + 180.0) % 360.0, 90.0 - dip
|
||
|
||
|
||
def intersection(a: dict, b: dict) -> tuple[float, float] | None:
|
||
"""두 면 교선(추향 · 선주각) — 나란하면 None."""
|
||
v = np.cross(
|
||
plane_normal(a["dip_direction"], a["dip"]), plane_normal(b["dip_direction"], b["dip"])
|
||
)
|
||
if np.linalg.norm(v) < 1e-9:
|
||
return None
|
||
return vector_line(v)
|
||
|
||
|
||
def apparent_dip(dip: float, dip_dir: float, trend: float) -> float:
|
||
"""면 (α · ψ) 의 trend 방향 겉보기 경사 — 반대쪽(90° 넘음)이면 음수."""
|
||
return math.degrees(
|
||
math.atan(math.tan(math.radians(dip)) * math.cos(math.radians(trend - dip_dir)))
|
||
)
|
||
|
||
|
||
def strike_to_dip_direction(strike: float) -> float:
|
||
"""오른손 법칙 주향 → 경사방향."""
|
||
return (strike + 90.0) % 360.0
|
||
|
||
|
||
# ---- Markland 판정 -------------------------------------------------------------
|
||
def kinematic(
|
||
slope: dict,
|
||
joints: list[dict],
|
||
lateral_planar: float = LATERAL_PLANAR_DEG,
|
||
lateral_toppling: float = LATERAL_TOPPLING_DEG,
|
||
) -> dict[str, list]:
|
||
"""slope = {dip_direction, dip} · joints = [{id, dip_direction, dip, phi}, …].
|
||
|
||
평면 = Δ(αf, αp) ≤ 측방 · φ < ψp < 겉보기 사면 경사 / 쐐기 = 교선 Δ < 90 · φ < 선주각 <
|
||
겉보기 사면 경사(φ = 두 면 중 작은 값) /
|
||
전도 = Δ(αf + 180, αp) ≤ 측방 · ψp > 90 − ψf + φ (Goodman).
|
||
"""
|
||
af, pf = float(slope["dip_direction"]), float(slope["dip"])
|
||
planar, toppling, wedge = [], [], []
|
||
for j in joints:
|
||
aj, pj, phi = float(j["dip_direction"]), float(j["dip"]), float(j["phi"])
|
||
if ang_diff(af, aj) <= lateral_planar and phi < pj < apparent_dip(pf, af, aj):
|
||
planar.append(j["id"])
|
||
if ang_diff(af + 180.0, aj) <= lateral_toppling and pj > 90.0 - pf + phi:
|
||
toppling.append(j["id"])
|
||
for i, a in enumerate(joints):
|
||
for b in joints[i + 1 :]:
|
||
line = intersection(a, b)
|
||
if line is None:
|
||
continue
|
||
trend, plunge = line
|
||
phi = min(float(a["phi"]), float(b["phi"]))
|
||
if ang_diff(af, trend) < 90.0 and phi < plunge < apparent_dip(pf, af, trend):
|
||
wedge.append(
|
||
{
|
||
"ids": [a["id"], b["id"]],
|
||
"trend": round(trend, 2),
|
||
"plunge": round(plunge, 2),
|
||
}
|
||
)
|
||
return {"planar": planar, "wedge": wedge, "toppling": toppling}
|
||
|
||
|
||
# ---- 평면 파괴 (Hoek-Bray · 2D 단위 폭) -----------------------------------------
|
||
def _plane_block(H: float, pf: float, pp: float, ps: float, crack_x: float | None):
|
||
"""사면 끝(0,0) · x = 산 쪽 수평 · z = 위. 덩이 다각형 · 활동면 끝 x · 균열 깊이.
|
||
|
||
윗면 = 마루(H cotψf, H)에서 ψs 로 오름. 균열 없음 = 활동면이 윗면과 만나는 점까지."""
|
||
tf, tp, ts = (math.tan(math.radians(v)) for v in (pf, pp, ps))
|
||
xc = H / tf # 마루
|
||
|
||
def top(x: float) -> float:
|
||
return x * tf if x <= xc else H + (x - xc) * ts
|
||
|
||
if crack_x is None:
|
||
if ts >= tp - 1e-12:
|
||
return None # 활동면이 윗면과 안 만남 = 덩이 안 생김
|
||
xe = (H - xc * ts) / (tp - ts)
|
||
poly = [(0.0, 0.0), (xc, H), (xe, xe * tp)]
|
||
return poly, xe, 0.0
|
||
if crack_x < xc or top(crack_x) - crack_x * tp <= 1e-9:
|
||
return None
|
||
poly = [(0.0, 0.0), (xc, H), (crack_x, top(crack_x)), (crack_x, crack_x * tp)]
|
||
return poly, crack_x, top(crack_x) - crack_x * tp
|
||
|
||
|
||
def _area(poly) -> float:
|
||
xs, zs = np.array(poly).T
|
||
return 0.5 * abs(float(np.dot(xs, np.roll(zs, -1)) - np.dot(zs, np.roll(xs, -1))))
|
||
|
||
|
||
def plane_fs(
|
||
*,
|
||
H: float,
|
||
slope_dip: float,
|
||
plane_dip: float,
|
||
c: float,
|
||
phi: float,
|
||
gamma: float,
|
||
upper_dip: float = 0.0,
|
||
crack: dict | None = None,
|
||
water_fill: float = 0.0,
|
||
kh: float = 0.0,
|
||
bolt_t: float = 0.0,
|
||
bolt_plunge: float = 0.0,
|
||
gamma_w: float = GAMMA_W,
|
||
) -> dict[str, Any] | None:
|
||
"""평면 파괴 안전율 — 덩이가 안 생기면 None.
|
||
|
||
crack = None | {"mode": "depth", "z": m} | {"mode": "critical"}(안전율 최소 자리 찾음) ·
|
||
물 = 균열 있음: zw = fill·z · U = ½γw zw A · V = ½γw zw² / 균열 없음: 가운데 최대 수두
|
||
fill·h/2(h = 활동면 오름) · U = ½γw (fill·h/2) A (Hoek · RocPlane) ·
|
||
kh = 바깥 수평 · 볼트 = 수동력(산 쪽 · 수평 아래 bolt_plunge°)."""
|
||
if not 0 < plane_dip < slope_dip < 90 or H <= 0 or upper_dip >= slope_dip:
|
||
return None
|
||
mode = (crack or {}).get("mode", "none")
|
||
if mode == "critical":
|
||
return _critical_crack(locals())
|
||
tp = math.tan(math.radians(plane_dip))
|
||
xc = H / math.tan(math.radians(slope_dip))
|
||
crack_x = None
|
||
if mode == "depth":
|
||
z = float(crack["z"])
|
||
crack_x = _crack_x_for_depth(H, slope_dip, plane_dip, upper_dip, z)
|
||
if crack_x is None:
|
||
return None
|
||
elif mode == "x":
|
||
crack_x = xc + float(crack["b"])
|
||
block = _plane_block(H, slope_dip, plane_dip, upper_dip, crack_x)
|
||
if block is None:
|
||
return None
|
||
poly, x_end, z = block
|
||
psi = math.radians(plane_dip)
|
||
A = x_end / math.cos(psi)
|
||
W = gamma * _area(poly)
|
||
if crack_x is not None:
|
||
zw = water_fill * z
|
||
U, V = 0.5 * gamma_w * zw * A, 0.5 * gamma_w * zw * zw
|
||
else:
|
||
zw = water_fill * (x_end * tp) / 2.0
|
||
U, V = 0.5 * gamma_w * zw * A, 0.0
|
||
d = np.array([-math.cos(psi), -math.sin(psi)]) # 아래로 미끄러짐(바깥 = −x)
|
||
n = np.array([-math.sin(psi), math.cos(psi)]) # 면에서 덩이 쪽
|
||
active = np.array([-kh * W - V, -W]) + U * n
|
||
drive = float(active @ d)
|
||
if drive <= 0:
|
||
return None
|
||
bolt = bolt_t * np.array(
|
||
[math.cos(math.radians(bolt_plunge)), -math.sin(math.radians(bolt_plunge))]
|
||
)
|
||
N = max(0.0, -float((active + bolt) @ n))
|
||
resist = c * A + N * math.tan(math.radians(phi)) - float(bolt @ d)
|
||
return {
|
||
"F": resist / drive,
|
||
"W": W,
|
||
"A": A,
|
||
"U": U,
|
||
"V": V,
|
||
"z": z,
|
||
"zw": zw,
|
||
"crack_b": None if crack_x is None else crack_x - xc,
|
||
"N": N,
|
||
"drive": drive,
|
||
"resist": resist,
|
||
"polygon": [[round(x, 3), round(y, 3)] for x, y in poly],
|
||
}
|
||
|
||
|
||
def _crack_x_for_depth(H, pf, pp, ps, z) -> float | None:
|
||
"""윗면 아래 균열 깊이 z 가 되는 자리(마루 뒤) — 깊이가 마루 틈보다 크면 None."""
|
||
tf, tp, ts = (math.tan(math.radians(v)) for v in (pf, pp, ps))
|
||
xc = H / tf
|
||
gap = H - xc * tp # 마루에서 윗면 − 활동면
|
||
if z <= 0 or z >= gap:
|
||
return None
|
||
if ts >= tp: # 틈이 커지기만 함 = 뒤로 갈수록 깊어짐
|
||
return xc + (z - gap) / (ts - tp) if ts > tp else None
|
||
return xc + (gap - z) / (tp - ts)
|
||
|
||
|
||
def plane_critical_depth_hoek(H: float, slope_dip: float, plane_dip: float) -> float:
|
||
"""Hoek-Bray 임계 인장균열 깊이(윗면 수평 · 건조) z = H(1 − √(cotψf tanψp))."""
|
||
return H * (
|
||
1.0 - math.sqrt(math.tan(math.radians(plane_dip)) / math.tan(math.radians(slope_dip)))
|
||
)
|
||
|
||
|
||
def _critical_crack(kw: dict) -> dict | None:
|
||
"""균열 자리를 마루 뒤로 옮겨 안전율 최소 — 균열 없음보다 크면 균열 없음."""
|
||
args = {k: kw[k] for k in kw if k not in ("crack", "mode")}
|
||
tf, tp, ts = (math.tan(math.radians(kw[k])) for k in ("slope_dip", "plane_dip", "upper_dip"))
|
||
xc = kw["H"] / tf
|
||
x_far = (kw["H"] - xc * ts) / (tp - ts) if ts < tp else xc + 10 * kw["H"]
|
||
span = x_far - xc
|
||
|
||
def fs_at(b: float) -> float:
|
||
r = plane_fs(**args, crack={"mode": "x", "b": b})
|
||
return r["F"] if r else math.inf
|
||
|
||
best = minimize_scalar(fs_at, bounds=(span * 1e-3, span * (1 - 1e-3)), method="bounded")
|
||
found = plane_fs(**args, crack={"mode": "x", "b": float(best.x)})
|
||
plain = plane_fs(**args, crack=None)
|
||
if found is None or (plain is not None and plain["F"] <= found["F"]):
|
||
return plain
|
||
return found
|
||
|
||
|
||
# ---- 쐐기 파괴 (Hoek-Bray · SWedge 꼴 사면체) ------------------------------------
|
||
def _meet(planes: list[tuple[np.ndarray, float]]) -> np.ndarray | None:
|
||
"""세 면 n·x = d 의 교점."""
|
||
m = np.array([p[0] for p in planes])
|
||
if abs(np.linalg.det(m)) < 1e-10:
|
||
return None
|
||
return np.linalg.solve(m, np.array([p[1] for p in planes]))
|
||
|
||
|
||
def _tri_area(a, b, c) -> float:
|
||
return 0.5 * float(np.linalg.norm(np.cross(b - a, c - a)))
|
||
|
||
|
||
def wedge_geometry(j1: dict, j2: dict, slope: dict, upper: dict, H: float) -> dict | None:
|
||
"""사면체 쐐기 — 꼭짓점 O(두 절리 · 사면 = 사면 끝) · A(1 · 사면 · 윗면) · B(2 · 사면 · 윗면) ·
|
||
C(1 · 2 · 윗면). 사면 높이 H = 사면 위 꼭짓점 A · B 중 높은 쪽의 O 위 높이(SWedge 와 같음 —
|
||
Priest 예제 부피 45.201 m³ 재현). 교선이 사면 · 윗면으로 안 나오면 None(쐐기 안 생김)."""
|
||
n1 = plane_normal(j1["dip_direction"], j1["dip"])
|
||
n2 = plane_normal(j2["dip_direction"], j2["dip"])
|
||
nf = plane_normal(slope["dip_direction"], slope["dip"])
|
||
nu = plane_normal(upper["dip_direction"], upper["dip"])
|
||
# 윗면 = 사면 위로 H 오른 점을 지남(사면 경사방향 반대로 H cotψf 들어간 점)
|
||
back = -line_vector(slope["dip_direction"], 0.0) * (H / math.tan(math.radians(slope["dip"])))
|
||
crest = back + np.array([0.0, 0.0, H])
|
||
du = float(nu @ crest)
|
||
O = np.zeros(3) # noqa: E741 — 꼭짓점 이름(SWedge 그림과 같음)
|
||
A = _meet([(n1, 0.0), (nf, 0.0), (nu, du)])
|
||
B = _meet([(n2, 0.0), (nf, 0.0), (nu, du)])
|
||
C = _meet([(n1, 0.0), (n2, 0.0), (nu, du)])
|
||
if A is None or B is None or C is None:
|
||
return None
|
||
# C 는 사면 뒤(산 쪽) · 위 — 교선이 사면으로 나와 윗면에 닿음
|
||
if C[2] <= 1e-9 or float(nf @ C) >= -1e-9 or A[2] <= 1e-9 or B[2] <= 1e-9:
|
||
return None
|
||
scale = H / max(A[2], B[2]) # 면이 모두 O 를 지나 닮은꼴 — 높이만 맞춤
|
||
A, B, C = A * scale, B * scale, C * scale
|
||
volume = abs(float(np.linalg.det(np.array([A - O, B - O, C - O])))) / 6.0
|
||
if volume <= 1e-9:
|
||
return None
|
||
centroid = (O + A + B + C) / 4.0
|
||
inward = [n if float(n @ centroid) > 0 else -n for n in (n1, n2)]
|
||
return {
|
||
"vertices": {"O": O, "A": A, "B": B, "C": C},
|
||
"volume": volume,
|
||
"areas": [_tri_area(O, A, C), _tri_area(O, B, C)],
|
||
"face_area": _tri_area(O, A, B),
|
||
"upper_area": _tri_area(A, B, C),
|
||
"normals": inward,
|
||
"Hw": float(C[2]),
|
||
"line_length": float(np.linalg.norm(C)),
|
||
}
|
||
|
||
|
||
def _unit(v: np.ndarray) -> np.ndarray:
|
||
n = float(np.linalg.norm(v))
|
||
return v / n if n > 1e-12 else v * 0.0
|
||
|
||
|
||
def _sliding(A: np.ndarray, n1: np.ndarray, n2: np.ndarray, W: np.ndarray):
|
||
"""SWedge 미끄럼 방향 — (갈래, 방향). 갈래 = lift · 1 · 2 · 12 · stable."""
|
||
s1 = _unit(np.cross(np.cross(n1, A), n1))
|
||
s2 = _unit(np.cross(np.cross(n2, A), n2))
|
||
c12 = np.cross(n1, n2)
|
||
s12 = _unit(c12) * np.sign(float(c12 @ A))
|
||
if A @ n1 > 0 and A @ n2 > 0:
|
||
return "lift", _unit(A)
|
||
if A @ n1 <= 0 and s1 @ n2 > 0:
|
||
return "1", s1
|
||
if A @ n2 <= 0 and s2 @ n1 > 0:
|
||
return "2", s2
|
||
if s1 @ n2 <= 0 and s2 @ n1 <= 0:
|
||
return "12", s12
|
||
return "stable", None
|
||
|
||
|
||
def _normals(mode: str, F: np.ndarray, n1: np.ndarray, n2: np.ndarray) -> tuple[float, float]:
|
||
if mode == "1":
|
||
return -float(F @ n1), 0.0
|
||
if mode == "2":
|
||
return 0.0, -float(F @ n2)
|
||
if mode == "12":
|
||
c12, c21 = np.cross(n1, n2), np.cross(n2, n1)
|
||
return (
|
||
-float(np.cross(F, n2) @ c12) / float(c12 @ c12),
|
||
-float(np.cross(F, n1) @ c21) / float(c21 @ c21),
|
||
)
|
||
return 0.0, 0.0
|
||
|
||
|
||
def wedge_fs(
|
||
*,
|
||
j1: dict,
|
||
j2: dict,
|
||
slope: dict,
|
||
upper: dict,
|
||
H: float,
|
||
gamma: float,
|
||
water: dict | None = None,
|
||
kh: float = 0.0,
|
||
seismic_trend: float | None = None,
|
||
external: list[dict] | None = None,
|
||
bolt: dict | None = None,
|
||
gamma_w: float = GAMMA_W,
|
||
) -> dict[str, Any] | None:
|
||
"""쐐기 안전율(SWedge 이론서 꼴) — 쐐기가 안 생기면 None.
|
||
|
||
j = {dip_direction, dip, c, phi} · water = {"fill": p}(중력식 U = p³γw Hw a/6) |
|
||
{"pressure": [u1, u2]}(평균 수압 kPa) · kh = 수평(기본 교선 추향) · external = [{force, trend,
|
||
plunge}] · bolt = {force, trend, plunge}(수동력). FS = max(지지 전 · 지지 뒤 · 들림)."""
|
||
geo = wedge_geometry(j1, j2, slope, upper, H)
|
||
if geo is None:
|
||
return None
|
||
n1, n2 = geo["normals"]
|
||
a1, a2 = geo["areas"]
|
||
weight = gamma * geo["volume"]
|
||
Wv = np.array([0.0, 0.0, -weight])
|
||
active = Wv.copy()
|
||
water = water or {}
|
||
if "pressure" in water:
|
||
u1, u2 = (float(u) for u in water["pressure"])
|
||
else:
|
||
p = float(water.get("fill", 0.0))
|
||
u1 = u2 = p**3 * gamma_w * geo["Hw"] / 6.0
|
||
Uf = [u1 * a1 * n1, u2 * a2 * n2]
|
||
active += Uf[0] + Uf[1]
|
||
line = vector_line(np.cross(n1, n2))
|
||
if kh:
|
||
trend = line[0] if seismic_trend is None else seismic_trend
|
||
active += kh * weight * line_vector(trend, 0.0)
|
||
for item in external or []:
|
||
active += float(item["force"]) * line_vector(item["trend"], item["plunge"])
|
||
passive = np.zeros(3)
|
||
if bolt and bolt.get("force"):
|
||
passive = float(bolt["force"]) * line_vector(bolt["trend"], bolt["plunge"])
|
||
mode, s = _sliding(active, n1, n2, Wv)
|
||
out = {
|
||
"volume": geo["volume"],
|
||
"W": weight,
|
||
"areas": [a1, a2],
|
||
"face_area": geo["face_area"],
|
||
"upper_area": geo["upper_area"],
|
||
"Hw": geo["Hw"],
|
||
"line_length": geo["line_length"],
|
||
"line": [round(line[0], 2), round(line[1], 2)],
|
||
"water_force": [float(np.linalg.norm(Uf[0])), float(np.linalg.norm(Uf[1]))],
|
||
"mode": mode,
|
||
"vertices": {k: [round(float(x), 3) for x in v] for k, v in geo["vertices"].items()},
|
||
}
|
||
if mode == "stable":
|
||
return {**out, "F": math.inf, "N": [0.0, 0.0], "drive": 0.0, "resist": 0.0}
|
||
drive = float(active @ s)
|
||
if mode == "lift":
|
||
fs = -float(passive @ s) / drive if drive > 0 else math.inf
|
||
return {**out, "F": fs, "N": [0.0, 0.0], "drive": drive, "resist": -float(passive @ s)}
|
||
joints = (j1, j2)
|
||
|
||
def resisting(F: np.ndarray) -> tuple[float, list[float]]:
|
||
N = _normals(mode, F, n1, n2)
|
||
total = 0.0
|
||
for k, (n_k, a_k) in enumerate(zip(N, (a1, a2))):
|
||
if str(k + 1) not in mode or n_k <= 0:
|
||
continue # 떨어진 면 = 저항 없음(SWedge 와 같음)
|
||
tau = float(joints[k]["c"]) + n_k / a_k * math.tan(
|
||
math.radians(float(joints[k]["phi"]))
|
||
)
|
||
total += tau * a_k
|
||
return total, [max(0.0, n) for n in N]
|
||
|
||
j_u, n_u = resisting(active)
|
||
fs_u = j_u / drive if drive > 0 else math.inf
|
||
fs, n_used, resist = fs_u, n_u, j_u
|
||
if np.any(passive):
|
||
j_s, n_s = resisting(active + passive)
|
||
fs_s = (-float(passive @ s) + j_s) / drive if drive > 0 else math.inf
|
||
fs_lift = -float(passive @ _unit(active)) / float(np.linalg.norm(active))
|
||
if fs_s > fs:
|
||
fs, n_used, resist = fs_s, n_s, -float(passive @ s) + j_s
|
||
fs = max(fs, fs_lift)
|
||
return {**out, "F": fs, "N": n_used, "drive": drive, "resist": resist}
|