一个零件受力之后，屏幕上的红色区域到底代表什么？它是位移大、应力大，还是只是色标被调低了？如果模型被放大变形十倍，探针读数是否也应该乘十？这些问题必须沿着 CAD 几何 → 数学模型 → 求解结果 → 表面映射 → GPU 绘图 逐层回答。渲染器只消费结果，可靠的云图必须保留结果的单位、位置和计算条件。
本文先用一个可以独立校核的悬臂零件跑通完整链路，再对照已有的泵支架模态数据。悬臂案例在浏览器中实际执行有限元静力求解；泵支架读取上游预计算结果。两者对应不同的物理问题。
先把零件、载荷和坐标说清楚
设想一条矩形钢制悬臂，长度 L=200 mm，截面沿 Y 高 h=8 mm、沿 Z 宽 b=20 mm。X 轴沿长度，X=0 的端面固定，X=L 的截面中心受向下的集中力 F=−100 N。Y 正方向向上，因此挠度 w 为负。这里的集中力作用于梁模型的自由端节点；实际夹具、螺栓接触与力的分布需要另建模型。
CAD 用 CadQuery 2.6.1 参数化建模。cantilever.json 保存几何契约，建模脚本读取它生成实体与 STEP；求解器与显示表面也读取同一个契约。脚本检查单一实体、几何有效性、体积与 STEP 重导入尺寸，避免“图看起来差不多，但坐标系或单位已经变了”。
# examples/stage-02-engineering/user_code/build_cad.py
# X: 长度；Y: 弯曲高度；Z: 截面宽度。CAD 参数单位为 mm。
part = cq.Workplane("XY").box(200, 8, 20,
centered=(False, True, True))
cq.exporters.export(part, "cantilever.step")
# 单个封闭实体；6 个平面；体积 32000 mm³。本例采用 Euler–Bernoulli 梁模型：截面在弯曲后保持平面且垂直于中性轴，忽略剪切变形。L/h=25，适合作为细长梁的小挠度基准。截面惯性矩 I=b·h³/12=853.333 mm⁴，它描述截面对弯曲的几何抵抗能力；I 不是质量惯性矩。求解全程转换为 SI：m、N、Pa。弹性模量 E=210 GPa 描述钢在弹性阶段的应力与应变比例，和零件几何共同决定 EI 弯曲刚度。
CadQuery 负责参数化实体建模；TU Delft 的梁单元推导给出本文采用的 Hermite 插值与刚度矩阵。尺寸、载荷、材料数值是本站的基准案例设定。
CAD 表面为什么不能直接拿来做有限元
STEP 描述的是连续几何与拓扑，渲染用的是三角表面，求解用的是离散的自由度与单元。三个对象服务的任务不同。画得更细的三角形，不会自动给求解器增加物理自由度；反过来，求解器输出的内部节点也不一定是表面可见顶点。
本例沿 X 把中性轴划分为 n 个梁单元。每个节点有横向位移 w 与截面转角 θ 两个自由度，共 2(n+1) 个。截面厚度通过 I 进入计算，求解器没有给实体内部生成四面体。显示则把 CAD 六个平面按长度方向分成 80 段，保留尖锐边界的独立顶点记录。这是针对矩形实体的结构化表面采样，尺寸由已经重导入验证的 CAD 契约给出；复杂 STEP 的实时三角化见另一篇 CAD 渲染实践。
查看 Mermaid 源码
flowchart TD
Spec[CAD 尺寸契约 / mm] --> CAD[CadQuery 实体与 STEP]
CAD --> Check[重导入 / 尺寸与体积校核]
Spec --> Beam[中心线梁单元 / m]
Beam --> Solve[组装 K 与 f / 施加固定约束]
Solve --> Result[节点 w 与 θ / 反力]
Spec --> Surface[六个 CAD 平面的结构化采样]
Result --> Map[Hermite 插值 / 位移与应力恢复]
Surface --> Map
Map --> GPU[position 与 color 缓冲 / Three.js]对于真实复杂零件，通常需要 CAD 清理、物理区域定义、体网格生成、单元质量检查、网格收敛验证，再从求解网格提取边界。圆角、薄壁和孔洞会改变应力分布，不能为了方便画图擅自删除。展示网格若另行简化，必须保存它与求解网格的插值映射。
从一个梁单元组装到 Ku=f
梁单元的局部未知量是 d=[w1, θ1, w2, θ2]。用三次 Hermite 形函数插值挠度，可以同时维持单元交界处的位移和斜率连续。令单元长度为 ℓ，ξ=x/ℓ，挠度 w(x)=N(x)d。二阶导数给出曲率 κ=w″，截面轴向应变为 εxx=−yκ，应力为 σxx=Eεxx。
刚度矩阵来自弯曲应变能：对单元积分 E·I·(N″)TN″，得到下面的矩阵。θ 是无量纲的弧度，w 是米，所以 K 的子块具有不同量纲；不能简单把整个矩阵统一叫作 N/m。
┌ 12 6ℓ −12 6ℓ ┐
Ke = EI/ℓ³ · │ 6ℓ 4ℓ² −6ℓ 2ℓ² │
│ −12 −6ℓ 12 −6ℓ │
└ 6ℓ 2ℓ² −6ℓ 4ℓ² ┘共享节点的矩阵项累加进全局 K，端部力写入 f 的末端 w 分量。固定端 w0=θ0=0，把对应的两行两列消去，求解剩余的自由自由度。原始完整 K 和 f 仍保留：求出 u 后计算 r=Ku−f，受约束的两个分量是固定端反力与反力矩，自由分量应接近零。
// 正式求解器：src/components/engineering/beam-solver.ts
const result = solveBeam({
length: 0.2, width: 0.02, height: 0.008,
young: 210e9, force: -100, elements: 20,
});
const section = result.sample(0.1); // 输入位置 x，单位 m
// sample 用所在单元的四个已求解自由度作 Hermite 插值。
console.log(result.tip * 1000); // −1.488095… mm
console.log(section.stress(0.004) / 1e6);// 上表面 +46.875 MPa本站小矩阵采用带部分主元选择的消元，计算规模控制在 80 个单元以内。工程中的大型三维模型应使用稀疏装配、适当的直接/迭代求解器与预条件器，并记录求解容差、迭代历史和收敛状态。若漏掉固定约束，整体刚体运动会令刚度矩阵奇异；用一个极小数字替换主元，会隐藏建模错误。
先校核，再看三维结果
这组均匀梁、端部集中力的解析解是 w(x)=F·x²(3L−x)/(6EI)，自由端 w(L)=FL³/(3EI)，弯矩 M(x)=F(L−x)。本例在这一简单载荷下的真实解恰好是三次多项式，因此一个三次梁单元就能表达它。这不能推广为“所有模型一个单元就足够”；分布载荷、变截面、复杂边界都要重新验证。
| 校核量 | 本例 F=−100 N、E=210 GPa | 说明 |
|---|---|---|
| 自由端位移 | −1.488095 mm | 与解析解比较，保留方向 |
| 固定端竖向反力 | +100 N | 和端部外力平衡 |
| 固定端反力矩 | +20 N·m | 和外力关于固定端的力矩平衡 |
| 根部上表面应力 | +93.75 MPa | 拉伸；下表面为 −93.75 MPa |
| 中性轴轴向应力 | 0 MPa | y=0，所以弯曲应变为零 |
下面把载荷从 100 N 改为 200 N：位移与应力都应加倍。把钢换成 70 GPa 的铝，位移变为三倍；在同一几何与静定端力下，弯曲应力保持不变。再改变色标和变形倍率，探针数值应完全相同。