4. 二维平面应力三角形单元刚度矩阵组装
问题对象与离散化
取一块 200 mm × 100 mm、中心带 φ40 mm 圆孔的平板，厚度 t = 5 mm，材料 E = 210 GPa、ν = 0.3，左端固定，右端受均布拉力 q = 100 N/mm。采用常应变三角形单元（CST）离散，每个节点有 u、v 两个自由度。CST 单元内应变恒定，因此网格越密，孔边应力集中越接近真实值。
单元刚度矩阵由下式给出：
k = t · A · Bᵀ · D · B其中 A 为三角形面积，B 为 3×6 应变-位移矩阵，D 为平面应力弹性矩阵：
D = E/(1-ν²) · [[1, ν, 0],
[ν, 1, 0],
[0, 0, (1-ν)/2]]单元几何与 B 矩阵
设三角形节点按逆时针编号 1、2、3，坐标为 (x₁,y₁)、(x₂,y₂)、(x₃,y₃)。定义：
b₁ = y₂ - y₃, b₂ = y₃ - y₁, b₃ = y₁ - y₂
c₁ = x₃ - x₂, c₂ = x₁ - x₃, c₃ = x₂ - x₁
A = ½ |b₁c₂ - b₂c₁|B 矩阵为：
B = 1/(2A) · [[b₁, 0, b₂, 0, b₃, 0],
[0, c₁, 0, c₂, 0, c₃],
[c₁, b₁, c₂, b₂, c₃, b₃]]组装整体刚度矩阵
整体刚度矩阵 K 的维度为 2N × 2N（N 为节点总数）。对每个单元，将其 6×6 的 k 按节点自由度映射到 K 的对应位置并累加。以右端受拉为例，载荷向量在右端节点 v 方向施加 q·t·(节点间距) 的等效节点力。
边界条件：左端所有节点 u = v = 0。采用置大数法或行列划去法处理，求解 K·U = F 得到节点位移，再由 σ = D·B·uᵉ 计算每个单元的应力。
关键实现要点
- 节点编号与自由度映射：节点 i 的 u 对应全局自由度 2i-1，v 对应 2i。单元节点 (n₁,n₂,n₃) 的 6 个局部自由度依次映射到 [2n₁-1, 2n₁, 2n₂-1, 2n₂, 2n₃-1, 2n₃]。
- 面积符号：用行列式计算 2A = b₁c₂ - b₂c₁，取绝对值保证面积为正；若节点逆时针排列，行列式为正。
- 组装累加：K[I,J] += k[a,b]，其中 I、J 为全局自由度编号，a、b 为局部编号。多个单元共享节点时刚度自动叠加。
- 边界条件：左端节点 u=v=0，采用行列划去法将对应行列从 K 中移除后再求解，避免置大数法引入数值误差。
- 应力恢复：每个单元内 σ = D·B·uᵉ 为常值，孔边单元应力显著高于远场，需加密网格捕捉梯度。
CST 单元在孔边弯曲或应力梯度大的区域精度有限，建议在孔周布置更密的三角形网格，或改用后邻例中的四边形等参单元。