本页目录

前沿 I · 有限元与结构优化

层次：本科接触 + 业界日常工具 + 研究生方向 ｜ 🔗 前置：能量法（第五页）、线性代数与数值方法（数学站）。 有限元法（FEM）把第一至五页的固体力学变成计算机能解的问题。它是机械工程师最常用的分析工具，也是最容易被误用的工具——因为它总能给出漂亮的彩色云图，无论输入是否合理。本页讲原理、讲怎么用对、以及它如何延伸到自动生成结构的拓扑优化。

学习层：先审计一根杆，再相信有限元结果

具体工程情境：变截面拉杆的刚度校核

把一个左端固定、右端受拉的轻量化拉杆简化成一维轴向杆。默认 \(L=1.2\,\mathrm{m}\)、左端面积 \(A_0=400\,\mathrm{mm^2}\)，面积沿长度按 \(A(x)=A_0(1+0.5x/L)\) 增长，\(E=200\,\mathrm{GPa}\)，右端载荷 \(F=10\,\mathrm{kN}\)，用 \(n=4\) 个两节点杆单元近似。这个 toy model 足够小，可以同时看到刚度矩阵组装、边界条件、反力平衡和解析误差，而不被软件云图遮住。

必须提交的预测

  1. 把 \(n\) 加倍后，变截面杆的节点位移误差会下降、上升，还是完全不变？
  2. 左端固定反力与右端 \(F\) 的方向关系是什么？网格加密会不会改变整体平衡？
  3. 杆单元刚度矩阵来自什么局部结构，组装时如何把相邻单元共享节点的贡献相加？

只有完成预测后才展开矩阵、SVG 和收敛表；展开后可以改网格与 taper，观察预测在哪些参数下仍成立。

正式公式桥：从形函数到可核对的解析解

线性两节点杆单元在长度 \(\ell_e\) 内取

\[ u(x)=N_1u_i+N_2u_j,\qquad B=\begin{bmatrix}-1/\ell_e&1/\ell_e\end{bmatrix}. \]

若单元内取平均面积 \(A_e\)，则

\[ k_e=\int_0^{\ell_e}B^T E A_e B\,dx =\frac{EA_e}{\ell_e} \begin{bmatrix}1&-1\\-1&1\end{bmatrix}. \]

把每个 \(k_e\) 按全局节点号累加得到 \(K\)，施加 \(u(0)=0\) 后解

\[ K_{ff}u_f=F_f, \]

再用 \(R=Ku-F\) 读出约束反力。对本页的线性变截面杆，轴力处处为 \(F\)，所以连续模型的解析位移是

\[ u_{\mathrm{exact}}(x)=\frac{F}{E}\int_0^x\frac{d\xi}{A(\xi)} =\frac{FL}{EA_0\,\mathrm{taper}} \ln\!\left(1+\mathrm{taper}\frac{x}{L}\right), \]

当 \(\mathrm{taper}=0\) 时取极限 \(u=Fx/(EA_0)\)。因此这里的误差是可直接计算的 \(u_h-u_{\mathrm{exact}}\)，而不是“看云图觉得差不多”。对恒截面杆，线性杆单元在节点上会恰好精确；变截面设置则让网格收敛有可见内容。

误区与模型边界

  • 反力平衡是必要检查，不是网格收敛证明；一个边界条件错误的模型也可能在错误问题上精确平衡。
  • 位移由方程直接解出，单元应力还要对位移求导，精度与收敛行为不同；不能用位移漂亮就宣称局部应力可信。
  • 本 lab 只有一维轴向、小变形、线弹性、理想固定端和末端集中载荷；弯曲、接触、屈曲、塑性、三维应力、点载荷奇异和真实装配柔度都被排除在模型边界外。
  • 解析误差只对应明确写出的 \(A(x)\)、载荷和边界；换成复杂几何或非线性材料后，仍需独立的手算量级、反力检查与网格无关性证据。

对应 lab 的静态 fallback

无 JavaScript 时仍可核对：默认 \(L=1.2\,\mathrm{m}\)、\(A_0=4.0\times10^{-4}\,\mathrm{m^2}\)、\(\mathrm{taper}=0.5\)、\(E=200\,\mathrm{GPa}\)、\(F=10000\,\mathrm{N}\)、选中网格 \(n=4\)。

\(n\) 最大节点误差 (m) 相对 tip 误差 \(u_h(L)\) (mm) \(R(0)\) (N)
1 \(1.6395\times10^{-6}\) \(1.3479\%\) 0.120000 -10000
2 \(4.2741\times10^{-7}\) \(0.3514\%\) 0.121212 -10000
4 \(1.0808\times10^{-7}\) \(0.0889\%\) 0.121531 -10000
8 \(2.7100\times10^{-8}\) \(0.0223\%\) 0.121612 -10000
16 \(6.7800\times10^{-9}\) \(0.0056\%\) 0.121633 -10000
32 \(1.6953\times10^{-9}\) \(0.0014\%\) 0.121638 -10000

连续解析端位移为 \(u_{\mathrm{exact}}(L)=0.121640\,\mathrm{mm}\)（四舍五入），所以网格加密的方向是向它靠近；所有网格的 \(R(0)+F\) 都应在舍入误差内为零。脚本可用时还会显示 \(K\) 左上角组装预览和 FE/解析位移形状。

一、原理：从能量法到刚度方程

思路（承第五页最小势能原理）：

  1. 把连续体离散成有限个单元，单元由节点连接；
  2. 在每个单元内，用形函数把位移场表示为节点位移的插值：\(\mathbf{u} = \mathbf{N}\mathbf{u}_e\)；
  3. 由几何关系得应变 \(\boldsymbol{\varepsilon}=\mathbf{B}\mathbf{u}_e\)，由本构得应力；
  4. 写出总势能，令其对节点位移取极小 → 单元刚度矩阵：
\[\mathbf{k}_e = \int_V \mathbf{B}^T\mathbf{D}\mathbf{B}\,dV\]
  1. 组装成整体方程并施加边界条件：
\[\mathbf{K}\mathbf{u} = \mathbf{F}\]

核心认识：有限元是能量法的数值化（第五页说过）。它不是"新理论"，而是把已有的力学原理用分片近似 + 数值积分实现。

求解后：位移 → 应变 → 应力。注意精度递降——位移是直接解出的，应力是由位移微分得到的，精度天然低一个量级。这是"应力结果要打折看"的数学理由。

二、单元类型与常见陷阱

单元 用途 注意
杆/梁 桁架、框架 高效，但不给局部应力细节
壳 薄壁结构（钣金、容器） 厚度方向假设，厚壁不适用
实体（四面体/六面体） 一般三维 一阶四面体过刚，慎用
接触单元 装配、过盈、摩擦 非线性、收敛难

几个必须知道的数值陷阱：

最后一条极其重要：看到"网格加密后应力还在涨"，先怀疑是奇异点而非真实应力。 处理办法是加圆角（承第一页：真实零件本来就有圆角）、用分布载荷代替点载荷、或改用断裂力学方法。

三、用对有限元：一份检查清单

这是本页最实用的部分。有限元不出错的关键不在软件，而在建模判断：

① 边界条件是最大的误差来源。 真实结构的支撑既不是完全固定也不是完全自由。过约束会使结构显得过硬、应力被低估或转移。建模时先问：这个零件在真实装配中到底被什么约束着？

② 网格无关性验证：至少算两套网格，看关注量是否收敛。不做这一步的结果不可信。

③ 用简化模型交叉验证：先用手算或材料力学公式估一个量级，与有限元对比。数量级不符时，几乎总是建模错了，而不是手算错了。

④ 检查平衡：反力总和是否等于外载？

⑤ 变形形状比应力云图更可靠：先看变形动画是否符合物理直觉——变形方向不对，说明边界条件或载荷有误。

⑥ 关注量决定网格策略：算刚度可以用粗网格，算局部应力必须在关注区细化。

⑦ 材料与非线性：是否需要塑性、大变形、接触？线性分析在接触与失稳问题上会给出完全错误的结果。

一句总结：有限元不会告诉你模型错了，它只会认真地解一个错误的问题。 这与 🔗 微电子站"抽象会漏"、自动化站"模型永远是错的"是同一种警觉。

四、常见分析类型

五、结构优化：让计算机来设计

三个层次（自由度递增）：

① 尺寸优化：厚度、截面参数； ② 形状优化：边界形状（如优化圆角轮廓以降低应力集中）； ③ 拓扑优化：连材料该不该存在都由算法决定。

拓扑优化的基本形式（SIMP 方法）：

\[\min_{\rho} \ \mathbf{u}^T\mathbf{K}(\rho)\mathbf{u} \quad \text{s.t.}\ \sum \rho_i v_i \le V^*,\ \ \rho_i\in[0,1]\]

把每个单元的"存在与否"松弛为连续密度 \(\rho\)，用惩罚指数抑制中间密度，用伴随法高效求梯度（🔗 grad-math 优化线；伴随法与深度学习的反向传播是同一思想）。

结果往往像骨骼或树枝——因为它们解决的是同一个问题：在给定材料下最有效地传递载荷。（🔗 生物站：骨小梁沿主应力方向排列，是演化"算"出的拓扑优化解。）

工程现实：

六、仿真在产品开发中的位置

与试验的关系：仿真不能取代试验，二者互补。

数字孪生（第二十二页）是这个思路的延伸：让模型与实物持续同步。

七、要点

下一页：本站的最后一页正文——这个行业正在往哪里走，以及本硕博的实际去向。