4節点四面体要素(TET4)
理论与物理
TET4元素是什么
老师，TET4是三维分析中最基本的元素吗？
虽然是最基本的，但它是不应该在实际工作中使用的元素。我先说清楚这一点。理由稍后会详细解释，但首先我们来理解理论。
一上来就是“不要用”宣言啊…。
学习TET4的理论很重要。因为它是理解FEM基础的最佳元素。但是，如果实际工作中想要有精度的结果，就应该使用TET10。
形状函数
TET4是拥有4个节点的四面体，每个节点有3个自由度（$u, v, w$）。总计12个自由度。
TET4是拥有4个节点的四面体，每个节点有3个自由度（$u, v, w$）。总计12个自由度。
形状函数用体积坐标（barycentric coordinates）$L_1, L_2, L_3, L_4$ 表示：
位移的插值：
因为形状函数是线性（一次）的，所以位移在单元内也是线性变化的呢。
是的。位移是线性的，这意味着应变在单元内是恒定的。这是TET4的致命弱点。
恒定应变的问题
应变恒定有什么问题呢？
在应力急剧变化的部位（应力集中部位、存在弯曲应力梯度的部位），完全无法表现应力的变化。
例如，用TET4对承受弯曲的梁进行建模时：
- 理论上，截面上面受拉，下面受压，呈线性分布
- 而TET4在每个单元内应力恒定，所以只能得到阶梯状的粗略近似
- 要得到精确的弯曲应力，需要极其精细的网格
TET10（二次元）的话，应变在单元内线性变化，所以也能表现弯曲。
没错。TET10有10个节点，每条边上有中间节点，使用二次形状函数。应变在单元内线性变化，因此可以在一个单元内表现弯曲的应力分布。
TET4的刚度矩阵
TET4的单元刚度矩阵是12×12，由于B矩阵是常数矩阵，所以可以无需数值积分，以闭合形式计算：
$V_e$ 是单元体积，$[B]$ 是应变-位移矩阵（常数），$[D]$ 是本构关系矩阵。
不需要数值积分是个优点呢。
计算速度快。但是精度低，所以最终需要细化网格，导致自由度数量增加，计算速度的优势被抵消。
定量的精度比较
能给我看看TET4和TET10精度差异的具体数值吗？
悬臂梁（$L/h = 10$）的尖端挠度与理论值比较：
| 元素类型 | 网格 | DOF | 挠度误差 |
|---|---|---|---|
| TET4 | 粗糙（100个单元） | 200 | -40%（过小） |
| TET4 | 中等（1,000个单元） | 1,500 | -15% |
| TET4 | 精细（10,000个单元） | 12,000 | -5% |
| TET10 | 粗糙（100个单元） | 1,200 | -2% |
| TET10 | 中等（1,000个单元） | 9,000 | -0.3% |
在应力精度方面，差距会进一步扩大。因为TET4的应力在单元内恒定，所以应力集中的峰值无论网格细化到什么程度，都很难接近理论值。
应该使用TET4的场景
有正当理由使用TET4吗？
几乎没有。硬要说的话：
- 显式解法（explicit）的冲击分析 — TET4的质量矩阵简单，有时对显式解法的稳定时间增量有利（LS-DYNA的某些情况）
- 与流体网格的耦合 — 当CFD代码基于TET4时，有时结构侧也使用TET4以保持一致性
- 前处理器的限制 — 一些旧的自动网格生成器只能输出TET4
但这些都是“不得已的情况”，如果能使用TET10，就应该始终使用TET10。
总结
我来整理一下TET4的理论。
要点：
- 4节点，线性形状函数，单元内应变恒定 — 最基本的3D FEM元素
- 精度低 — 要表现弯曲或应力集中，需要比TET10多5~10倍的自由度
- 实际工作中不使用 — TET10作为替代品总是更优
- 最适合理论学习 — 理解FEM原理的教科书式元素
- “不要用TET4”是FEM的铁则 — 仅凭知道这一点，分析质量就能显著提高
为了理解“为什么不用”而学习，这很有趣呢。
了解TET4的局限性对于理解FEM的精度非常重要。能够解释为什么TET4不行的工程师，才是能够正确使用FEM的工程师。
TET4的线性形状函数与恒定应变特性
4节点四面体元素（TET4）拥有使用体积坐标L1~L4的线性形状函数，导致单元内应变恒定（Constant Strain Tetrahedron, CST）。1943年，库兰特（R. Courant）提出了三角形元素的概念，1956年，特纳等人将其扩展到三维四面体。作为最简单的3D实体元素，它易于实现，但由于应变恒定，在再现弯曲和应力集中时需要大量单元，这是其缺点。
各项的物理意义
- 惯性项（质量项）：$\rho \ddot{u}$，即“质量×加速度”。您有过急刹车时身体被向前甩出去的经历吗？那种“被带走的感觉”正是惯性力。物体越重越难启动，一旦启动也越难停止。地震时建筑物摇晃，也是因为地面突然移动，而建筑物的质量“跟不上”。静力分析中此项设为零，这是假设“因为缓慢施力，所以加速度可以忽略”。对于冲击载荷或振动问题，此项绝对不能省略。
- 刚度项（弹性恢复力）：$Ku$ 或 $\nabla \cdot \sigma$。拉弹簧时会感觉到“想要恢复原状的力”吧？那就是胡克定律 $F=kx$，也是刚度项的本质。那么提问——铁棒和橡皮筋，用相同的力拉，哪个伸得更长？当然是橡皮筋。这种“不易伸长性”就是杨氏模量 $E$，它决定了刚度。常见的误解：“刚度高=强度高”是不对的。刚度是“不易变形的程度”，强度是“不易破坏的程度”，是不同的概念。
- 外力项（载荷项）：体积力 $f_b$（重力等）和表面力 $f_s$（压力、接触力等）。可以这样想——桥上的卡车重量是“作用在整个内部上的力”（体积力），轮胎压路面的力是“只作用在表面上的力”（表面力）。风压、水压、螺栓的紧固力…全都是外力。这里容易犯的错误：弄错载荷方向。本想“拉伸”却变成了“压缩”——听起来像笑话，但在3D空间中坐标系发生旋转时，确实会发生。
- 阻尼项：瑞利阻尼 $C\dot{u} = (\alpha M + \beta K)\dot{u}$。弹一下吉他的弦试试。声音会一直响下去吗？不，会逐渐变小。这是因为振动能量通过空气阻力或弦的内部摩擦变成了热。汽车的减震器也是同样的原理——故意吸收振动能量来提高乘坐舒适性。如果阻尼为零会怎样？建筑物在地震后会一直摇晃下去。实际上不会这样，所以设定适当的阻尼很重要。
假设条件与适用范围
量纲分析与单位制
| 变量 | SI单位 | 注意点·换算备忘 |
|---|---|---|
| 位移 $u$ | m（米） | 输入mm时，载荷·弹性模量也要统一为MPa/N系 |
| 应力 $\sigma$ | Pa（帕斯卡）= N/m² | MPa = 10⁶ Pa。与屈服应力比较时注意单位制不一致 |
| 应变 $\varepsilon$ | 无量纲（m/m） | 注意工程应变与对数应变的区别（大变形时） |
| 弹性模量 $E$ | Pa | 钢：约210 GPa，铝：约70 GPa。注意温度依赖性 |
| 密度 $\rho$ | kg/m³ | mm系中是tonne/mm³（钢为 = 10⁻⁹ tonne/mm³） |
| 力 $F$ | N（牛顿） | mm系用N，m系也用N统一 |
数值解法与实现
TET4的实现细节
TET4的实现比其他元素简单吗？
是最简单的。因为B矩阵是常数，所以不需要数值积分。作为FEM编程的入门是最合适的。
B矩阵的计算
根据4个节点的坐标 $(x_i, y_i, z_i)$ 计算体积坐标的梯度，构成B矩阵。TET4的B矩阵是6×12的常数矩阵：
这里 $b_i, c_i, d_i$ 是从节点坐标计算出的常数。$V$ 是四面体的体积。
$V$ 怎么计算？
根据4个节点的坐标：
如果 $V > 0$，则节点顺序正确（右手系）。如果 $V < 0$，则单元翻转了。
各求解器对应的元素名称
| 求解器 | 元素名称 | 备注 |
|---|---|---|
| Nastran | CTETRA（4节点版） | PSOLID属性 |
| Abaqus | C3D4 | Standard/Explicit通用 |
| ANSYS | SOLID285 | 旧版本是SOLID45 |
| LS-DYNA | TYPE 13（四面体） | 主要用于显式冲击分析 |
| Calculix | C3D4 | 与Abaqus兼容 |
不同求解器名称不同，但理论是相同的吧？
是的。理论相同，但实现细节（如质量矩阵的集中化方法）可能略有不同。
实现代码示例
用Python（NumPy）实现TET4刚度矩阵计算的示例：
import numpy as np
def tet4_stiffness(nodes, E, nu):
"""计算TET4的刚度矩阵"""
# 节点坐标 (4x3)
x1, y1, z1 = nodes[0]
x2, y2, z2 = nodes[1]
x3, y3, z3 = nodes[2]
x4, y4, z4 = nodes[3]
# 计算体积
V = np.abs(np.linalg.det([
[x2-x1, y2-y1, z2-z1],
[x3-x1, y3-y1, z3-z1],
[x4-x1, y4-y1, z4-z1]
])) / 6.0
# 计算b, c, d系数
b1 = (y2-y4)*(z3-z4) - (y3-y4)*(z2-z4)
b2 = (y3-y4)*(z1-z4) - (y1-y4)*(z3-z4)
b3 = (y1-y4)*(z2-z4) - (y2-y4)*(z1-z4)
b4 = -b1 - b2 - b3
c1 = (z2-z4)*(x3-x4) - (z3-z4)*(x2-x4)
c2 = (z3-z4)*(x1-x4) - (z1-z4)*(x3-x4)
c3 = (z1-z4)*(x2-x4) - (z2-z4)*(x1-x4)
c4 = -c1 - c2 - c3
d1 = (x2-x4)*(y3-y4) - (x3-x4)*(y2-y4)
d2 = (x3-x4)*(y1-y4) - (x1-x4)*(y3-y4)
d3 = (x1-x4)*(y2-y4) - (x2-x4)*(y1-y4)
d4 = -d1 - d2 - d3
# 构建B矩阵 (6x12)
B = np.zeros((6, 12))
coef = 1.0 / (6.0 * V)
b = [b1, b2, b3, b4]
c = [c1, c2, c3, c4]
d = [d1, d2, d3, d4]
for i in range(4):
idx = i * 3
B[0, idx] = b[i] * coef
B[1, idx+1] = c[i] * coef
B[2, idx+2] = d[i] * coef
B[3, idx] = c[i] * coef
B[3, idx+1] = b[i] * coef
B[4, idx+1] = d[i] * coef
B[4, idx+2] = c[i] * coef
B[5, idx] = d[i] * coef
B[5, idx+2] = b[i] * coef
# 弹性矩阵D (6x6)
mu = E / (2.0 * (1.0 + nu))
lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
D = np.array([
[lam+2*mu, lam, lam, 0, 0, 0],
[lam, lam+2*mu, lam, 0, 0, 0],
[lam, lam, lam+2*mu, 0, 0, 0],
[0, 0, 0, mu, 0, 0],
[0, 0, 0, 0, mu, 0],
[0, 0, 0, 0, 0, mu]
])
# 刚度矩阵 Ke = V * B^T * D * B
Ke = V * B.T @ D @ B
return Ke确实很简单。没有数值积分循环。
这就是TET4的优点。但别忘了，精度低是致命的。
网格质量的影响
TET4对网格质量敏感吗？
非常敏感。特别是形状比（aspect ratio）差（细长）的单元，精度会急剧下降。
理想的TET4是正四面体。但实际形状中，建议：
- 形状比 < 5
- 内角 > 15度
- 雅可比比（Jacobian ratio） > 0.1
TET10对网格质量不那么敏感吗？
二次元对形状的适应性更强。但极端差的网格对任何元素都是问题。
总结
关于TET4的实现，我明白了以下几点：
相关主题
なった
詳しく
報告