PyNite DKMQ板单元揭秘：四边形板有限元公式推导详解-育师

PyNite DKMQ板单元揭秘：四边形板有限元公式推导详解

【免费下载链接】PyNiteA 3D structural engineering finite element library for Python.项目地址: https://gitcode.com/gh_mirrors/py/PyNite

PyNite 是一个用 Python 编写的 3D 结构工程有限元库，本文带你完整看懂其 DKMQ 四边形板单元的有限元公式推导过程：从自由度布置、双线性形状函数、离散 Kirchhoff 约束，到弯曲刚度矩阵的组装，每一步都讲透，新手也能轻松跟上 👇

一、为什么 DKMQ 板单元值得深究？

PyNite 提供两种板单元，它们定位不同（详见文档 docs/source/plate.rst 的说明）：

单元类型	几何要求	核心思想	适用场景
`Rect`矩形板	必须为矩形	12 项多项式弯曲函数	矩形网格、快速建模
`Quad`四边形板	任意四边形	DKMQ 等参数化公式	厚薄板通吃、扭曲网格

DKMQ（Discrete Kirchhoff Mindlin Quadrilateral）公式的精髓在于：把 Kirchhoff 板理论的高精度弯曲行为，与 Mindlin 理论的横向剪切变形"嫁接"在一起。这就让它既不失薄板的精度，又不会在厚板时出现剪切锁定——这也是 PyNite 官方示例中评价它"对厚板和薄板都能给出很准确结果"的原因。

实现代码位于 Pynite/Quad3D.py，文件头部列出了 4 篇经典参考文献（Katili 的 DKMQ/DSQ/MITC4 对比研究、Bathe、Logan、Gallagher），是学习四边形板有限元公式推导的宝藏起点。

二、自由度与局部坐标系：板单元的"关节" 🦴

一个 DKMQ 四边形板单元共有16 个自由度：

4 个角节点：每个节点带 3 个弯曲自由度——横向位移w和绕两条板面轴的转角βx、βy，合计 12 个；
4 个边中点节点：每个边中点带 1 个绕边法线的转角Δβs，合计 4 个。

边中点转角的存在，正是 DKMQ 区别于普通双线性单元的关键：它让板面斜率可以表达出二次变化，从而显著改善弯曲精度。

在结构组装层面，每个角节点有 6 个自由度（3 平移 + 3 转角），4 个节点共 24 个，单元刚度矩阵为 24×24 阶。其中绕板面法线的转动（即"钻孔自由度"）在纯弯曲理论中是无约束的，PyNite 通过弱旋转弹簧（刚度取其他转动刚度的 1/1000）来保证数值稳定，这一做法在 Pynite/Quad3D.py 的类说明和 Pynite/Plate3D.py 的ke_b方法中都有体现。

单元的局部坐标系由节点顺序决定：i → j方向为局部 x 轴，法向量由叉积确定 z 轴。理解这套局部坐标是读懂所有板单元推导的前提：

PyNite有限元成员局部坐标系与截面内力方向定义

三、形状函数：双线性角点 + 不完整二次边中点

DKMQ 在自然坐标系(ξ, η) ∈ [-1, 1]上插值横向位移w，形状函数分两类（对应 Pynite/Quad3D.py 中的N_i与P_k方法）：

角点双线性函数（i = 1~4）：

$$N_i = \frac{1}{4}(1 \pm \xi)(1 \pm \eta)$$

边中点不完整二次函数（k = 5~8）：

$$P_k = \frac{1}{2}(1 - \xi^2)(1 \mp \eta),\quad P_k = \frac{1}{2}(1 \pm \xi)(1 - \eta^2)$$

注意"不完整"二字：二次项里缺少ξη交叉项。这不是偷懒，而是刻意设计——去掉交叉项后，边中点转角Δβs与角点转角在单元内部通过离散 Kirchhoff 约束保持协调：在单元内的高斯积分点上，转角与位移之间的 Kirchhoff 条件（转角 = 位移斜率）被逐点强制成立，而无需像 C¹ 连续单元那样要求跨节点斜率连续。

这一机制的完整符号推导，可以在Derivations/DMKQ Quad Element.ipynb（Jupyter Notebook，需用 Sympy 逐格运行）中逐步查看。

四、弯曲刚度矩阵的推导链条 🔗

弯曲刚度的推导是 DKMQ 公式的心脏，核心链条如下（对应 Pynite/Quad3D.py 中的各矩阵方法）：

1. 剪切耦合系数 φk（phi_k方法）：

$$\phi_k = \frac{2}{\kappa(1-\nu)}\left(\frac{t}{L_k}\right)^2,\quad \kappa = \frac{5}{6}$$

其中Lk为第 k 条边的长度，t为板厚，ν为泊松比。板越薄，φk 越小，约束越接近刚性 Kirchhoff；板越厚，约束自动"松弛"以吸收横向剪切变形。

2. 转角协调矩阵：A_Delta_inv_DKMQ方法给出−3/2 · diag(1/(1+φk))，把边中点转角 Δβ 与角点转角 β 在高斯点上协调起来。

3. 应变-位移矩阵 B（B_b方法）：

$$\mathbf{B}b = \mathbf{B}{b(\beta)} + \mathbf{B}{b(\Delta\beta)},\mathbf{A}\Delta^{-1},\mathbf{A}_u$$

其中A_u由各边长与方向余弦（dir_cos方法）构成，A_gamma、N_gamma负责转角到截面斜率的映射。

4. 雅可比矩阵与高斯积分：J方法用局部坐标构造 2×2 雅可比，将参考坐标系求导转为物理坐标，再在积分点上做BᵀDbB·det(J)加权求和，得到 12×12 弯曲刚度子矩阵。

5. 叠加膜力刚度：弯曲之上再叠加一个等参数平面应力膜单元（4 节点双线性 + 2×2 高斯积分），扩张到 24×24 后直接相加，即得单元总刚度ke = ke_b + ke_m。本构矩阵Dm（面内）与Db（弯曲，含t³/12因子）在 Pynite/Plate3D.py 中实现，还支持kx_mod/ky_mod正交各向异性刚度折减——这对模拟开裂混凝土非常实用。

五、符号约定：读懂内力结果的关键 📐

有限元公式推导的最后一环，是内力结果的提取与符号约定。PyNite 的梁弯曲符号约定如下（y 向与 z 向各一张手绘图，直观标注了荷载图、变形形状与 M(x)、V(x)、δ(x) 的关系）：

板单元的结果提取同样讲究符号：Pynite/Plate3D.py 的moment()方法通过 12 项系数矩阵C和曲率矩阵Q求得任意点(x, y)的Mx、My、Mxy；shear()由弯矩对坐标求导合成Qx、Qy；membrane()则在 4 个高斯点算应力后，用外插形状函数H平移到目标位置。四边形板则直接在等参数坐标(ξ, η)上采样，注意与矩形板的局部长度坐标区分开。

六、实战检验：与 Timoshenko 经典解对比 🏆

推导是否正确，经典解说了算。示例 Examples/Rectangular Plate Bending - Qauds.py 用 1ft×1ft 的 Quad 网格建模一面 10ft×20ft 的四周固支墙体，施加均布面压：

model.add_rectangle_mesh('MSH1', mesh_size, width, height, t, 'Concrete', 1, 1, [0, 0, 0], 'XY', element_type='Quad') model.analyze(check_statics=True)

结果与 Timoshenko《板壳理论》表 35 的解析解对比，弯矩幅值非常接近（注意两者符号约定相反）；位移略偏大，正是因为 DKMQ 考虑了横向剪切变形——厚板理论使然，而非误差。

PyNite四边形板单元剪壁结构有限元分析结果云图示例

两个实用细节：渲染的弯矩云图默认做平滑处理（对汇聚于同节点的各单元角点应力取平均），比直接用角点应力更准确；网格对象还提供max_moment/min_moment方法直接提取极值。

七、推导资料清单：跟着官方一步步推 ✅

想亲手复现整个公式推导，按下面的路径走：

📓 Derivations/DMKQ Quad Element.ipynb：DKMQ 四边形板单元的完整符号推导（需 Jupyter + Sympy，运行全部单元格即可看到输出）；
📓 Derivations/MITC4 Quad Element.ipynb：对比学习 MITC4 另一种四边形板公式；
📓 Derivations/Rectangular Plate Element.ipynb：12 项多项式矩形板的推导；
📓 Derivations/Fixed End Reactions - Linear Distributed Load.ipynb：固定端反力推导，理解等效节点力的基础；
📖 教材参考：Katili (2015)、Bathe《Finite Element Procedures》、Logan《A First Course in the Finite Element Method》、Gallagher《Finite Element Analysis Fundamentals》（完整书目见 Pynite/Quad3D.py 文件头注释）。