有限元仿真分析：从原理理解到工程判断的完整方法

有限元分析（FEA）不是软件操作，而是一种将连续数学问题转化为离散代数问题的数值方法。理解有限元的基本原理，才能在遇到异常结果时做出正确判断——是网格不够密、材料模型不对、还是边界条件设置有误。

有限元方法的核心数学原理

有限元方法的本质是变分原理——将偏微分方程的边值问题转化为泛函极值问题。以弹性力学为例，平衡方程为：

∇·σ + f = 0

对应的变分形式（虚功原理）为：

∫_Ω σ:δε dV = ∫_Ω f·δu dV + ∫_Γ t·δu dA

左式是内力虚功，右式是外力虚功。有限元的离散化是将位移场u用节点位移d和形函数N表示：

u = N·d

代入虚功原理后得到线性方程组：

K·d = F

其中K是刚度矩阵（∫BᵀDB dV），F是节点力向量。求解这个方程组就得到节点位移，再通过应变-位移关系和应力-应变关系得到应变和应力。

理解这个推导过程的意义在于：刚度矩阵K的性态（条件数、带宽、正定性）直接决定求解的精度和效率。网格畸变会导致K的条件数变差，材料参数差异大（如钢-橡胶界面）会导致K的尺度跨越多个量级，这些都会表现为求解不收敛或结果异常。

网格的本质：形函数与精度

有限元解的精度取决于形函数的阶数和网格密度。

形函数阶数

二次单元在相同网格密度下精度远高于线性单元，因为形函数可以描述弯曲的位移场。对于应力分析，线性单元在弯曲主导的区域（如梁、壳）会表现”过刚”——位移偏小，应力偏低。二次单元可以缓解这个问题，但计算量增加约4-8倍。

网格收敛性测试

有限元解的网格无关性验证是结果可信的必要条件。具体操作：

1. 用三套网格（粗、中、细）分别计算，网格尺寸比例为2:1:0.5
2. 提取关心点的物理量（位移、应力、温度等）
3. 用Richardson外推法估计精确解：

   f_exact ≈ f_fine + (f_fine – f_medium) / (r^p – 1)

   其中r是网格尺寸比（2），p是收敛阶数（线性单元p≈1-2，二次单元p≈2-3）

4. 如果|f_fine – f_medium| / f_medium < 5%，认为网格收敛

单元类型选择指南

不同物理问题和几何特征需要选择不同单元类型：

结构力学

实体单元：三维块体结构。六面体优于四面体（同精度下网格数少3-5倍）。对于复杂几何无法扫掠的情况，用二次四面体（C3D10）替代线性四面体（C3D4），精度提升显著。

壳单元：薄壁结构（厚度<1/10特征尺寸）。壳单元将三维问题降为二维，计算效率高。ANSYS的SHELL181支持厚度变化和复合材料铺层定义。

梁单元：细长结构（长度>10倍截面尺寸）。BEAM188考虑剪切变形，适用于中等长细比；纯欧拉梁理论适用于极细长结构。

什么是”锁死”（Locking）：线性四面体和低阶实体单元在模拟近乎不可压缩材料（ν>0.48，如橡胶）时会出现”体积锁死”——刚度严重偏大，位移几乎为零。解决方法：用混合单元（u-p格式，位移-压力分离插值）或高阶单元。

传热学

传热单元的选择比结构力学简单——温度是标量场，只需要一个自由度。线性四面体在传热中表现可接受（不涉及体积锁死），但界面处温度梯度仍然需要加密网格。

电磁学

电磁场的有限元需要处理矢量场（电场、磁场），用Nédélec矢量单元（棱边单元）。节点单元在处理金属-介质界面时会产生伪解，矢量单元通过将自由度放在棱边而非节点上避免了这个问题。

边界条件的物理意义

有限元中边界条件分两类：

Dirichlet（本质边界）：指定场量的值，如固定位移、已知温度。这类条件直接修改刚度矩阵的对应行，是”硬约束”。

Neumann（自然边界）：指定场量的导数，如力载荷、热流密度。这类条件出现在载荷向量F中，是”软约束”。

Robin（混合边界）：场量和导数的线性组合，如对流换热（hT + q = hT_amb）。在有限元中自然满足，不需要特殊处理。

理解圣维南原理对边界条件设置很重要：如果力的施加方式（分布力vs集中力）只在局部产生影响，远处的应力分布不受影响。这意味着工程中用集中力代替螺栓预紧力、用均匀压力代替非均匀接触压力，在远离施力区域的应力分析中是合理的。

结果验证的工程方法

有限元结果需要通过多个层次验证才能信任：

与解析解对比

简单问题有解析解（如悬臂梁端部位移δ = PL³/(3EI)），作为基准验证模型设置是否正确。如果连解析解都对不上，说明模型有根本性错误。

能量守恒检查

对于热分析，输入总热流应该等于输出总热流加储热率。在ANSYS中可以通过APDL宏计算各边界的热流积分，检查是否平衡。不平衡通常意味着边界条件遗漏或网格不匹配。

对称性检查

如果几何和载荷对称，结果也应该对称。非对称结果暗示网格不对称或求解器数值噪声。对于对称问题，可以利用对称性只建一半或四分之一模型，但对称面上的边界条件必须正确（法向位移为零、切向力为零）。

物理直觉检验

应力分布是否符合载荷传递路径？最大应力出现在预期位置？变形方向与载荷方向一致？这些”常识检查”虽然不能量化精度，但能快速发现模型设置错误——如约束反方向、载荷单位错误、材料各向异性方向设反等。

非线性分析的收敛控制

非线性分析（塑性、大变形、接触）需要迭代求解。Newton-Raphson迭代的收敛准则：

• 力收敛：|ΔF| / |F_applied| < 0.5%（默认）
• 位移收敛：|Δu| / |u| < 5%
• 能量收敛：|Δu·ΔF| / (u·F) < 5%

三个准则同时满足才算收敛。如果只有力收敛而位移不收敛，说明刚度矩阵在迭代中变化太大，需要减小载荷步长。

弧长法（Arc-Length Method）在处理结构失稳（屈曲、后屈曲）时比牛顿-拉夫逊更稳定，因为它追踪载荷-位移曲线而不是固定载荷逐步求解。ANSYS的ARCLEN命令和ABAQUS的RIKS方法是弧长法的实现。