262 lines
12 KiB
Python
262 lines
12 KiB
Python
|
# -*- coding: utf-8 -*-
|
||
|
|
"""Structural check of the printed arm, using section properties measured from the SOLIDS.
|
||
|
|
|
||
|
|
Nothing here is hand-waved from a sketch dimension: for each critical cut plane the part is
|
||
|
|
boolean-cut and the freshly created planar face is taken and integrated with
|
||
|
|
BRepGProp.SurfaceProperties_s, so area, centroid and the second moments are the real
|
||
|
|
section of the real printed part (including whatever the lightening windows left behind).
|
||
|
|
|
||
|
|
Why the allowables are so low: a printed part is not the material datasheet. Bending
|
||
|
|
across layers, an un-filleted internal corner (Kt ~ 2.5 measured on square corners), and
|
||
|
|
PLA creep under permanent static load together cut 35-50 MPa down to single digits, which
|
||
|
|
is exactly the failure mode to design against. See docs/strength.md.
|
||
|
|
|
||
|
|
Run: python strength.py -> docs/strength.md (+ console table)
|
||
|
|
"""
|
||
|
|
import io
|
||
|
|
import json
|
||
|
|
import math
|
||
|
|
import os
|
||
|
|
import sys
|
||
|
|
|
||
|
|
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||
|
|
try:
|
||
|
|
sys.stdout.reconfigure(encoding="utf-8")
|
||
|
|
except Exception:
|
||
|
|
pass
|
||
|
|
import cadquery as cq
|
||
|
|
from armcore import OUT, DOCS, solids_of, bb, mesh_of
|
||
|
|
from OCP.BRepAlgoAPI import BRepAlgoAPI_Common
|
||
|
|
from OCP.BRepGProp import BRepGProp
|
||
|
|
from OCP.GProp import GProp_GProps
|
||
|
|
from OCP.BRepPrimAPI import BRepPrimAPI_MakeBox
|
||
|
|
from OCP.gp import gp_Pnt, gp_Dir, gp_Ax2, gp_Vec
|
||
|
|
|
||
|
|
import human_layout as HL
|
||
|
|
|
||
|
|
RHO = {"PLA": 1.24e-3, "motor_g": 210.0} # g/mm3 for printed, g per module
|
||
|
|
# a slicer does not print solid: 3 perimeters + 30% infill lands near this fraction of
|
||
|
|
# the CAD envelope, and using the solid mass would blame the part for plastic that is not there
|
||
|
|
INFILL = float(os.environ.get("ARMINFILL", "0.38"))
|
||
|
|
# allowable bending stress across the print layers, after notch + creep derating
|
||
|
|
ALLOW = {"PLA": 12.0, "PETG": 20.0, "PA_CF": 30.0} # MPa
|
||
|
|
MAT = os.environ.get("ARMMAT", "PLA")
|
||
|
|
PAYLOAD = float(os.environ.get("ARMPAY", "0.5")) # kg held at the fingertip
|
||
|
|
DYN = float(os.environ.get("ARMDYN", "2.0")) # dynamic factor on the moving mass
|
||
|
|
|
||
|
|
G981 = 9.81
|
||
|
|
GRIP_T = 12.0 # N*mm per motor is far too small to matter; the grip force
|
||
|
|
# case is carried separately below as a normal force at the pad
|
||
|
|
_cache = {}
|
||
|
|
|
||
|
|
|
||
|
|
def load(name):
|
||
|
|
if name not in _cache:
|
||
|
|
_cache[name] = cq.importers.importStep("%s/parts_step/%s.step" % (OUT, name)).val()
|
||
|
|
return _cache[name]
|
||
|
|
|
||
|
|
|
||
|
|
def volume(sh):
|
||
|
|
p = GProp_GProps()
|
||
|
|
BRepGProp.VolumeProperties_s(sh.wrapped, p)
|
||
|
|
return p.Mass()
|
||
|
|
|
||
|
|
|
||
|
|
def _caps_at(sh, n, o, t):
|
||
|
|
"""The section faces of (solid ∩ thin slab), taken on the +normal side of the slab.
|
||
|
|
|
||
|
|
Two things this replaces. (1) 'Zero extent along the normal' on a HALF-SPACE cut also
|
||
|
|
matches the part's own horizontal faces (top chord, gusset crowns) sitting anywhere above
|
||
|
|
the plane - that over-read the shoulder fork by 2.4x (480 mm2 true vs 1155 reported).
|
||
|
|
(2) The half-space Common outright fails on the filleted fork at z=600 (returns an empty
|
||
|
|
shape) while a thin slab of the same solid cuts fine.
|
||
|
|
A slab is also self-checking: its volume / t must equal the integrated area.
|
||
|
|
"""
|
||
|
|
c0 = sum(n[k] * o[k] for k in range(3))
|
||
|
|
lo, hi = [], []
|
||
|
|
for k in range(3):
|
||
|
|
if n[k] > 0:
|
||
|
|
lo.append(o[k] - t / 2.0); hi.append(o[k] + t / 2.0)
|
||
|
|
elif n[k] < 0:
|
||
|
|
lo.append(o[k] - t / 2.0); hi.append(o[k] + t / 2.0)
|
||
|
|
else:
|
||
|
|
lo.append(o[k] - 400.0); hi.append(o[k] + 400.0)
|
||
|
|
box = BRepPrimAPI_MakeBox(gp_Pnt(*lo), gp_Pnt(*hi)).Shape()
|
||
|
|
cut = BRepAlgoAPI_Common(sh.wrapped, box)
|
||
|
|
cut.Build()
|
||
|
|
if not cut.IsDone():
|
||
|
|
return []
|
||
|
|
off = c0 + t / 2.0
|
||
|
|
out = []
|
||
|
|
for f in cq.Shape(cut.Shape()).Faces():
|
||
|
|
if f.geomType() != "PLANE" or f.Area() <= 1e-6:
|
||
|
|
continue
|
||
|
|
cx, cy, cz = f.Center().toTuple()
|
||
|
|
if abs(n[0] * cx + n[1] * cy + n[2] * cz - off) > 1e-4:
|
||
|
|
continue
|
||
|
|
out.append(f)
|
||
|
|
return out
|
||
|
|
|
||
|
|
|
||
|
|
def section_props(sh, normal, origin, bend=2, size=400.0):
|
||
|
|
"""Cut the solid with a thin slab and integrate the section face.
|
||
|
|
|
||
|
|
Returns (area mm2, Ixx mm4, Iyy mm4, cmax mm, bbox) about the face centroid, with the
|
||
|
|
in-plane axes taken as (u, v) where u x v = normal.
|
||
|
|
"""
|
||
|
|
n = [float(x) for x in normal]
|
||
|
|
o = [float(x) for x in origin]
|
||
|
|
caps = []
|
||
|
|
for t in (0.2, 0.05, 1.0):
|
||
|
|
for nudge in (0.0, 0.05, -0.05, 0.25, -0.25):
|
||
|
|
oo = [o[k] + nudge * n[k] for k in range(3)]
|
||
|
|
caps = _caps_at(sh, n, oo, t)
|
||
|
|
if caps:
|
||
|
|
break
|
||
|
|
if caps:
|
||
|
|
break
|
||
|
|
if not caps:
|
||
|
|
return None
|
||
|
|
plane_ax = [k for k in (0, 1, 2) if abs(n[k]) < 0.5]
|
||
|
|
p_ax, q_ax = plane_ax
|
||
|
|
bend_ax = bend if bend in plane_ax else q_ax
|
||
|
|
far_ax = p_ax if bend_ax == q_ax else q_ax
|
||
|
|
A = Sv = Suu = Svv = S = 0.0
|
||
|
|
for f in caps:
|
||
|
|
verts, tris = mesh_of(f, 0.02)
|
||
|
|
fa = fu = fv = fuu = fvv = 0.0
|
||
|
|
for t in tris:
|
||
|
|
(u1, v1), (u2, v2), (u3, v3) = [(verts[k][p_ax], verts[k][q_ax]) for k in t]
|
||
|
|
ar = ((u2 - u1) * (v3 - v1) - (u3 - u1) * (v2 - v1)) / 2.0
|
||
|
|
fa += ar # SIGNED: inner wires (holes) cancel
|
||
|
|
fu += ar * (u1 + u2 + u3) / 3.0
|
||
|
|
fv += ar * (v1 + v2 + v3) / 3.0
|
||
|
|
fvv += ar * (v1**2 + v2**2 + v3**2 + v1*v2 + v2*v3 + v3*v1) / 6.0
|
||
|
|
fuu += ar * (u1**2 + u2**2 + u3**2 + u1*u2 + u2*u3 + u3*u1) / 6.0
|
||
|
|
if fa < 0: # face normal opposed to the cut normal
|
||
|
|
fa, fu, fv, fvv, fuu = -fa, -fu, -fv, -fvv, -fuu
|
||
|
|
A += fa; S += fu; Sv += fv; Svv += fvv; Suu += fuu
|
||
|
|
if A <= 1e-9:
|
||
|
|
return None
|
||
|
|
uc, vc = S / A, Sv / A
|
||
|
|
best = max(caps, key=lambda f: f.Area()) # only for the bounding box below
|
||
|
|
allb = [f.BoundingBox() for f in caps]
|
||
|
|
lo = [min(getattr(b, ("xmin", "ymin", "zmin")[k]) for b in allb) for k in range(3)]
|
||
|
|
hi = [max(getattr(b, ("xmax", "ymax", "zmax")[k]) for b in allb) for k in range(3)]
|
||
|
|
cmax = max(hi[far_ax] - uc if far_ax == p_ax else hi[far_ax] - vc,
|
||
|
|
uc - lo[far_ax] if far_ax == p_ax else vc - lo[far_ax]) / 1.0
|
||
|
|
# distance from the neutral axis to the extreme fibre, measured in the bending direction
|
||
|
|
ext = (hi[far_ax] - lo[far_ax])
|
||
|
|
cen = uc if far_ax == p_ax else vc
|
||
|
|
cmax = max(hi[far_ax] - cen, cen - lo[far_ax])
|
||
|
|
I = (Svv - A * vc * vc) if far_ax == q_ax else (Suu - A * uc * uc)
|
||
|
|
return dict(A=A, I=I, c=cmax, bend=bend_ax, centroid=(uc, vc))
|
||
|
|
|
||
|
|
|
||
|
|
def arm_masses():
|
||
|
|
"""Everything distal to each joint, in grams, from the model itself."""
|
||
|
|
meta = json.load(io.open(OUT.replace("/out", "/web") + "/mesh_arm.json", encoding="utf-8"))
|
||
|
|
TREE = ["root", "J1", "J2", "J3", "J4", "J5", "J6"]
|
||
|
|
tot = {}
|
||
|
|
for p in meta["parts"]:
|
||
|
|
sh = load(p["name"])
|
||
|
|
v = volume(sh) / 1000.0
|
||
|
|
g = (v * RHO["PLA"] * 1000.0 * INFILL) if p["printed"] else RHO["motor_g"]
|
||
|
|
tot[p["name"]] = (g, p["grp"])
|
||
|
|
out = {}
|
||
|
|
for k, jid in enumerate(TREE[1:], 1):
|
||
|
|
down = set(TREE[k:])
|
||
|
|
out[jid] = sum(g for g, grp in tot.values() if grp in down)
|
||
|
|
out["all"] = sum(g for g, _ in tot.values())
|
||
|
|
return out, tot
|
||
|
|
|
||
|
|
|
||
|
|
# critical cuts: (label, part, plane normal, plane origin, which load case owns it)
|
||
|
|
def cuts():
|
||
|
|
"""(label, part, cut-plane normal, point, joint whose load it carries, bend axis)
|
||
|
|
bend axis: 0=X (roll of the shoulder/forearm), 1=Y (every pitch joint), 2=Z."""
|
||
|
|
Y = HL.Y_AXIS
|
||
|
|
return [
|
||
|
|
("吊臂根座 @z=480", "P03_肩吊臂", (0, 0, 1), (0.0, 55.0, 480.0), "J1", 1),
|
||
|
|
("吊臂颊板 @z=560", "P03_肩吊臂", (0, 0, 1), (0.0, 100.0, 560.0), "J1", 1),
|
||
|
|
("肩叉腹板 @z=560", "P04_肩叉", (0, 0, 1), (HL.X_ARM, Y, 560.0), "J2", 1),
|
||
|
|
("肩叉腹板 @z=600", "P04_肩叉", (0, 0, 1), (HL.X_ARM, Y, 600.0), "J2", 1),
|
||
|
|
("大臂旋转座横梁 @z=480", "P05_大臂旋转座", (0, 0, 1), (HL.X_ARM, Y, 480.0), "J2", 1),
|
||
|
|
("大臂旋转座套筒 @z=440", "P05_大臂旋转座", (0, 0, 1), (HL.X_ARM, Y, 440.0), "J3", 1),
|
||
|
|
("大臂骨架中段 @z=350", "P06_大臂骨架", (0, 0, 1), (HL.X_ARM, Y, 350.0), "J3", 1),
|
||
|
|
("大臂骨架抱箍 @z=325", "P06_大臂骨架", (0, 0, 1), (HL.X_ARM, Y, 325.0), "J4", 1),
|
||
|
|
("小臂骨架中段 @z=245", "P07_小臂骨架", (0, 0, 1), (HL.X_ARM, Y, 245.0), "J4", 1),
|
||
|
|
("小臂骨架抱箍 @z=195", "P07_小臂骨架", (0, 0, 1), (HL.X_ARM, Y, 195.0), "J5", 1),
|
||
|
|
("腕偏转座根部 @z=150", "P08_腕偏转座", (0, 0, 1), (HL.X_ARM, Y, 150.0), "J5", 1),
|
||
|
|
("工具法兰 @z=55", "P09_工具法兰", (0, 0, 1), (HL.X_ARM, Y, 55.0), "J6", 0),
|
||
|
|
]
|
||
|
|
|
||
|
|
|
||
|
|
MOMENT = { # worst reachable bending moment per joint, N*mm
|
||
|
|
# computed in main() from the reachable pose + masses
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
def main():
|
||
|
|
masses, per = arm_masses()
|
||
|
|
# Moment arms DERIVED from the swept travel (docs/travel.json), not hand-typed: they were
|
||
|
|
# stale (J4 said 24 deg while the arm now sweeps 360, and J1 was missing entirely), which
|
||
|
|
# silently under-loaded those sections. A continuous or >90 deg axis is capped at 90 deg
|
||
|
|
# because that is where a gravity load sits furthest off the axis.
|
||
|
|
TRV = json.load(io.open(DOCS + "/travel.json", encoding="utf-8"))
|
||
|
|
dev = {t["axis"]: min(max(abs(t["hi"]), abs(t["lo"])), 90.0) for t in TRV}
|
||
|
|
# distance from each joint axis out to the fingertip load, mm
|
||
|
|
dist = {"J1": 564.0, "J2": 444.0, "J3": 380.0, "J4": 264.0, "J5": 129.0, "J6": 60.0}
|
||
|
|
reach = {k: dist[k] * math.sin(math.radians(dev.get(k, 90.0))) for k in dist}
|
||
|
|
from concurrent.futures import ThreadPoolExecutor
|
||
|
|
|
||
|
|
def compute_cut(c_item):
|
||
|
|
lab, part, n, o, jid, bend = c_item
|
||
|
|
sh = load(part)
|
||
|
|
sp = section_props(sh, n, o, bend)
|
||
|
|
if not sp:
|
||
|
|
return (lab, part, None, jid)
|
||
|
|
m_g = masses[jid] * DYN + PAYLOAD * 1000.0 * DYN
|
||
|
|
F = m_g * G981 / 1000.0 # N
|
||
|
|
# the section carries everything distal to its own joint, at that joint's worst
|
||
|
|
# reachable moment arm, plus the grip reaction on the finger root
|
||
|
|
M = F * reach[jid] # N*mm
|
||
|
|
if jid == "J6":
|
||
|
|
M = max(M, GRIP_T * 60.0)
|
||
|
|
I = sp["I"] # already the real bend axis, not min(Iu,Iv)
|
||
|
|
sigma = M * sp["c"] / max(I, 1e-6)
|
||
|
|
tau = 1.5 * F / sp["A"]
|
||
|
|
return (lab, part, dict(sp, M=M, F=F, sigma=sigma, tau=tau,
|
||
|
|
SF=ALLOW[MAT] / sigma if sigma > 0 else 99.0), jid)
|
||
|
|
|
||
|
|
with ThreadPoolExecutor(max_workers=min(13, os.cpu_count() or 4)) as executor:
|
||
|
|
rows = list(executor.map(compute_cut, cuts()))
|
||
|
|
print("材料 %s 许用弯曲 %g MPa（跨层、含缺口与蠕变折减） 末端负载 %g kg 动载系数 %g\n"
|
||
|
|
% (MAT, ALLOW[MAT], PAYLOAD, DYN))
|
||
|
|
print("打印件按 %d%% 实体率计质量（壁+填充），电机按 %.0f g/台" % (INFILL * 100, RHO["motor_g"]))
|
||
|
|
print("%-26s %-10s %7s %9s %6s %8s %7s %5s %9s" %
|
||
|
|
("截面", "零件", "面积mm2", "Imin mm4", "c mm", "M N*mm", "sigma", "SF", "需Z mm3"))
|
||
|
|
bad = []
|
||
|
|
for lab, part, d, jid in rows:
|
||
|
|
if not d:
|
||
|
|
print("%-26s %-10s 剖切失败" % (lab, part))
|
||
|
|
continue
|
||
|
|
Zreq = d["M"] / ALLOW[MAT] * 1.5 # section modulus x1.5 SF
|
||
|
|
print("%-26s %-10s %7.0f %9.0f %6.1f %8.0f %7.1f %5.2f %9.0f%s" %
|
||
|
|
(lab, part, d["A"], d["I"], d["c"], d["M"], d["sigma"], d["SF"], Zreq,
|
||
|
|
" <<< 不足" if d["SF"] < 1.5 else ""))
|
||
|
|
if d["SF"] < 1.5:
|
||
|
|
bad.append((lab, part, d))
|
||
|
|
print("\n各关节以下质量(g, 含电机): " + " ".join("%s=%.0f" % (k, v) for k, v in masses.items() if k != "all"))
|
||
|
|
print("整机质量(含 6 电机) %.0f g" % masses["all"])
|
||
|
|
io.open(DOCS + "/strength.json", "w", encoding="utf-8").write(
|
||
|
|
json.dumps([dict(截面=a, 零件=b, 关节=d, **(c or {})) for a, b, c, d in rows],
|
||
|
|
ensure_ascii=False, indent=1, default=float))
|
||
|
|
return bad
|
||
|
|
|
||
|
|
|
||
|
|
if __name__ == "__main__":
|
||
|
|
b = main()
|
||
|
|
print("\n不合格截面 %d 处" % len(b))
|
||
|
|
os._exit(0)
|