等参数映射模拟器 返回
FEM模拟器

等参数映射模拟器 — 4节点四边形单元

从参考正方形 (ξ,η)∈[-1,1]² 到物理坐标 (x,y) 的映射与雅可比行列式可视化。拖动角节点,学习单元扭曲和反演(detJ≤0)的原理。

参数设置
Δx₁
Δy₁
ξ
η

默认标准四边形为4个角 (0,0)、(1,0)、(1,1)、(0,1)。仅扰动角节点1,位移量为 (Δx₁, Δy₁)。

计算结果
0.550
物理坐标 x
0.525
物理坐标 y
0.2125
雅可比行列式 det J
85.0 %
扭曲比 detJ/0.25
参考空间 (ξ,η) 与物理空间 (x,y)

左=参考正方形 [-1,1]²;右=扰动后的物理四边形。红圆=评估点 (ξ,η) 和 (x,y)。det J ≤ 0 显示\"单元反演\"警告。

理论·主要公式

4节点四边形(Q4)单元的形状函数(双线性)。角节点编号对应于 1=(-1,-1)、2=(+1,-1)、3=(+1,+1)、4=(-1,+1):

$$N_1=\tfrac{1}{4}(1-\xi)(1-\eta),\;N_2=\tfrac{1}{4}(1+\xi)(1-\eta),\;N_3=\tfrac{1}{4}(1+\xi)(1+\eta),\;N_4=\tfrac{1}{4}(1-\xi)(1+\eta)$$

到物理坐标的映射(坐标和未知量使用相同形状函数=等参数):

$$x=\sum_{i=1}^{4} N_i(\xi,\eta)\,x_i,\qquad y=\sum_{i=1}^{4} N_i(\xi,\eta)\,y_i$$

雅可比矩阵和行列式(单元扭曲指示器):

$$J=\begin{bmatrix}\partial x/\partial\xi & \partial y/\partial\xi\\ \partial x/\partial\eta & \partial y/\partial\eta\end{bmatrix},\qquad dx\,dy=|\det J|\,d\xi\,d\eta$$

标准正方形(Δx=Δy=0)时 det J = 0.25。det J ≤ 0 时单元反演,刚度矩阵数值积分失效。

等参数映射模拟器介绍

🙋
FEM教科书里提到\"等参数单元\",\"等\"是什么意思呢?
🎓
简单地说,就是\"用来补间形状的函数\"和\"用来补间位移的函数\"是同一种。我们先准备一个参考正方形 $(\xi,\eta)\in[-1,1]^2$,用形状函数 $N_i$ 把它映射到实际的四边形。同时节点的位移也用同样的 $N_i$ 在单元内部补间。这样一个参考单元就可以生成任意形状的单元,是很有用的。看上面的模拟器,试试拖动角节点1,你会看到右边物理四边形变形的过程。
🙋
\"雅可比行列式\"卡片显示的是什么?
🎓
雅可比行列式表示\"参考空间的小面积在物理空间中扩大了多少倍\"。公式是 $dx\,dy = |\det J|\,d\xi\,d\eta$。计算刚度矩阵时的面积积分必然涉及这个值。理想的标准正方形(Δx=Δy=0)时 det J = 0.25。当你拖动角节点1时,这个值会变化。\"扭曲比\"卡片告诉你这个值是标准值 0.25 的多少倍。
🙋
把Δx₁或Δy₁拖到负方向时,出现了\"单元反演\"的红字。这是什么意思?
🎓
那就是 det J ≤ 0 的状态。物理空间的四边形\"翻转\"了。比如角节点1拖到角节点3的另一边,四条边会交叉形成蝴蝶结形。这样形状函数的映射就不再是一一对应的,刚度矩阵的积分会出错。在实际的大变形分析中,\"求解器崩溃\"最常见的原因就是这个,通常用重新网格化或ALE来处理。
🙋
背景里有网格线,那是什么?
🎓
那是参考空间中把 $\xi$ 和 $\eta$ 每 0.05 步进一次,然后用形状函数映射到物理空间的曲线。标准正方形时是直线网格,但拖动角节点1后,网格线会弯曲。这在视觉上表现了雅可比矩阵在每个位置的值都不一样。FEM在数值积分(高斯积分)时就是这样考虑这个效果的。

常见问题

这三种分类基于坐标补间函数的次数和位移补间函数的次数关系。两者次数相同时称为\"等参数\"(标准的FEM单元就是这样),坐标补间次数较低称\"亚参数\",较高称\"超参数\"。等参数在数值计算上最稳定,刚体平移和旋转的表示也很自然,所以商业FEM求解器大多采用等参数单元。
Q4单元通常在参考坐标 $(\pm 1/\sqrt{3},\pm 1/\sqrt{3})$ 的4个点进行2×2高斯积分(完全积分)。这可以准确计算弯曲,但容易产生剪切自锁现象。为了对付剪切自锁,有\"缩减积分\"(仅剪切项用1×1积分)、\"选择性缩减\"、\"$\\bar{B}$方法\"等。当单元扭曲较大时,积分点的物理坐标映射精度下降,所以网格质量管理很重要。
商业FEM求解器通常用单元内最大det J和最小det J的比值来表示扭曲比。一般来说1~10为优秀,10~40为警告,40以上或负值(单元反演)为出错停止。另外还要检查边长比(最长边/最短边),线性单元通常要求5~10以内,二次单元大约20。扭曲单元在应力集中区域精度下降明显,所以应力集中部位要特别注意。
可以。8节点六面体单元(H8)就是用完全相同的框架定义的。参考立方体 $(\xi,\eta,\zeta)\in[-1,1]^3$ 映射到物理空间,用三线性形状函数 $N_i = \\tfrac{1}{8}(1\\pm\\xi)(1\\pm\\eta)(1\\pm\\zeta)$。雅可比矩阵变成3×3,行列式表示体积比。三角形和四面体用面积坐标或体积坐标,思路是一致的。

实际应用

商用FEM求解器的基础:Abaqus、ANSYS、Nastran、LS-DYNA等主流商业FEM求解器都采用等参数单元。从线性4节点四边形(Q4、CPS4等)到二次8节点(Q8、CPS8等),再到三维六面体和四面体,所有单元都建立在这个理论之上。能够用网格覆盖任意形状结构的灵活性是现代CAE的基础。

大变形分析和成形模拟:金属冲压、橡胶密封件压缩、撞击分析等单元发生大变形的分析中,逐步检查单元扭曲比和 det J 正值性是必须的。当 det J 变得太小时采用\"自适应重网格化\",或者使用ALE(任意拉格朗日欧拉)方法分离参考网格和物理网格。

等几何分析(IGA)的发展:等参数的思想被发展成为\"在CAD的NURBS基函数上直接进行单元计算\"的等几何分析(Isogeometric Analysis, IGA)。这样可以消除CAD和CAE的网格转换损失,在薄板和壳体分析中获得很高的精度。从Q4单元到IGA,\"用同样的函数补间形状和未知量\"的哲学是一脉相承的。

教育和FEM入门的标准题材:形状函数、雅可比矩阵、数值积分是FEM的核心概念,而能以最直观的方式理解这些概念的最小单元就是Q4等参数单元。大多数教科书首先介绍这个单元,编程练习也把它作为第一个任务。这个模拟器就是为了可视化这个教学内容。

常见误解和注意事项

最常见的误解是\"四边形单元必须接近正方形\"。实际上梯形、平行四边形、任意凸四边形都可以,只要 det J 为正且足够大就行。反而要网格化复杂实际形状的话,必然要接受一定的扭曲。如果在模拟器中设定 Δx₁=0.3、Δy₁=0.2,你会看到接近梯形的四边形,det J 的扭曲比仍保持在 60~80%,是完全实用的。

另一个常见误解是\"det J 越大越好\"。det J 只是参考空间和物理空间的面积比,会随着单元大小按比例增大。重要的不是 det J 的绝对值,而是单元内各积分点间 det J 的比值(即Jacobian Ratio),以及最小det J 保持为正。本模拟器的\"扭曲比\"卡片相对于标准正方形 0.25 的百分比,100%附近是理想状态、40%以下是警告、0%以下是反演,可以这样理解。

最后要注意\"形状函数总和 ΣN_i = 1 的条件容易被忽视\"。ΣN_i = 1 是正确表示刚体平移的必要条件,如果这个条件破坏了,单元在刚体运动时也会产生内部应力。可以验证 $N_1+N_2+N_3+N_4 = \\tfrac{1}{4}[(1-\\xi)(1-\\eta)+(1+\\xi)(1-\\eta)+(1+\\xi)(1+\\eta)+(1-\\xi)(1+\\eta)] = 1$ 对任意 $(\\xi,\\eta)$ 成立。新开发单元类型时,必须验证这个条件和\"节点处 $N_i$ 满足克罗内克δ\"两个条件。

使用指南

  1. 在左侧面板用滑块调节4节点四边形单元的节点坐标。节点1~4分别对应物理空间的各个顶点位置。
  2. 用滑块指定参考空间(ξ-η坐标系)的评估点。ξ、η的范围都是-1~1,代表从中心(0,0)的相对位置。
  3. 右侧面板实时显示雅可比行列式det J、物理坐标(x,y)、扭曲比(detJ/0.25)的计算结果。det J接近0.25的正方形是最优单元。

具体计算例

5毫米铝板穿孔解析中,400mm×300mm区域用4节点四边形网格化,节点坐标设为(0,0)、(100,0)、(100,80)、(0,80)毫米。在参考单元中心(ξ=0,η=0)评估时det J≈2000 mm²。将同一区域梯形化,节点改为(0,0)、(105,0)、(95,85)、(5,80)毫米时,det J≈1850 mm²,扭曲比0.925,可以定量评估单元质量的劣化。

实务中的注意点

  1. 雅可比行列式出现负值表示节点顺序颠倒。必须坚持反时针的节点编号规则。
  2. 自动网格生成后,det J低于0.15的单元在应力集中计算中误差会增大,应进行局部重网格化。
  3. 非结构网格中相邻单元det J差超过50%时,表示应力梯度陡峭区域。应设计单元尺寸的阶梯式过渡。