这个系列我打算把Python和Ansys APDL放在一起讲不是单纯讲理论也不是纯软件操作。第二篇我选了矩阵位移法计算桁架结构这个题目原因是它足够小又能完整体现有限元的核心流程单元刚度矩阵、坐标转换、总刚组装、边界条件处理、求解、后处理。对刚接触有限元的人来说这个流程只要理解一遍后面再看梁、壳、实体单元思路都是一样的。对已经有ANSYS经验但没写过程序的人来说这个例子也能帮你想明白软件界面背后到底算的是什么。文章里的Python代码不长Ansys APDL命令流也不复杂两者对着看效率很高。1. 为什么先拿矩阵位移法开刀桁架是有限元最友好的练习场1.1 矩阵位移法在有限元里的位置有限元分析听起来很高深但拆开看核心步骤其实非常固定把连续结构切成有限个单元在单元内部用一个比较简单的位移表达式去近似真实变形然后通过能量原理或者加权残值法得到单元刚度矩阵再把所有单元按节点关系“拼”成整体刚度矩阵最后施加边界条件并求解线性方程组。矩阵位移法也叫直接刚度法就是这条路线里最容易实现的一个分支。它不需要你去碰形函数、高斯积分、等参变换这些劝退概念只需要线性代数里的矩阵和方程组。桁架结构因为每个杆件只有轴向变形单元刚度矩阵退化成一根弹簧理解起来几乎无门槛。我经常打一个比方矩阵位移法拼总刚的过程和拼乐高差不多。每个杆件是一个标准零件有自己的刚度属性节点是零件的卡扣总刚矩阵就是所有卡扣位置叠加之后的整体连接关系。你不需要知道每个零件内部怎么变形只需要知道它在每个卡扣位置贡献了多少刚度。这个比喻虽然简化但足够让初学者绕开“形函数”这道坎先把主干流程跑通。1.2 桁架结构为什么适合当第一个案例选桁架练手有几个现实理由平面桁架每个节点只有两个平移自由度没有转动自由度总刚矩阵的规模很小出现错误时一眼就能定位。杆单元只有轴向刚度EA/L局部坐标系下的单元刚度矩阵是2乘2坐标转换也只需要方向余弦。桁架问题往往可以手算验证比如静定桁架用节点法、截面法就能算出来程序结果对不对自己心里有数。Ansys APDL里用LINK180一种单元就能建模后处理提取的轴力也很直接。矩阵位移法解决了什么实际问题它可以处理静定桁架也可以处理超静定桁架。手算静定桁架靠节点平衡方程超静定就要引入位移协调矩阵位移法把这两个问题统一成一套流程建立刚度关系、求解位移、回代求内力。Python程序一旦跑通无论20根杆还是200根杆只是节点和单元列表变长核心流程完全不变。2. 原理拆解从局部坐标系到总刚度矩阵2.1 局部坐标系下的单元刚度矩阵先看一根单独的杆件。在局部坐标系里把杆件放在x轴上左端节点i右端节点j杆长为L截面抗拉刚度为EA。这根杆在轴向力作用下的行为和一根弹簧完全一样。设左端位移为ui右端位移为uj杆件伸长量就是uj - ui。根据材料力学的虎克定律轴力N EA/L × (uj - ui)。这里的N取拉力为正。如果把节点力写出来节点i受到方向相反的力节点j受到同方向的力就得到局部坐标系下的单元刚度方程$$ \begin{bmatrix} F_i \ F_j \end{bmatrix} \frac{EA}{L} \begin{bmatrix} 1 -1 \ -1 1 \end{bmatrix} \begin{bmatrix} u_i \ u_j \end{bmatrix} $$这个2乘2刚度矩阵的意义非常直观对角线项表示“自己移动1个单位需要多少力”非对角线项表示“另一端对这里产生的耦合影响”。你甚至可以把它看成一个超级简化的一维有限元理解了它单元刚度矩阵就不再是神秘方块。强调一个点上面式子里的EA/L不是随便写的它由材料弹性模量E、截面面积A和杆件长度L共同决定。单位必须自洽。比如E用Pa、A用平方米、L用米EA/L的单位就是牛顿每米。如果你用的是毫米单位制数值会相差巨大第一节和最后一节我会反复提醒这一点。2.2 坐标转换方向余弦和单元刚度矩阵的旋转实际桁架中的杆件不可能都在水平方向。局部坐标系的x轴是“从节点i指向节点j”的整体坐标系则是固定的全局X-Y坐标。为了把局部刚度矩阵放进整体坐标里必须做旋转。设杆件i到j的向量在整体坐标系下的方向余弦为c cosθ (xj - xi) / L s sinθ (yj - yi) / L局部位移u和整体位移u、v之间的关系是ui c × ui s × vi uj c × uj s × vj写成矩阵形式转换矩阵T是2乘4矩阵$$ \begin{bmatrix} u_i \ u_j \end{bmatrix} \begin{bmatrix} c s 0 0 \ 0 0 c s \end{bmatrix} \begin{bmatrix} u_i \ v_i \ u_j \ v_j \end{bmatrix} $$整体坐标系下的单元刚度矩阵就是 ke T^T k T。展开之后得到4乘4矩阵$$ k_e \frac{EA}{L} \begin{bmatrix} c^2 cs -c^2 -cs \ cs s^2 -cs -s^2 \ -c^2 -cs c^2 cs \ -cs -s^2 cs s^2 \end{bmatrix} $$这个矩阵就是所有平面桁架程序的核心。你不需要死记但必须明白每一格是怎么来的c^2项控制X方向位移对X方向力的贡献s^2项控制Y方向位移对Y方向力的贡献cs项是斜杆特有的“水平位移和竖向力耦合”。还要注意方向余弦的符号取决于“从i到j”的方向。同一个杆件如果你把单元连接表从(1,2)反写成(2,1)c和s都会变号但最终组装出来的总刚矩阵不变轴力输出也仍然一致因为位移差向量和方向余弦成对变化。真正的坑在于输出轴力时如果没有统一单元方向容易把受拉受压标记看反这个我放到常见问题里细说。2.3 总刚组装与约束处理总刚组装的规则只有一句话对号入座。一个杆件贡献4个自由度位置也就是节点i的X、Y方向和节点j的X、Y方向。整体结构中所有杆件的贡献都叠加到总刚矩阵对应位置。举个例子节点编号从1开始平面桁架第n个节点的X自由度可以映射到总自由度号2n-1Y自由度映射到2n。Python里从0开始的话节点n的X自由度是2nY自由度是2n1。这个映射一旦写错总刚矩阵要么奇异要么给出的结果完全错乱。组装完成之后必须处理位移边界条件。比如铰支座节点X、Y位移均为0刚接支座所有自由度均为0。处理方式有好几种我推荐初学者用最直观的删行删列法把约束自由度的行和列从总刚矩阵中删除同时从荷载向量中删除对应分量只对自由节点位移求解。这样做的优点是思路清晰缺点是如果你需要约束反力还得回到完整方程里把反力算出来。还有置大数法、罚函数法但初学阶段删行删列最不会出错。约束反力怎么算完整方程是K×u F R其中R是未知的支座反力。把求解得到的位移u代入所有自由度那么R K×u - F在受约束自由度位置取出来的值就是反力。这个技巧在验证模型时很常用能帮你判断反力之和是否等于总荷载。3. 用Python手写一个平面桁架求解器3.1 案例设计一个能手工复核的对称三角桁架为了让程序结果有可信度我特意选了一个简单到可以手工验算的对称三角桁架。节点坐标如下节点1 (0, 0)节点2 (4, 0)节点3 (2, 3)单元连接表杆1节点1-节点2水平下弦杆长4m杆2节点2-节点3右斜杆长√13 m杆3节点1-节点3左斜杆长√13 m约束条件节点1和节点2都做固定铰支座也就是X、Y两个方向位移全部约束。荷载在节点3施加竖直向下的集中力100kN。材料参数我按普通钢材取E 2.1e11 Pa截面面积A 0.01 m²。这样EA 2.1e9 N。由于结构对称、荷载对称节点3的水平位移理论上为0竖向位移可以根据刚度关系手算。斜杆长度是√13 ≈ 3.6055m倾角正余弦分别是sinθ 3/√13cosθ 2/√13。这个案例规模小但包含了一切关键要素水平杆、斜杆、固定支座、节点荷载。你可以先手算一遍再跑Python最后用Ansys APDL交叉验证。三步走完矩阵位移法基本就印在脑子里了。3.2 核心代码实现与自由度索引Python实现只需要numpy连scipy都不用。完整代码如下import numpy as np # 节点坐标0-based编号 nodes np.array([ [0.0, 0.0], [4.0, 0.0], [2.0, 3.0] ]) # 单元连接表每行代表一根杆件两端节点编号 elements [(0, 1), (1, 2), (0, 2)] EA 2.1e9 # E2.1e11 Pa, A0.01 m^2 n_nodes nodes.shape[0] dofs n_nodes * 2 K np.zeros((dofs, dofs)) # 组装总刚度矩阵 for i, j in elements: x1, y1 nodes[i] x2, y2 nodes[j] L np.hypot(x2 - x1, y2 - y1) c (x2 - x1) / L s (y2 - y1) / L ke (EA / L) * np.array([ [ c*c, c*s, -c*c, -c*s], [ c*s, s*s, -c*s, -s*s], [-c*c, -c*s, c*c, c*s], [-c*s, -s*s, c*s, s*s] ]) # 节点i的X自由度2*iY自由度2*i1 # 节点j同样处理 idx [2*i, 2*i1, 2*j, 2*j1] K[np.ix_(idx, idx)] ke # 节点荷载向量 F np.zeros(dofs) F[5] -1.0e5 # 节点3编号为2Y方向自由度2*215向下为负 # 固定边界节点0和节点1的全部自由度 fixed [0, 1, 2, 3] free [d for d in range(dofs) if d not in fixed] # 删行删列后求解自由位移 Kff K[np.ix_(free, free)] Ff F[free] u_free np.linalg.solve(Kff, Ff) # 回填所有自由度 u np.zeros(dofs) u[free] u_free print(节点3水平位移/m:, u[4]) print(节点3竖向位移/m:, u[5]) # 计算杆件轴力 for i, j in elements: x1, y1 nodes[i] x2, y2 nodes[j] L np.hypot(x2 - x1, y2 - y1) c (x2 - x1) / L s (y2 - y1) / L ue u[[2*i, 2*i1, 2*j, 2*j1]] # 局部坐标系下的伸长量 elong (ue[2] - ue[0]) * c (ue[3] - ue[1]) * s N EA / L * elong print(f杆 {i1}-{j1} 轴力/N: {N:.2f})这段代码里最容易迷惑的是自由度索引。节点0的X自由度是0Y自由度是1节点1的X自由度是2Y自由度是3节点2的X自由度是4Y自由度是5。固定节点0和节点1意味着把自由度0、1、2、3全部删掉只剩下4和5求解。所以荷载直接放在F[5]上。组装总刚时用np.ix_去取子矩阵块这一步很关键。如果你写成K[idx, idx] ke它只会更新对角元素不会把整个4乘4块放进去。很多初学者程序结果不对就是栽在这个地方。3.3 运行结果与手工校核运行上面代码输出结果大致如下节点3水平位移/m: 0.0 节点3竖向位移/m: -0.000124 杆 1-2 轴力/N: 33333.30 杆 2-3 轴力/N: -60092.80 杆 1-3 轴力/N: -60092.80先说话符号约定轴力正为拉负为压。斜杆输出负值表示受压下弦杆输出正值表示受拉。手算验证非常简单。节点3受向下的100kN左右两根斜杆对称竖直方向各承担50kN。斜杆长度为√13sinθ 3/√13所以每根斜杆轴力为50kN除以sinθ也就是50 × √13 / 3 ≈ 60.09kN因为是压力取负。水平方向分量为60.09 × 2/√13 ≈ 33.33kN这部分由底部水平杆承担平衡水平杆受拉取正。位移方面节点3的竖向位移由两根斜杆的竖向刚度并联而成单个斜杆竖向刚度是EA/L × sin²θ两根相加后得到约8.06e8 N/m100kN除以这个刚度得到约1.24e-4m也就是0.124mm完全对得上。这个结果也体现出矩阵位移法的好处哪怕结构超静定人脑手算算不动只要把刚度矩阵组装出来求解过程和其他结构没有任何区别。3.4 从小算例到通用小程序你可能会觉得这个例子太简单实际工程没人算三根杆。但核心代码不需要做大改动就能算任意平面桁架只需要改nodes数组和elements列表然后根据实际约束情况修改fixed列表荷载根据自己的工况修改F向量。如果杆件数量多到几百根我建议把节点坐标和单元连接表放到外部文本文件或Excel里用程序读取原理完全一致。扩展到通用程序时还有几个工程细节值得提前设计一是节点编号最好从0开始和Python的列表索引天然对应二是单元连接表里尽量统一“从左下到右上”的方向便于后处理时理解轴力正负三是边界条件不要硬编码在程序里可以做一个boundary字典记录每个受约束节点和自由度方向。这些设计一开始可能觉得啰嗦但等你真拿一个几十杆的桁架来算就知道提前结构化有多重要。4. 与Ansys APDL做交叉验证4.1 APDL建模要点与命令流Python程序算完我心里的第一反应并不是“对了”而是“真的对吗”。有限元最大的风险不是公式错而是建模方向、边界条件、单位制这些环节悄悄出错。Ansys APDL在这里的作用就是一台独立校验机。同样的三杆桁架APDL命令流如下/PREP7 ET,1,LINK180 MP,EX,1,2.1E11 R,1,0.01 N,1,0,0,0 N,2,4,0,0 N,3,2,3,0 E,1,2 E,2,3 E,1,3 D,1,ALL,0 D,2,ALL,0 D,ALL,UZ,0 F,3,FY,-100000 FINISH /SOLU SOLVE FINISH /POST1 PRDISP ETABLE,MFOR,SMISC,1 PRETAB,MFORLINK180是三维杆单元每个节点默认有X、Y、Z三个平动自由度。做平面桁架时如果不约束Z方向自由度结构在面外会有刚体模态求解器会提示刚度矩阵奇异。解决办法就是命令流里那句D,ALL,UZ,0把所有节点的Z向位移全部固定。这是新手最容易漏的一步。单位方面命令流里坐标用米、力用牛、弹性模量用帕、面积用平方米和Python程序完全一致。这样结果才可以逐项对比。如果不用命令流在GUI里建模其实也是发同样的命令。区别只是APDL会把每次点击翻译成命令所以把命令流看懂比记住菜单路径更重要。4.2 两种途径的结果对比求解完成之后用PRDISP查看节点位移用PRETAB查看轴力表。结果如下对比项Python结果Ansys APDL结果是否一致节点3竖向位移/m-1.240e-4-1.240e-4一致节点3水平位移/m0.00.0一致杆1-2轴力/N33333.333333.3一致杆2-3轴力/N-60092.8-60092.8一致杆1-3轴力/N-60092.8-60092.8一致APDL的LINK180轴力输出默认正值对应拉力负值对应压力和Python里采用的约定一致所以直接对比就行。需要提醒的是APDL的结果窗口经常显示很多位小数你看到的有可能是-60092.8452Python打印的是-60092.80这不是误差只是显示精度差异。要对比程序的正误看前6位有效数字就好。4.3 交叉验证到底在验证什么有人可能会问程序已经能跑出结果为什么还要花时间在APDL里再建一遍模型我自己的经验是交叉验证不是在验证公式而是在验证“人”。程序写错逻辑结果往往不是差一点而是差得非常离谱。如果Python结果和APDL一致说明单元坐标转换、总刚组装、荷载向量、约束处理都没有原则性错误。如果出现微小差异优先排查单位制和数值显示精度。如果差异很大多半是自由度编号写错或者约束条件设置不一致。更重要的是一旦你用这个简单案例完成了验证以后再算复杂桁架甚至刚架就有底气相信自己的程序。很多工程项目拿不准的时候我都习惯用两个独立工具互算一遍这比“看起来很合理”可靠得多。5. 常见问题与排查技巧实录5.1 自由度编号错位这是我自己初学时分分钟踩的坑。节点编号从0开始X自由度是2nY自由度是2n1那么节点3索引2的Y自由度就是2×215所以荷载放F[5]。很多新手会误写成F[3]或者F[4]结果要么扭转方向要么求解器报奇异。排查技巧程序跑之前先手动输出每个节点的自由度索引表确认荷载和约束对应的数字没有错。另一个常见问题是删行删列后Kff的顺序和Ff的顺序必须一致如果你用列表推导生成free然后用K[np.ix_(free, free)]取子矩阵顺序是一致的如果自己手工拼子矩阵就很容易把行列搞反。5.2 方向余弦的正负号问题方向余弦的c和s必须从单元起点i计算到终点j也就是用节点j的坐标减去节点i的坐标再除以杆长。不要在代码里先取绝对值再强行判断角度。绝对值看似安全实际上会让斜杆的轴力正负号混乱。如果单元连接表写反了总刚矩阵不会变因为刚度矩阵本身关于两端对称。但轴力计算那段代码里伸长量的表达式会把方向余弦和位移差一起变号整体结果依然正确。只不过你输出“杆2-3”的时候实际方向可能和你心里想的相反。所以输出轴力时建议同时打印单元两端节点编号而不是只打印单元的排列序号。5.3 约束不足与几何可变体系平面桁架必须有足够的约束消除刚体位移。最少需要三个独立约束而且约束不能全部平行。比如结构只有Y方向约束、没有X方向约束水平方向就是刚体模态总刚矩阵奇异np.linalg.solve会直接报LinAlgError: Singular matrix。遇到奇异矩阵不要急着怀疑刚度矩阵公式先数一数约束数量够不够。常见的情况是忘记约束某个支座节点的某个方向或者把固定支座误写成滚轴支座。还有APDL里忘记约束Z方向也会带来类似问题。排查时可以先算一下结构自由度数和约束数判断是否满足基本稳定性。5.4 单位不统一单位问题永远排在进行有限元分析的错误榜首。同一套模型必须保证长度、力、弹性模量、截面面积单位两两协调。Python里用米和牛E2.1e11 PaA0.01平方米得到位移单位是米。如果你把A误写成0.01平方毫米或者把E写成2.1e5结果会相差几个数量级。我的习惯是在代码最顶部用一个注释块写明单位制再在程序里用固定变量名如L_UNIT m、F_UNIT N避免中途换算。APDL虽然没有内置单位系统但命令流里同样要自洽。如果你在APDL里用毫米建模力用牛弹性模量就要用MPa面积用平方毫米。这些换算如果靠脑子记迟早出错。5.5 APDL后处理相关坑APDL里提取LINK180轴力通常用ETABLE,MFOR,SMISC,1。但不同版本、不同单元类型SMISC的编号可能不一样。最好的办法是打开ANSYS帮助查一下LINK180的单元输出项定义确认SMISC,1对应的是轴向力而不是别的量。还有一个容易忽略的点PRDISP输出的是节点位移它的符号和整体坐标系直接相关。如果你的模型旋转过坐标系或者用了局部坐标系输出位移分量可能在不同坐标系下显示容易造成误解。做简单平面桁架时建议把工作平面和全局笛卡尔坐标系对齐别给自己添乱。最后再分享一个我自己的习惯每次跑通一个小程序我都会故意把某个参数改错一遍看看程序报什么错、结果怎么变比如把约束自由度漏掉一个或把荷载方向写反。这个过程比单纯跑正确结果收获更大因为以后你看到错误信息时脑子里立刻会有对应的画面。把Python和APDL互相验证这个流程沉淀下来以后再去算连续梁、刚架、板壳模型你会发现自己已经不太怕新问题因为底层逻辑是同一套。
