Python实现有限元作业代码：CST单元单刚组装与求解全流程-北京尧图网络

简介浙江大学2020至2021学年春夏学期有限元方法课程作业代码包面向需要将有限元理论转化为可运行程序的学生与研究人员内容覆盖网格生成、弱形式建立、插值函数选择、矩阵组装、线性系统求解及后处理等核心编程环节。压缩包共16个文件核心C源码hpp/cpp配Gmsh几何文件geo实现网格划分结果图片png展示单元剖分与数值解LaTeX与Markdown文档、PDF报告便于对照阅读与编译复现整体仅956KB已有307人学习下载。代码采用模块化的Point、Mesh、Element、Geometry、Function等类设计从弱形式构造到全局矩阵组装与后处理可视化均有清晰实现通过分析并运行可掌握有限元编程的一般框架并了解内存优化、非线性与接触问题等进阶技巧适合数值计算课程设计或科研入门者系统参考。1. 有限元方法课程作业代码把纸面推导变成可信数值结果的那一步有限元方法这门课里推导单刚矩阵是纸上的功夫但2020-2021学年春夏学期浙江大学有限元方法课程作业把问题推到了下一层网格给好、材料给好、载荷和约束给好你要输出节点位移和单元应力。没有代码这些数值就只能靠手算死磕有了代码你能在几分钟内改网格密度、换材料参数、甚至把平面应力改成平面应变重跑一遍。这个标题要解决的正是从“理解公式”到“能跑出结果”之间的那一步。适合正在补交类似作业的本科生、需要带学生的研究生以及在工作里想快速验证有限元概念、又不愿意为一个小问题开一整套商业软件的工程师。下文代码基于 Python 3.8依赖 numpy、scipy、matplotlib一个普通笔记本就能跑完所有算例。2. 有限元方法课程作业代码的核心单刚矩阵与总刚组装顺序写代码之前先确定整个求解流程。有限元平面问题的稳态线性求解本质上就是四步逐单元算单刚、按自由度编号组装总刚、引入载荷与约束后解方程、再回代求应力。课程作业代码写不顺的绝大多数死在第二步编号上而不是死在公式上。下面的示例全部用常应变三角形单元CST原因是它在平面问题作业里出现频率最高、也最容易手算验证理解完这一套骨架换任意单元只是替换单刚计算函数组装和求解流程不动。2.1 平面应力CST单元的B矩阵与D矩阵两个矩阵决定一切CST 单元的形函数是线性的单元内应变处处相等所以单元刚度矩阵不需要数值积分直接用三角形面积乘被积函数。实现只需要两个矩阵几何矩阵 B 把节点位移映射成单元应变弹性矩阵 D 把应变映射成应力两者乘积的二次型就是单刚。平面应力假设下D 矩阵为D E/(1-v²) · [[1, v, 0], [v, 1, 0], [0, 0, (1-v)/2]]这里应变量排列顺序是 (εx, εy, γxy)D 的行列必须跟这个顺序一一对应。如果题目要的是平面应变E 和 v 要按 E/(1-v²) 与 v/(1-v) 替换。作业里最容易在这里翻车同样代码、同样网格只是条件从 plane stress 改成 plane strain位移数值就整体变掉先别怀疑程序先看题目写的是哪种条件。B 矩阵由形函数对坐标的偏导数组成。对单元内第 i 个节点b_i (y_j - y_m) / (2A)c_i (x_m - x_j) / (2A)其中 A 是三角形面积下标 j、m 是单元内另外两个节点顺序必须跟单元连接表一致。计算中每个量的检查点可以直接记成一张表量计算式检查点单元面积A 叉积 / 2必须大于零否则节点顺时针排列b_i(y_j - y_m) / (2A)下标与节点号严格对应c_i(x_m - x_j) / (2A)差一个负号是常见事故应变排列εx, εy, γxy与D矩阵的行顺序一致2.2 Python 实现单刚计算逆时针检查与单位统一放在入口单刚函数不应当只返回一个数组它还应该承担输入校验。下面这段代码用 NumPy 完成全部运算入参只有坐标和材料参数没有全局状态方便单独验证import numpy as np def cst_stiffness(coord, E, nu, t1.0): 平面应力常应变三角形单元单刚矩阵。 coord: (3, 2) 数组三节点坐标逆时针排列。 E: 弹性模量nu: 泊松比t: 厚度。 x coord[:, 0] y coord[:, 1] # 叉积的两倍就是有向面积 A2 (x[1]-x[0])*(y[2]-y[0]) - (x[2]-x[0])*(y[1]-y[0]) if A2 0: raise ValueError(节点必须逆时针排列且面积不能为零) # 形函数导数b_i (y_j - y_m)/2Ac_i (x_m - x_j)/2A bx np.array([y[1]-y[2], y[2]-y[0], y[0]-y[1]]) / A2 cy np.array([x[2]-x[1], x[0]-x[2], x[1]-x[0]]) / A2 B np.zeros((3, 6)) for i in range(3): B[0, 2*i] bx[i] # εx ∂u/∂x B[1, 2*i1] cy[i] # εy ∂v/∂y B[2, 2*i] cy[i] # γxy 中 ∂u/∂y 项 B[2, 2*i1] bx[i] # γxy 中 ∂v/∂x 项 # 平面应力弹性矩阵行序与应变排列一致 D E / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1-nu)/2] ]) ke (t * A2 / 2) * B.T D B return ke逻辑说明A2 是两倍面积如果为负或为零说明单元节点顺序反了或三点共线这类单元一旦进入组装总刚必然奇异bx、cy 各三个值分别对应单元内三个节点它们和 B 矩阵列顺序的对应关系是这段代码里最需要盯住的地方节点 0 的自由度占第 0、1 列节点 1 占第 2、3 列节点 2 占第 4、5 列。参数说明t 默认 1.0平面应力问题里厚度只做整体缩放不影响各点位移的相对分布E 和 nu 与几何、载荷共用一套单位制常见的是 mm-N-MPa 或 m-Pa混用会导致结果整体偏差 10^3 或 10^6 倍这是作业里最隐蔽的一类错误。验证这个函数是否写对有一个不用手算 6×6 矩阵的土办法令 E1、nu0给任意三角形坐标检查 ke 每行之和应当接近零这对应刚体平动再取 (u, v)(-y, x) 作为节点位移向量ke 乘上它应当也接近零这对应刚体转动。组装之后未加约束的总刚应该有 3 个接近零的特征值多于 3 个就要怀疑网格里有悬空节点。2.3 组装总刚自由度编号是唯一的对应关系总刚组装是一个“散开再叠加”的过程。单元刚度矩阵是 6×6散到总刚里时每个数值要放到全局自由度索引指定的位置上。这里唯一的对应关系就是节点编号到自由度编号的映射。下面这段代码用固定规则节点 i 的 x 方向自由度是 2iy 方向是 2i1。def assemble(K_list, elems, nnodes): 把单元单刚叠加成总刚。 K_list: 每个单元的单刚矩阵列表。 elems: (nelem, 3) 整数数组每行是单元的三个节点编号。 nnodes: 节点总数。 n 2 * nnodes K np.zeros((n, n)) for ke, e in zip(K_list, elems): idx np.array([ 2*e[0], 2*e[0]1, 2*e[1], 2*e[1]1, 2*e[2], 2*e[2]1, ]) K[np.ix_(idx, idx)] ke return K参数说明K_list 的顺序必须和 elems 的行顺序一致zip 按位置配对两个列表排序不一致时所有单元刚度都会装错位置elems 里节点编号约定从 0 开始跟网格生成函数的编号方式统一常见事故是从 MATLAB 习惯改到 Python 时忘记减一nnodes 只用来确定总刚尺寸真正的索引逻辑全在 idx 里。叠加用而不是因为公共节点被多个单元共享贡献必须累加。组装完成后跑一句 np.allclose(K, K.T) 检查对称性再检查对角线元素全部大于零这两个检查能挡住九成以上的索引错误。提示总刚对角线出现零几乎都是网格中存在未被任何单元引用的孤立节点。检查单元连接表里每个节点号是否都被用到。3. 载荷向量与边界条件决定解长什么样的两个入口刚度矩阵本身只决定“解空间的样子”把问题变成唯一解的是右侧载荷向量和位移约束。作业代码在这里最容易变乱有人把约束硬编码在主脚本里有人把集中力写死进网格函数。我的建议是先规划两个接口一个负责把边分布力和集中力转成节点力一个负责处理给定位移保证主流程里只出现一次组装和一次求解。3.1 集中力、体力和边分布力的节点等效集中力最好处理力作用在某个节点上直接在全局载荷向量 F 的对应自由度位置加值。麻烦的是作用在单元边界上的分布力。对 CST 单元形函数在边上线性分布边上均布面力的等效节点力积分结果就是两端节点各承担一半。下面这个函数把边上的均布载荷转成节点力def add_edge_load(F, nodes, node_a, node_b, qx, qy): 向全局载荷向量添加边界均布力。 F: 全局载荷向量。 nodes: 全部节点坐标。 node_a, node_b: 受力边两端节点编号。 qx, qy: 该边上均布力在 x、y 方向的单位长度合力分量。 xa, ya nodes[node_a] xb, yb nodes[node_b] L np.hypot(xb - xa, yb - ya) if L 1e-12: return F F[2*node_a] qx * L / 2 F[2*node_a 1] qy * L / 2 F[2*node_b] qx * L / 2 F[2*node_b 1] qy * L / 2 return F逻辑说明线性形函数在边上的积分结果是固定的 1/2 分配这个结论只在单元边上有效如果载荷作用在单元内部就要改用 f_e ∫ Nᵀ b dACST 单元的数值结果是把体力总量按 1/3 分摊给三个节点。参数说明qx、qy 的单位要和弹性模量配套用 mm-N-MPa 体系时是 N/mm用 m-Pa 体系时是 N/m。单位不一致导致的错误有个明显特征位移结果整体偏大或偏小同样的倍数且与网格密度无关只有量纲问题会这样“错得均匀”。3.2 给定位移约束删行删列、罚函数与划零置一的取舍位移约束有三种常见实现不需要都写选一个贯彻到底。第一种删行删列把被约束的自由度从总刚和载荷中移除得到缩减后的 Kff 与 Ff解出来就是自由位移之后再把完整位移向量补回去。数值最稳定条件数最好我一般默认用它。第二种罚函数把被约束自由度对应的主对角元素乘一个大数再把载荷项置为约束值与罚值的乘积实现简单但罚值取 1e8 还是 1e12 会影响精度取太大又恶化条件数适合演示工程上不推荐。第三种划零置一对应行和列清零、对角线置 1实现简单但破坏矩阵结构节点数上千后求解明显变慢。删行删列的完整实现def apply_prescribed(K, F, prescribed): 删行删列处理给定位移。 prescribed: dict{全局自由度编号: 给定位移值}。返回缩减后的 Kff、Ff 和自由自由度列表 free。 fixed list(prescribed.keys()) free [i for i in range(K.shape[0]) if i not in fixed] Kff K[np.ix_(free, free)] Ff F[free].copy() for dof, val in prescribed.items(): if val ! 0.0: # 给定位移对自由自由度的耦合贡献搬到右端 Ff - K[np.ix_(free, [dof])].flatten() * val return Kff, Ff, free逻辑说明prescribed 的键是被约束的自由度编号值是对应位移零位移最常见。非零位移时被约束自由度与自由自由度之间的耦合项 K[np.ix_(free, [dof])] 乘上位移值要从右端减去这是最容易漏的一步。free 列表后续用来回填新建一个长度为 2*nnodes 的零数组free 位置放 u_freefixed 位置放 prescribed 里对应的值就得到完整位移向量。3.3 求解环节稠密求解、稀疏求解与奇异排查求解本身不难选对方式就行。自由度数低于一千时直接用 np.linalg.solve再往上Kff 是稀疏的却用稠密存储内存和耗时都涨得快三五千自由度以后换成 scipy 稀疏求解from scipy.sparse import csr_matrix from scipy.sparse.linalg import spsolve u_free spsolve(csr_matrix(Kff), Ff)选用 spsolve 后存储从 (n, n) 稠密数组变成只存非零元的 CSR 矩阵典型平面网格上求解速度快一到两个数量级。课程作业动辄几十乘几十的网格自由度过万很常见换这一行就能把时间从分钟级降到秒级。下面这个经验表可以直接当起点自由度数推荐路径备注 1000np.linalg.solve最稳省去稀疏格式转换1000 ~ 10000spsolve(csr_matrix(Kff))稀疏收益明显 10000先检查网格质量再考虑局部细化或换单元求解前先确认约束足以消除刚体位移。约束不够时求解器常常不直接报错而是返回巨大数值。看到结果异常先回头检查 prescribed 是否覆盖了所有刚体位移再谈算法。4. 完整算例平面应力悬臂梁的网格、求解与应力输出前面三章的函数全部就位后主流程就是组合调用。这一章用一个矩形悬臂梁算例把这些代码串起来左端固支右端面中点受向下集中力材料用结构钢参数。这个算例能跟梁理论解对照每一处参数错误都会在位移量级和分布形状上露出马脚。4.1 矩形网格生成横竖切分与单元方向控制网格生成函数只需要做一件事把 [0, width] × [0, height] 矩形区域按 nx × ny 切分再用对角线把每个小矩形切成两个三角形。三角形方向必须是逆时针否则 2.2 节的面积检查会直接报错。下面函数直接返回节点数组和单元连接表def rect_mesh(nx, ny, width, height): 生成矩形区域的三角形网格。 nx, ny: x、y 方向的分段数。 width, height: 矩形宽度和高度。返回 nodes: (n, 2) 坐标数组elems: (m, 3) 单元连接表。 nno (nx 1) * (ny 1) nodes np.zeros((nno, 2)) for j in range(ny 1): for i in range(nx 1): nodes[j*(nx1) i] [i*width/nx, j*height/ny] elems [] for j in range(ny): for i in range(nx): a j*(nx1) i b a 1 c a (nx1) d c 1 elems.append([a, b, c]) # 左下、右下、左上逆时针 elems.append([b, d, c]) # 右下、右上、左上逆时针 return nodes, np.array(elems)参数说明nx 和 ny 决定网格粗细。悬臂梁以弯曲为主长度方向分密些、厚度方向分少些也能接受但应力集中区域要局部加密。每个小矩形切成的两个三角形共享一条对角线对角线方向一致时应力云图上可能出现斜向纹理这是网格各向异性的表现不是代码错。追求更均匀结果的做法是每行交错改变对角线方向但课程作业里先保持一致方向也够用。4.2 主程序流程五次调用串起完整求解主程序把所有函数按顺序串起来长度控制在五十行以内# 参数单位制用 m - N - Pa E, nu, t 210e9, 0.3, 0.01 # 结构钢厚度 10mm L, H 0.8, 0.1 # 梁长 800mm高 100mm nx, ny 40, 10 nodes, elems rect_mesh(nx, ny, L, H) nnodes len(nodes) # 1. 组装总刚 K np.zeros((2*nnodes, 2*nnodes)) for e in elems: ke cst_stiffness(nodes[e], E, nu, t) idx np.array([2*e[0], 2*e[0]1, 2*e[1], 2*e[1]1, 2*e[2], 2*e[2]1]) K[np.ix_(idx, idx)] ke # 2. 载荷右端面中点施加向下集中力 1000 N P -1000.0 load_node (ny // 2) * (nx 1) nx # inx, jny//2 F np.zeros(2*nnodes) F[2*load_node 1] P # 3. 约束左端 x0 上所有节点 6 个自由度全固定 prescribed {} for j in range(ny 1): node j * (nx 1) prescribed[2*node] 0.0 prescribed[2*node 1] 0.0 # 4. 求解 Kff, Ff, free apply_prescribed(K, F, prescribed) u_free np.linalg.solve(Kff, Ff) # 4000 自由度以内直接用 u np.zeros(2*nnodes) u[free] u_free # 5. 输出端部加载点的竖向位移 print(端部加载点竖向位移:, u[2*load_node 1])主流程一共五步组装、载荷、约束、求解、取结果。组装处的 nodes[e] 用 NumPy 花式索引一次取出单元三个节点的坐标矩阵正好是 cst_stiffness 期望的 (3,2) 输入载荷节点必须是网格内部节点左端约束必须同时固定 x、y 两个方向否则梁会刚体转动np.linalg.solve 在自由度数约一千以内安全更大就把求解替换成前面的 spsolve。这个算例中用集中力直接跟梁理论对比会比较粗糙因为集中力会在加载点附近产生应力集中弹性解在那里奇异。更标准的做法是把右端面所有节点各分一部分力近似成均匀剪力或者把加载点往里挪一个截面高度。课程作业如果要跟解析解对标优先用均布端面力替代集中力。4.3 参数调试与失败排查症状、原因和排查顺序症状可能原因排查顺序矩阵奇异或解出 1e15 量级位移刚体位移未约束悬空节点先数约束个数再查连通性位移分布的对称性被破坏载荷方向写错约束不对称检查 F 符号检查 prescribed 的键应力云图跨单元跳变明显网格太粗CST是常应变单元加密网格或换高次单元结果与解析解整体差一个倍数单位制不统一对 E、t、载荷做量纲核对排查顺序的固定经验先单位再约束再载荷方向最后才怀疑刚度矩阵。组装完成后打印 np.linalg.norm(K)约束后打印 np.linalg.cond(Kff)能快速定位问题出在哪一步。如果把单元换成四节点等参元单刚里会多出雅可比矩阵和 2×2 高斯积分但组装、约束、求解流程完全一致积分阶数取低时会出现零能模式表现为位移结果带着网格尺度的锯齿这一点在做四边形单元作业时提前留意。5. 收敛性验证与结果可视化让作业结果禁得起追问5.1 用解析解做定量收敛对照悬臂梁端部受均布端面剪力时欧拉-伯努利梁理论的端部位移为 w P L³/(3EI)其中 I t H³/12。梁理论是近似但细长比大于 10 时平面应力解和它差别很小足以当量级尺子。把 4.2 的主程序放进循环nx 依次取 10、20、40、80记录端部位移与解析解误差。这时会看到一个现象误差随网格加密下降但降到某个程度后不再降停在梁理论与平面弹性解的模型差上。这不是代码 bug而是解析解模式本身的限制在报告里写清楚这一条比多加密一次网格更说明问题。CST 单元的位移按能量范数是一阶收敛误差大致随单元尺寸 h 线性下降。四组数据画成 log-log 图斜率应接近 1。斜率明显小于 1优先怀疑网格里有畸形单元曲线完全不下降回 4.3 查约束和单位。5.2 用 matplotlib 输出变形图与应力云图matplotlib 自带三角网格绘制不需要额外库。变形图把节点坐标加上位移乘以放大系数再画应力云图用 tricontourfimport matplotlib.pyplot as plt import matplotlib.tri as mtri # 变形图放大系数 scale 让变形肉眼可见 scale 100.0 deformed nodes scale * u.reshape(-1, 2) plt.triplot(deformed[:, 0], deformed[:, 1], elems, lw0.3) # 应力云图stress_nodal 是节点平均应力 triang mtri.Triangulation(nodes[:, 0], nodes[:, 1], elems) plt.tricontourf(triang, stress_nodal, levels20, cmapviridis)坑点有两个。位移量级常在 0.01mm 以下不放大直接画会和原网格重叠scale 至少要取 100tricontourf 需要节点上的值而 CST 应力是单元常数直接传单元应力会导致相邻单元色块跳变处理办法是先把单元应力算出来再按共享节点取平均得到节点值。平均会抹平角点应力报告里要说明这一点。5.3 整理作业代码时的两个习惯第一个习惯是把参数集中放在文件开头E、nu、t、几何尺寸、载荷数值、网格分段全部排成一组常量或一个 dict注释里写明单位。这样换算例只改一块区域不用在函数堆里翻。第二个习惯是给每个函数留一句 docstring写明输入的单位和顺序尤其是 cst_stiffness 的 coord 行序和 assemble 的 elems 编号起始。两周后把这种注释补全的代码翻出来5 分钟内跑通并解释每一行比记住任何算法结论都更快进入下一个算例。本文还有配套的精品资源点击获取

Python实现有限元作业代码：CST单元单刚组装与求解全流程

相关推荐

长沙跨境电商静态页开发：HTML5语义化+CSS响应式+本地JS交互

Python+OpenCV图像处理实战：从基础算法到应用案例

Vue3+TS舆情监测系统模板：实时词云、权限路由与数据清洗实战

更多资讯

Multisim 14.3 安装配置指南：数据库报错与元件库修复全攻略

JavaWeb购物商城项目：原生Servlet+JSP+MySQL全链路实战

Keil MDK 5.39嵌入式开发实操指南：稳定适配STM32与国产MCU

IAR嵌入式开发全链路解析：编译、链接与SWO调试深度实践

免疫算法求解配送中心选址问题的MATLAB实现

拆解大众ID. Era 9X座舱与智驾域控制器，读懂德系稳健的电子电气架构

今日更新

51单片机开环控制磁阻传感器的硬件匹配与代码实现

AI Agent安全操作数据库：SQL注入防御与参数化查询实战

Android逆向实战：还原激励类App签到接口的签名算法

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？