简介面向具备材料力学基础与一定 Python 编程能力的高校学生、科研人员及工程技术人员这份悬臂梁变形分析资料聚焦受力条件下梁体挠度、转角、弯矩与剪力的数值求解与可视化呈现。内容围绕悬臂梁力学模型的建立展开结合材料力学基本理论可采用有限差分法或解析法求解微分方程并通过 Python 编程实现仿真输出变形曲线与内力分布便于对照理论解验证算法正确性。压缩包共 1 个 docx 文件、约 28KB以文字与代码说明为主篇幅精简却完整覆盖建模思路、离散化处理与结果解读方便直接嵌入教学演示或个人学习笔记。目前已有 51 人学习下载。读者可据此掌握从力学建模到编程仿真的完整流程并尝试调整梁长、载荷大小与边界条件观察挠度与转角的变化规律为后续复杂结构分析与课程设计积累可复现的代码框架与排错经验。1. 从一根“算不准”的悬臂梁说起这套 Python 方案到底解决什么问题悬臂梁变形分析几乎是结构力学课程设计的标配题目但真正动手做的时候绝大多数人卡在同一个地方自由端集中力有现成公式均布荷载叠加上去也能推可一旦要同时给出挠曲线、转角、弯矩、剪力图或者导师临时把梁长、截面、边界条件改一个手推的那套就全废了。这套方案做的事情很具体——用 Python 搭一条可复现的数值仿真流程把悬臂梁在集中力、均布荷载、端部力偶三类典型工况下的变形算清楚把四张图一次出全并且用解析解把数值误差钉死在明面上。它适合正在做结构力学课设、力学方向毕设的同学也适合需要快速搭一个梁单元计算原型的工程技术人员。科学计算的门槛不在公式在于你知不知道自己算出来的数可不可信这套代码就是围绕“可信”两个字组织的。2. 力学模型与计算路线解析式打底、有限差分兜底的双通道设计先把基本盘说清楚。一根等截面悬臂梁的小变形问题本质是解四阶常微分方程 EI·w q(x)或者退一步先拿到弯矩分布 M(x)再解二阶方程 EI·w M(x)。两条路都能走通但适用面不一样。方案里我把解析解和有限差分都保留下来不是为了炫技而是验收的时候手里必须有一把尺子数值解算出来是多是少得有个东西对照。只留数值解等于把黑匣子交给导师答辩时一问误差来源就露怯。2.1 符号约定与边界条件先把“正方向”钉死符号约定是这类仿真里最不值钱、也最容易翻车的一环。我采用的约定是x 轴从固定端指向自由端挠度 w 向下为正弯矩 M 以使梁产生正曲率下凹为正。在这个约定下自由端受向下集中力 P 时截面弯矩写成 M(x) P(L − x)积分两次得到的挠曲线 w(x) Px²(3L − x)/(6EI) 在自由端取正值和“向下挠”的物理直觉一致。如果换成“向上为正”整套公式里的符号要跟着翻一遍代码里任何一处漏改画出来的挠曲线就是反的。边界条件同样要提前写死。固定端 x 0 处位移和转角同时为零即 w(0) 0、w(0) 0自由端 x L 处弯矩和剪力为零即 M(L) 0、V(L) 0。前者是本质边界条件必须在方程组里显式施加后者是自然边界条件在二阶方程里通过弯矩分布本身体现不需要额外加方程。这个区分很关键下一章组装矩阵时会反复用到。提示把符号约定写成代码顶部的一行注释任何人改代码之前先看这一行。我见过太多方案公式推导和代码实现各用一套正方向最后挠度图上下颠倒排查花掉一晚上。2.2 解析解能覆盖的工况和它覆盖不了的部分三类标准工况的解析式先列成表后面既当基准也当教学对照。工况弯矩 M(x)挠度 w(x)自由端挠度自由端集中力 PP(L − x)Px²(3L − x)/(6EI)PL³/(3EI)全跨均布荷载 qq(L − x)²/2qx²(6L² − 4Lx x²)/(24EI)qL⁴/(8EI)自由端集中力偶 M₀M₀M₀x²/(2EI)M₀L²/(2EI)解析解的好处是干净、可手算、适合写进报告做理论推导。但它的边界很清楚荷载分段、截面沿轴线变化、需要同时输出内力图、或者后续要扩展到超静定结构时逐段推公式的工作量会迅速失控。课程设计里常见的一种误用是拿单工况公式去硬套多工况叠加叠加本身没错错在忘了检查每段荷载的适用区间。数值方法接手的正是这些解析式不方便展开的地方而解析式反过来负责给数值解做标定两者是配合关系不是替代关系。2.3 有限差分离散三对角方程组的组装逻辑数值路线选的是二阶中心差分理由很实际二阶方程只需要两个边界条件未知量少矩阵规模小边界处理比四阶差分干净。四阶差分要引入 w₋₁、w₋₂ 和 w₍N1₎、w₍N2₎ 四个虚节点固定端和自由端各要处理两层代码量和出错概率都翻倍。把梁均分为 N 段步长 h L/N节点 xᵢ i·hi 0…N。内部节点 i 1…N−1 上二阶导数用中心差分近似EI·(wᵢ₊₁ − 2wᵢ wᵢ₋₁)/h² Mᵢ固定端节点 0 处w₀ 0 已知转角为零这个条件用虚节点处理令 w₋₁ w₁代入节点 0 的差分方程得到 2w₁ − 2w₀ h²M₀/EI由于 w₀ 0化简为 2w₁ h²M₀/EI。这样固定端不是简单粗暴地令 w₁ 0而是保留了二阶精度收敛阶不会掉到一阶。未知量取 w₁…w_N 共 N 个方程取节点 0 到 N−1 共 N 个正好闭合。自由端节点 N 不作为独立方程出现它的值由节点 N−1 的方程给出弯矩为零这个条件在离散意义下近似满足。组装出来的系数矩阵是一个下三角加一条次对角线的三对角结构用 numpy 的线性求解器直接解规模到几千个节点也毫秒级完成。3. Python 实现参数、求解、后处理、批量扫描四个模块代码组织上不追求大而全按职责切成四块参数与荷载、求解器、后处理、网格扫描。每块单独可测改工况的时候只动第一块排查问题的时候能快速定位到具体环节。3.1 参数与荷载模块把工况做成一个可配置的数据结构from dataclasses import dataclass dataclass class BeamCase: 悬臂梁单工况参数长度单位统一为 m力统一为 N L: float 2.0 # 梁长 E: float 210e9 # 弹性模量钢约 210 GPa width: float 0.05 # 矩形截面宽 height: float 0.01 # 矩形截面高 P: float 0.0 # 自由端集中力 q: float 0.0 # 全跨均布荷载 M0: float 0.0 # 自由端集中力偶 N: int 200 # 分段数 property def I(self) - float: 矩形截面惯性矩单位 m^4 return self.width * self.height ** 3 / 12.0 property def EI(self) - float: 抗弯刚度单位 N·m^2 return self.E * self.I用 dataclass 而不是散落的全局变量好处是每个工况就是一次构造批量扫描时直接改字段、重新实例化不会互相污染。I 和 EI 写成 property避免手动计算惯性矩时把 height 的三次方写成平方——这类笔误在截面尺寸跨数量级时会被放大成几个数量级的偏差靠肉眼很难发现。字段命名上刻意避开 h 这个字母因为它在数值部分通常指步长混用会让后面读代码的人分不清是截面高还是网格间距。3.2 求解器解析解与有限差分走同一个接口两个求解函数返回完全相同的字典结构后处理只认这个结构不关心数据是推出来的还是算出来的。import numpy as np def solve_analytical(case: BeamCase) - dict: 三类荷载线性叠加的解析解 x np.linspace(0, case.L, case.N 1) EI, L case.EI, case.L w np.zeros_like(x) if case.P: w case.P * x**2 * (3 * L - x) / (6 * EI) if case.q: w case.q * x**2 * (6 * L**2 - 4 * L * x x**2) / (24 * EI) if case.M0: w case.M0 * x**2 / (2 * EI) theta np.gradient(w, x) # 转角一阶中心差分 M EI * np.gradient(theta, x) # 弯矩 V np.gradient(M, x) # 剪力 return {x: x, w: w, theta: theta, M: M, V: V} def solve_fd(case: BeamCase) - dict: 有限差分求解未知量 w_1..w_N N, L, EI case.N, case.L, case.EI x np.linspace(0, L, N 1) dx L / N # 节点弯矩的解析表达式集中力 均布 端部力偶 M (case.P * (L - x) case.q * (L - x) ** 2 / 2.0 case.M0 * np.ones_like(x)) A np.zeros((N, N)) # 系数矩阵N 个未知量对应 N 个方程 rhs np.zeros(N) # 右端项 # 节点 0固定端虚节点 w_{-1} w_1w_0 0 A[0, 0] 2.0 rhs[0] dx**2 * M[0] / EI # 节点 i 1..N-1w_{i-1} - 2w_i w_{i1} h^2 M_i / EI for i in range(1, N): if i 2: A[i, i - 2] 1.0 A[i, i - 1] -2.0 A[i, i] 1.0 rhs[i] dx**2 * M[i] / EI w_inner np.linalg.solve(A, rhs) w np.concatenate(([0.0], w_inner)) theta np.gradient(w, x) M_num EI * np.gradient(theta, x) V_num np.gradient(M_num, x) return {x: x, w: w, theta: theta, M: M_num, V: V_num}solve_fd里弯矩 M 直接采用解析式而不是由数值解反推这是有意的弯矩分布是荷载的已知结果没必要再引入一层数值误差。矩阵 A 的第 0 行是[2, 0, …]从第 1 行开始是[-2, 1, 0…]、[1, -2, 1, 0…]这样的三对角格式np.linalg.solve对这种结构非常稳定。参数 N 控制离散精度200 段对一根 2 米左右的梁已经能压到 10⁻⁶ 量级的相对误差日常课设够用做收敛性研究时可以拉到 320 或 640。两个函数接口一致后处理里可以并排调用也可以只调其一。这也是排查问题的便利之处数值解和解析解对不上时先看解析解自身是否写错再看数值边界处理路径是清晰的。3.3 后处理一张画布四张子图挠度转角弯矩剪力全出import matplotlib.pyplot as plt def plot_result(res: dict, title: str ) - None: x, w, theta, M, V res[x], res[w], res[theta], res[M], res[V] fig, axes plt.subplots(2, 2, figsize(11, 7)) fig.suptitle(title) axes[0, 0].plot(x, w * 1e3, b-) axes[0, 0].set_ylabel(挠度 w / mm) axes[0, 0].invert_yaxis() # 向下为正翻转纵轴更符合直觉 axes[0, 1].plot(x, np.degrees(theta), r-) axes[0, 1].set_ylabel(转角 θ / °) axes[1, 0].plot(x, M, g-) axes[1, 0].set_ylabel(弯矩 M / (N·m)) axes[1, 1].plot(x, V, k-) axes[1, 1].set_ylabel(剪力 V / N) for ax in axes.flat: ax.set_xlabel(x / m) ax.grid(alpha0.3) plt.tight_layout() plt.show()四张子图的纵轴单位分别换成毫米、度、牛米、牛报告里直接截图就能用。invert_yaxis()那一行是为了对齐工程直觉w 向下为正图上把纵轴翻转后梁向下弯的形态和肉眼看到的一致。转角用np.degrees转成度比弧度更适合作图展示。弯矩和剪力图直接画数值解的结果和解析式放在同一张图上对比时两条线应该几乎重合肉眼可见的偏离通常意味着荷载或边界条件写错了。3.4 网格收敛性扫描用误差曲线判断解有没有算对def grid_study(case: BeamCase, Ns(10, 20, 40, 80, 160, 320)) - list: 以解析解自由端挠度为基准返回 (N, w_end, 相对误差) 列表 ref solve_analytical(case) w_ref ref[w][-1] rows [] for N in Ns: case.N N res solve_fd(case) rel_err abs(res[w][-1] - w_ref) / abs(w_ref) rows.append((N, res[w][-1], rel_err)) return rows这个函数是整个方案里最有说服力的一段。二阶中心差分的理论收敛阶是 O(h²)网格加倍时相对误差应该降到大约四分之一。把grid_study的输出画成 log-log 曲线斜率接近 −2 就说明离散格式和边界处理都是对的如果斜率只有 −1 左右基本可以断定固定端用了 w₁ 0 这种一阶近似或者自由端多加了不该加的约束。答辩时这张图比任何口头解释都管用。4. 避坑与排查悬臂梁数值仿真里最容易翻车的五个地方这一章按“现象 → 原因 → 解决”组织每条都是实际项目里反复出现的。数值仿真有一半时间花在排查上提前知道坑在哪比事后 debug 省事得多。4.1 符号约定不统一挠度图上下颠倒现象集中力明明向下画出来的挠曲线却朝上或者挠度图和弯矩图的凹凸方向对不上看起来像两条互不相关的曲线。原因w 的正方向和 M 的正方向各自定义没有形成闭环。比如推导时用“向上为正”代码里画图时又按“向下为正”翻转中间积分常数的符号就翻了。更隐蔽的一种情况是弯矩公式写成 M −P(L−x)但挠曲线仍用正号积分结果数值全对、符号全反。解决在代码顶部用一行注释钉死约定所有公式围绕这一行推导。解析解在自由端取值的符号要和物理直觉一致作为第一道校验。画图前统一处理正方向要么翻转纵轴要么在轴标签上明确标注“向下为正”。我一般会在plot_result里加一句断言检查w[-1]的符号是否和荷载方向一致符号错了直接报错不进入画图环节。4.2 固定端用 w₁ 0 近似收敛阶从二阶掉到一阶现象网格从 40 加密到 160误差只降了一半而不是降到十六分之一log-log 收敛曲线的斜率在 −1 附近徘徊。原因固定端转角为零这个条件被简化成了“第二个节点位移为零”。这相当于用一阶精度的差分代替中心差分局部截断误差是 O(h)全局误差跟着降为一阶。很多教材示例为了讲解方便会这么写照搬到代码里就会拖累整个求解器。解决引入虚节点令 w₋₁ w₁代入节点 0 的中心差分方程得到 2w₁ h²M₀/EI。这样固定端保留了二阶精度收敛曲线斜率回到 −2 附近。改动只有一行矩阵元素和一个右端项代价很小收益是整条误差曲线好看一个量级。4.3 自由端缺方程矩阵奇异或末端出现锯齿振荡现象np.linalg.solve抛出LinAlgError: Singular matrix或者勉强解出来但自由端附近挠度剧烈振荡和解析解完全对不上。原因未知量个数和方程个数没对齐。一种情况是把 w₀ 也算进未知量又只写了 N−1 个内部方程另一种是误以为自由端也要加位移约束硬塞一个 w_N 0 进去把悬臂梁做成了简支梁。还有一种是弯矩数组长度和节点数差一位矩阵维度对不上。解决明确未知量是 w₁…w_N 共 N 个方程是节点 0 到 N−1 共 N 个自由端不附加任何位移条件弯矩为零通过 M(L) 0 在荷载分布里体现。组装完矩阵后打印一下A.shape和rhs.shape两者不一致就立刻停下来查。自由端振荡还有一种可能是弯矩数组在末端用了错误的索引检查M[-1]是否确实接近零。4.4 单位制混用EI 量级错出三个数量级现象算出来的自由端挠度是几米甚至几十米明显超出小变形的合理范围但代码逻辑检查不出毛病。原因弹性模量用 GPa 填写却没乘 10⁹截面尺寸用毫米填写却和米制的梁长混算惯性矩里 height 的三次方算成了平方。钢的弹性模量是 210 GPa写成 210 直接代入EI 会小九个数量级挠度相应放大九个数量级。截面从毫米换到米漏掉 10⁻³惯性矩就差 10⁻¹²同样致命。解决全部统一到国际单位制弹性模量写210e9而不是210截面尺寸用米。在BeamCase里加一个量纲自检方法打印 EI 的数量级矩形钢截面一般在 10⁴ 到 10⁶ N·m² 之间超出这个范围基本可以断定单位写错。挠度结果也做一次合理性检查小变形假设下自由端挠度不应超过梁长的十分之一超了先查单位再查公式。4.5 荷载突变点直接差分内力图出现尖刺现象挠曲线看着光滑正常但用np.gradient求弯矩或剪力时集中力作用位置附近冒出一个尖峰幅值明显不合理。原因集中力处剪力发生跳跃弯矩的一阶导数不连续中心差分跨过间断点时会把两侧的值平均产生虚假的尖刺。网格越粗尖刺越明显网格加密后尖刺变窄但不会消失。解决把集中力作用点设为网格节点保证间断点落在节点上而不是节点之间内力图优先用解析弯矩表达式绘制数值微分只作为辅助验证确实要对数值解求内力时在突变点两侧分段差分不要跨点。均布荷载的端点也有类似问题处理思路相同。5. 验收与答辩用三张表和一条收敛曲线把模型说清楚验收环节不需要讲太多理论把证据摆出来就行。我通常准备三张表加一张图。第一张是输入参数表梁长、截面宽高、弹性模量、惯性矩、三类荷载的取值全部标单位让评审一眼能复算。第二张是对比表三种工况下解析解与数值解的自由端挠度并排列出最后一列写相对误差正常应该落在 10⁻⁴ 到 10⁻⁶ 量级。工况解析解 w(L) / mm数值解 w(L) / mm相对误差集中力 P 500 N由 PL³/(3EI) 给出求解器输出目标 10⁻⁴均布 q 2 kN/m由 qL⁴/(8EI) 给出求解器输出目标 10⁻⁴端部力偶 M₀ 300 N·m由 M₀L²/(2EI) 给出求解器输出目标 10⁻⁴第三张是网格收敛表N 取 10、20、40、80、160、320每行记录自由端挠度和相对误差相邻两行的误差比值应该接近 4。把这张表画成双对数曲线斜率接近 −2离散格式的正确性就有了直接证据。答辩时被问到“你怎么知道算对了”指这张图比解释半天公式推导有效得多。还有三个问题几乎每次都会被问提前准备好答案为什么用二阶差分而不是四阶——未知量少、边界处理干净、二阶精度对小变形梁足够固定端怎么处理的——虚节点加中心差分收敛阶保持二阶误差来源有哪些——离散截断误差、边界条件近似、数值微分求内力的额外误差三条分开说别混成一团。从那以后我每次交方案前都强制跑一遍网格收敛哪怕解析解已经本文还有配套的精品资源点击获取 Python有限差分悬臂梁变形仿真与网格收敛验证
简介面向具备材料力学基础与一定 Python 编程能力的高校学生、科研人员及工程技术人员这份悬臂梁变形分析资料聚焦受力条件下梁体挠度、转角、弯矩与剪力的数值求解与可视化呈现。内容围绕悬臂梁力学模型的建立展开结合材料力学基本理论可采用有限差分法或解析法求解微分方程并通过 Python 编程实现仿真输出变形曲线与内力分布便于对照理论解验证算法正确性。压缩包共 1 个 docx 文件、约 28KB以文字与代码说明为主篇幅精简却完整覆盖建模思路、离散化处理与结果解读方便直接嵌入教学演示或个人学习笔记。目前已有 51 人学习下载。读者可据此掌握从力学建模到编程仿真的完整流程并尝试调整梁长、载荷大小与边界条件观察挠度与转角的变化规律为后续复杂结构分析与课程设计积累可复现的代码框架与排错经验。1. 从一根“算不准”的悬臂梁说起这套 Python 方案到底解决什么问题悬臂梁变形分析几乎是结构力学课程设计的标配题目但真正动手做的时候绝大多数人卡在同一个地方自由端集中力有现成公式均布荷载叠加上去也能推可一旦要同时给出挠曲线、转角、弯矩、剪力图或者导师临时把梁长、截面、边界条件改一个手推的那套就全废了。这套方案做的事情很具体——用 Python 搭一条可复现的数值仿真流程把悬臂梁在集中力、均布荷载、端部力偶三类典型工况下的变形算清楚把四张图一次出全并且用解析解把数值误差钉死在明面上。它适合正在做结构力学课设、力学方向毕设的同学也适合需要快速搭一个梁单元计算原型的工程技术人员。科学计算的门槛不在公式在于你知不知道自己算出来的数可不可信这套代码就是围绕“可信”两个字组织的。2. 力学模型与计算路线解析式打底、有限差分兜底的双通道设计先把基本盘说清楚。一根等截面悬臂梁的小变形问题本质是解四阶常微分方程 EI·w q(x)或者退一步先拿到弯矩分布 M(x)再解二阶方程 EI·w M(x)。两条路都能走通但适用面不一样。方案里我把解析解和有限差分都保留下来不是为了炫技而是验收的时候手里必须有一把尺子数值解算出来是多是少得有个东西对照。只留数值解等于把黑匣子交给导师答辩时一问误差来源就露怯。2.1 符号约定与边界条件先把“正方向”钉死符号约定是这类仿真里最不值钱、也最容易翻车的一环。我采用的约定是x 轴从固定端指向自由端挠度 w 向下为正弯矩 M 以使梁产生正曲率下凹为正。在这个约定下自由端受向下集中力 P 时截面弯矩写成 M(x) P(L − x)积分两次得到的挠曲线 w(x) Px²(3L − x)/(6EI) 在自由端取正值和“向下挠”的物理直觉一致。如果换成“向上为正”整套公式里的符号要跟着翻一遍代码里任何一处漏改画出来的挠曲线就是反的。边界条件同样要提前写死。固定端 x 0 处位移和转角同时为零即 w(0) 0、w(0) 0自由端 x L 处弯矩和剪力为零即 M(L) 0、V(L) 0。前者是本质边界条件必须在方程组里显式施加后者是自然边界条件在二阶方程里通过弯矩分布本身体现不需要额外加方程。这个区分很关键下一章组装矩阵时会反复用到。提示把符号约定写成代码顶部的一行注释任何人改代码之前先看这一行。我见过太多方案公式推导和代码实现各用一套正方向最后挠度图上下颠倒排查花掉一晚上。2.2 解析解能覆盖的工况和它覆盖不了的部分三类标准工况的解析式先列成表后面既当基准也当教学对照。工况弯矩 M(x)挠度 w(x)自由端挠度自由端集中力 PP(L − x)Px²(3L − x)/(6EI)PL³/(3EI)全跨均布荷载 qq(L − x)²/2qx²(6L² − 4Lx x²)/(24EI)qL⁴/(8EI)自由端集中力偶 M₀M₀M₀x²/(2EI)M₀L²/(2EI)解析解的好处是干净、可手算、适合写进报告做理论推导。但它的边界很清楚荷载分段、截面沿轴线变化、需要同时输出内力图、或者后续要扩展到超静定结构时逐段推公式的工作量会迅速失控。课程设计里常见的一种误用是拿单工况公式去硬套多工况叠加叠加本身没错错在忘了检查每段荷载的适用区间。数值方法接手的正是这些解析式不方便展开的地方而解析式反过来负责给数值解做标定两者是配合关系不是替代关系。2.3 有限差分离散三对角方程组的组装逻辑数值路线选的是二阶中心差分理由很实际二阶方程只需要两个边界条件未知量少矩阵规模小边界处理比四阶差分干净。四阶差分要引入 w₋₁、w₋₂ 和 w₍N1₎、w₍N2₎ 四个虚节点固定端和自由端各要处理两层代码量和出错概率都翻倍。把梁均分为 N 段步长 h L/N节点 xᵢ i·hi 0…N。内部节点 i 1…N−1 上二阶导数用中心差分近似EI·(wᵢ₊₁ − 2wᵢ wᵢ₋₁)/h² Mᵢ固定端节点 0 处w₀ 0 已知转角为零这个条件用虚节点处理令 w₋₁ w₁代入节点 0 的差分方程得到 2w₁ − 2w₀ h²M₀/EI由于 w₀ 0化简为 2w₁ h²M₀/EI。这样固定端不是简单粗暴地令 w₁ 0而是保留了二阶精度收敛阶不会掉到一阶。未知量取 w₁…w_N 共 N 个方程取节点 0 到 N−1 共 N 个正好闭合。自由端节点 N 不作为独立方程出现它的值由节点 N−1 的方程给出弯矩为零这个条件在离散意义下近似满足。组装出来的系数矩阵是一个下三角加一条次对角线的三对角结构用 numpy 的线性求解器直接解规模到几千个节点也毫秒级完成。3. Python 实现参数、求解、后处理、批量扫描四个模块代码组织上不追求大而全按职责切成四块参数与荷载、求解器、后处理、网格扫描。每块单独可测改工况的时候只动第一块排查问题的时候能快速定位到具体环节。3.1 参数与荷载模块把工况做成一个可配置的数据结构from dataclasses import dataclass dataclass class BeamCase: 悬臂梁单工况参数长度单位统一为 m力统一为 N L: float 2.0 # 梁长 E: float 210e9 # 弹性模量钢约 210 GPa width: float 0.05 # 矩形截面宽 height: float 0.01 # 矩形截面高 P: float 0.0 # 自由端集中力 q: float 0.0 # 全跨均布荷载 M0: float 0.0 # 自由端集中力偶 N: int 200 # 分段数 property def I(self) - float: 矩形截面惯性矩单位 m^4 return self.width * self.height ** 3 / 12.0 property def EI(self) - float: 抗弯刚度单位 N·m^2 return self.E * self.I用 dataclass 而不是散落的全局变量好处是每个工况就是一次构造批量扫描时直接改字段、重新实例化不会互相污染。I 和 EI 写成 property避免手动计算惯性矩时把 height 的三次方写成平方——这类笔误在截面尺寸跨数量级时会被放大成几个数量级的偏差靠肉眼很难发现。字段命名上刻意避开 h 这个字母因为它在数值部分通常指步长混用会让后面读代码的人分不清是截面高还是网格间距。3.2 求解器解析解与有限差分走同一个接口两个求解函数返回完全相同的字典结构后处理只认这个结构不关心数据是推出来的还是算出来的。import numpy as np def solve_analytical(case: BeamCase) - dict: 三类荷载线性叠加的解析解 x np.linspace(0, case.L, case.N 1) EI, L case.EI, case.L w np.zeros_like(x) if case.P: w case.P * x**2 * (3 * L - x) / (6 * EI) if case.q: w case.q * x**2 * (6 * L**2 - 4 * L * x x**2) / (24 * EI) if case.M0: w case.M0 * x**2 / (2 * EI) theta np.gradient(w, x) # 转角一阶中心差分 M EI * np.gradient(theta, x) # 弯矩 V np.gradient(M, x) # 剪力 return {x: x, w: w, theta: theta, M: M, V: V} def solve_fd(case: BeamCase) - dict: 有限差分求解未知量 w_1..w_N N, L, EI case.N, case.L, case.EI x np.linspace(0, L, N 1) dx L / N # 节点弯矩的解析表达式集中力 均布 端部力偶 M (case.P * (L - x) case.q * (L - x) ** 2 / 2.0 case.M0 * np.ones_like(x)) A np.zeros((N, N)) # 系数矩阵N 个未知量对应 N 个方程 rhs np.zeros(N) # 右端项 # 节点 0固定端虚节点 w_{-1} w_1w_0 0 A[0, 0] 2.0 rhs[0] dx**2 * M[0] / EI # 节点 i 1..N-1w_{i-1} - 2w_i w_{i1} h^2 M_i / EI for i in range(1, N): if i 2: A[i, i - 2] 1.0 A[i, i - 1] -2.0 A[i, i] 1.0 rhs[i] dx**2 * M[i] / EI w_inner np.linalg.solve(A, rhs) w np.concatenate(([0.0], w_inner)) theta np.gradient(w, x) M_num EI * np.gradient(theta, x) V_num np.gradient(M_num, x) return {x: x, w: w, theta: theta, M: M_num, V: V_num}solve_fd里弯矩 M 直接采用解析式而不是由数值解反推这是有意的弯矩分布是荷载的已知结果没必要再引入一层数值误差。矩阵 A 的第 0 行是[2, 0, …]从第 1 行开始是[-2, 1, 0…]、[1, -2, 1, 0…]这样的三对角格式np.linalg.solve对这种结构非常稳定。参数 N 控制离散精度200 段对一根 2 米左右的梁已经能压到 10⁻⁶ 量级的相对误差日常课设够用做收敛性研究时可以拉到 320 或 640。两个函数接口一致后处理里可以并排调用也可以只调其一。这也是排查问题的便利之处数值解和解析解对不上时先看解析解自身是否写错再看数值边界处理路径是清晰的。3.3 后处理一张画布四张子图挠度转角弯矩剪力全出import matplotlib.pyplot as plt def plot_result(res: dict, title: str ) - None: x, w, theta, M, V res[x], res[w], res[theta], res[M], res[V] fig, axes plt.subplots(2, 2, figsize(11, 7)) fig.suptitle(title) axes[0, 0].plot(x, w * 1e3, b-) axes[0, 0].set_ylabel(挠度 w / mm) axes[0, 0].invert_yaxis() # 向下为正翻转纵轴更符合直觉 axes[0, 1].plot(x, np.degrees(theta), r-) axes[0, 1].set_ylabel(转角 θ / °) axes[1, 0].plot(x, M, g-) axes[1, 0].set_ylabel(弯矩 M / (N·m)) axes[1, 1].plot(x, V, k-) axes[1, 1].set_ylabel(剪力 V / N) for ax in axes.flat: ax.set_xlabel(x / m) ax.grid(alpha0.3) plt.tight_layout() plt.show()四张子图的纵轴单位分别换成毫米、度、牛米、牛报告里直接截图就能用。invert_yaxis()那一行是为了对齐工程直觉w 向下为正图上把纵轴翻转后梁向下弯的形态和肉眼看到的一致。转角用np.degrees转成度比弧度更适合作图展示。弯矩和剪力图直接画数值解的结果和解析式放在同一张图上对比时两条线应该几乎重合肉眼可见的偏离通常意味着荷载或边界条件写错了。3.4 网格收敛性扫描用误差曲线判断解有没有算对def grid_study(case: BeamCase, Ns(10, 20, 40, 80, 160, 320)) - list: 以解析解自由端挠度为基准返回 (N, w_end, 相对误差) 列表 ref solve_analytical(case) w_ref ref[w][-1] rows [] for N in Ns: case.N N res solve_fd(case) rel_err abs(res[w][-1] - w_ref) / abs(w_ref) rows.append((N, res[w][-1], rel_err)) return rows这个函数是整个方案里最有说服力的一段。二阶中心差分的理论收敛阶是 O(h²)网格加倍时相对误差应该降到大约四分之一。把grid_study的输出画成 log-log 曲线斜率接近 −2 就说明离散格式和边界处理都是对的如果斜率只有 −1 左右基本可以断定固定端用了 w₁ 0 这种一阶近似或者自由端多加了不该加的约束。答辩时这张图比任何口头解释都管用。4. 避坑与排查悬臂梁数值仿真里最容易翻车的五个地方这一章按“现象 → 原因 → 解决”组织每条都是实际项目里反复出现的。数值仿真有一半时间花在排查上提前知道坑在哪比事后 debug 省事得多。4.1 符号约定不统一挠度图上下颠倒现象集中力明明向下画出来的挠曲线却朝上或者挠度图和弯矩图的凹凸方向对不上看起来像两条互不相关的曲线。原因w 的正方向和 M 的正方向各自定义没有形成闭环。比如推导时用“向上为正”代码里画图时又按“向下为正”翻转中间积分常数的符号就翻了。更隐蔽的一种情况是弯矩公式写成 M −P(L−x)但挠曲线仍用正号积分结果数值全对、符号全反。解决在代码顶部用一行注释钉死约定所有公式围绕这一行推导。解析解在自由端取值的符号要和物理直觉一致作为第一道校验。画图前统一处理正方向要么翻转纵轴要么在轴标签上明确标注“向下为正”。我一般会在plot_result里加一句断言检查w[-1]的符号是否和荷载方向一致符号错了直接报错不进入画图环节。4.2 固定端用 w₁ 0 近似收敛阶从二阶掉到一阶现象网格从 40 加密到 160误差只降了一半而不是降到十六分之一log-log 收敛曲线的斜率在 −1 附近徘徊。原因固定端转角为零这个条件被简化成了“第二个节点位移为零”。这相当于用一阶精度的差分代替中心差分局部截断误差是 O(h)全局误差跟着降为一阶。很多教材示例为了讲解方便会这么写照搬到代码里就会拖累整个求解器。解决引入虚节点令 w₋₁ w₁代入节点 0 的中心差分方程得到 2w₁ h²M₀/EI。这样固定端保留了二阶精度收敛曲线斜率回到 −2 附近。改动只有一行矩阵元素和一个右端项代价很小收益是整条误差曲线好看一个量级。4.3 自由端缺方程矩阵奇异或末端出现锯齿振荡现象np.linalg.solve抛出LinAlgError: Singular matrix或者勉强解出来但自由端附近挠度剧烈振荡和解析解完全对不上。原因未知量个数和方程个数没对齐。一种情况是把 w₀ 也算进未知量又只写了 N−1 个内部方程另一种是误以为自由端也要加位移约束硬塞一个 w_N 0 进去把悬臂梁做成了简支梁。还有一种是弯矩数组长度和节点数差一位矩阵维度对不上。解决明确未知量是 w₁…w_N 共 N 个方程是节点 0 到 N−1 共 N 个自由端不附加任何位移条件弯矩为零通过 M(L) 0 在荷载分布里体现。组装完矩阵后打印一下A.shape和rhs.shape两者不一致就立刻停下来查。自由端振荡还有一种可能是弯矩数组在末端用了错误的索引检查M[-1]是否确实接近零。4.4 单位制混用EI 量级错出三个数量级现象算出来的自由端挠度是几米甚至几十米明显超出小变形的合理范围但代码逻辑检查不出毛病。原因弹性模量用 GPa 填写却没乘 10⁹截面尺寸用毫米填写却和米制的梁长混算惯性矩里 height 的三次方算成了平方。钢的弹性模量是 210 GPa写成 210 直接代入EI 会小九个数量级挠度相应放大九个数量级。截面从毫米换到米漏掉 10⁻³惯性矩就差 10⁻¹²同样致命。解决全部统一到国际单位制弹性模量写210e9而不是210截面尺寸用米。在BeamCase里加一个量纲自检方法打印 EI 的数量级矩形钢截面一般在 10⁴ 到 10⁶ N·m² 之间超出这个范围基本可以断定单位写错。挠度结果也做一次合理性检查小变形假设下自由端挠度不应超过梁长的十分之一超了先查单位再查公式。4.5 荷载突变点直接差分内力图出现尖刺现象挠曲线看着光滑正常但用np.gradient求弯矩或剪力时集中力作用位置附近冒出一个尖峰幅值明显不合理。原因集中力处剪力发生跳跃弯矩的一阶导数不连续中心差分跨过间断点时会把两侧的值平均产生虚假的尖刺。网格越粗尖刺越明显网格加密后尖刺变窄但不会消失。解决把集中力作用点设为网格节点保证间断点落在节点上而不是节点之间内力图优先用解析弯矩表达式绘制数值微分只作为辅助验证确实要对数值解求内力时在突变点两侧分段差分不要跨点。均布荷载的端点也有类似问题处理思路相同。5. 验收与答辩用三张表和一条收敛曲线把模型说清楚验收环节不需要讲太多理论把证据摆出来就行。我通常准备三张表加一张图。第一张是输入参数表梁长、截面宽高、弹性模量、惯性矩、三类荷载的取值全部标单位让评审一眼能复算。第二张是对比表三种工况下解析解与数值解的自由端挠度并排列出最后一列写相对误差正常应该落在 10⁻⁴ 到 10⁻⁶ 量级。工况解析解 w(L) / mm数值解 w(L) / mm相对误差集中力 P 500 N由 PL³/(3EI) 给出求解器输出目标 10⁻⁴均布 q 2 kN/m由 qL⁴/(8EI) 给出求解器输出目标 10⁻⁴端部力偶 M₀ 300 N·m由 M₀L²/(2EI) 给出求解器输出目标 10⁻⁴第三张是网格收敛表N 取 10、20、40、80、160、320每行记录自由端挠度和相对误差相邻两行的误差比值应该接近 4。把这张表画成双对数曲线斜率接近 −2离散格式的正确性就有了直接证据。答辩时被问到“你怎么知道算对了”指这张图比解释半天公式推导有效得多。还有三个问题几乎每次都会被问提前准备好答案为什么用二阶差分而不是四阶——未知量少、边界处理干净、二阶精度对小变形梁足够固定端怎么处理的——虚节点加中心差分收敛阶保持二阶误差来源有哪些——离散截断误差、边界条件近似、数值微分求内力的额外误差三条分开说别混成一团。从那以后我每次交方案前都强制跑一遍网格收敛哪怕解析解已经本文还有配套的精品资源点击获取