Go back

CAE 与有限元分析基础(课程笔记)

Published:  at  02:00 PM Updated:  at  11:30 PM
阅读时间：33 分钟
编辑此页面

这是一门 CAE / 有限元法 课程的持续更新笔记。目标是把”CAE 基础”讲全、讲透、讲懂:既有公式推导,也有物理直觉。最新更新在顶部日期。

目录

〇、写在前面:CAE 是什么、为什么要学有限元

CAE(Computer-Aided Engineering,计算机辅助工程)就是用计算机仿真来分析工程问题,不用每个设计都造实物去试。它在产品研发流程里的位置:

CAD(设计)  ⟶  CAE(仿真分析)  ⟶  CAM(制造)\text{CAD(设计)} \;\longrightarrow\; \boxed{\text{CAE(仿真分析)}} \;\longrightarrow\; \text{CAM(制造)}

CAE 覆盖很多物理场:

类型解决的问题代表场景
FEA(结构)强度、刚度、变形、振动零件会不会断、变形多大
CFD(流体)流场、压力、阻力汽车风阻、机翼升力
热分析温度场、热应力芯片散热
电磁电磁场电机、天线
多物理场耦合(如热-力、流-固)涡轮叶片热应力

常用软件:ANSYS、Abaqus、COMSOL、HyperMesh(前处理网格)、Nastran、SolidWorks Simulation 等。

那为什么要学”有限元法(FEM)”? 因为几乎所有结构 CAE 软件的底层引擎都是有限元。软件是黑盒,但:

  • 懂原理 → 知道怎么建对模型(网格怎么画、约束怎么加)
  • 懂原理 → 能判断结果可不可信(是不是数值假象、收敛了没)
  • 懂原理 → 出问题时知道往哪查

这门课的主线就是:把 CAE 软件黑盒拆开,看清 Ku=FKu=F 这个方程是怎么从物理定律一步步变出来的。

一、有限元法(FEM)的核心思想

1.1 一句话思想

有限元法是一种数值方法:把一个连续的整体,切成有限个小块(单元),用每个小块上的简单近似,拼出整体解答。

关键词是离散化(discretization):连续体 → 有限个”单元(element)“通过”节点(node)“连接。

1.2 为什么这么做?

力学问题本质是偏微分方程 + 边界条件(比如弹性力学的平衡方程)。解析法要求方程在每一点都成立,复杂几何根本解不了。

有限元的思路是”降低要求、分而治之”:

  • 不要求每一点精确,只要整体上(加权)满足就行(弱形式)
  • 把结构切成单元,每个单元内用一个**简单函数(形函数)**近似真实位移
  • 单元越小越密,近似越准(收敛性)

1.3 求解的三大步

  1. 离散化:结构 → 网格(节点 + 单元)
  2. 单元分析:每个单元内,用形函数近似位移场 → 导出单元刚度矩阵
  3. 组装求解:所有单元刚度拼成全局刚度矩阵,加边界条件,解线性方程组

所有线性弹性问题最终都归约成同一个代数方程:

K u=FK\,u = F
符号含义直觉
KK全局刚度矩阵结构的”硬度”,多大位移需要多大力
uu节点位移向量待求 —— 每个节点移动了多少
FF节点载荷向量外力在每个节点上的等效

这就像弹簧 F=kuF=ku:力 = 刚度 × 位移。有限元只是把”一根弹簧”推广成了”成千上万根弹簧组成的网”,写成矩阵形式 Ku=FKu=F。

1.4 和其他方法的对比

方法思路优缺点
解析法求微分方程精确解准,但只对简单几何可行
有限差分法(FDM)在网格点用差商代替导数简单,但对复杂边界、非均匀网格不友好
有限元法(FEM)分片近似 + 变分/加权残值几何适应性强、理论基础牢,工程主流
边界元(BEM)只在边界离散单元少,但矩阵满、推导难

有限元赢在几何适应性(任意形状都能切网格)和理论完备(变分原理保证收敛),所以工程上几乎一统天下。

二、弹性力学基础:应力、应变、本构

有限元分析固体,底层是弹性力学。三个核心物理量:应力(力怎么分布)、应变(变形多大)、本构关系(力和变形的关系)。

2.1 应力 σ\sigma

应力是”内力在截面上的分布强度”,单位是 Pa(N/m²)。

想象一个受拉杆,在内部假想切一刀,截面上单位面积的拉力就是应力。三维下,一点的应力状态要用张量完整描述(因为不同方向截面上的应力不同)。柯西应力张量:

σ=[σxxτxyτxzτyxσyyτyzτzxτzyσzz]\sigma = \begin{bmatrix} \sigma_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \sigma_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \sigma_{zz} \end{bmatrix}
  • 对角线 σxx,σyy,σzz\sigma_{xx},\sigma_{yy},\sigma_{zz}:正应力(垂直于截面的拉/压)
  • 非对角 τxy\tau_{xy} 等:剪应力(沿截面的切向力)

剪应力互等定理:由力矩平衡可证 τxy=τyx\tau_{xy}=\tau_{yx}, τyz=τzy\tau_{yz}=\tau_{zy}, τzx=τxz\tau_{zx}=\tau_{xz},所以张量对称,9 个分量只有 6 个独立。

主应力:总可以旋转坐标系,让某斜面上剪应力为零,此时三个正应力叫主应力 σ1≥σ2≥σ3\sigma_1\geq\sigma_2\geq\sigma_3。它们是设计的关键依据(比如第四强度理论用它们判屈服)。

2.2 应变 ε\varepsilon

应变是”变形的度量”,无量纲。

最直观:一根长 LL 的杆拉长了 ΔL\Delta L,工程正应变就是:

ε=ΔLL\varepsilon = \frac{\Delta L}{L}

三维下,小变形假设(变形远小于尺寸),应变-位移关系(几何方程)为:

εij=12(∂ui∂xj+∂uj∂xi)\varepsilon_{ij} = \frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)

写成工程形式(正应变 + 剪应变):

εxx=∂u∂x,γxy=∂u∂y+∂v∂x,  …\varepsilon_{xx}=\frac{\partial u}{\partial x},\quad \gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x},\;\dots

(注意工程剪应变 γxy=2εxy\gamma_{xy}=2\varepsilon_{xy},差个系数 2,初学容易混)

直觉:应变是位移梯度的对称部分。εxx=∂u/∂x\varepsilon_{xx}=\partial u/\partial x 就是 x 方向单位长度的伸长率。它把”位移场”(运动学)和”变形”(几何)联系起来。

2.3 本构关系:广义胡克定律

应力和应变怎么联系起来?靠材料属性。

一维最简单:胡克定律 σ=E ε\sigma = E\,\varepsilon,EE 是弹性模量(杨氏模量)。

三维线弹性各向同性材料,广义胡克定律:

σ=D ε\sigma = D\,\varepsilon

DD 是弹性矩阵(由 EE 和泊松比 ν\nu 决定)。这是有限元里把”应变”翻译成”应力”的关键。

平面应力情形(薄板受面内力,常用):

D=E1−ν2[1ν0ν10001−ν2]D = \frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \dfrac{1-\nu}{2} \end{bmatrix}

平面应变情形(厚结构沿 z 方向约束):

D=E(1−ν)(1+ν)(1−2ν)[1ν1−ν0ν1−ν10001−2ν2(1−ν)]D = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)}\begin{bmatrix} 1 & \dfrac{\nu}{1-\nu} & 0 \\ \dfrac{\nu}{1-\nu} & 1 & 0 \\ 0 & 0 & \dfrac{1-2\nu}{2(1-\nu)} \end{bmatrix}

三个材料参数的关系:

G=E2(1+ν)G = \frac{E}{2(1+\nu)}
  • EE:弹性模量,越大越”硬”
  • ν\nu:泊松比,横向收缩比(拉长时变细),金属约 0.3,橡胶近 0.5,不可压缩材料 0.5
  • GG:剪切模量

平面应力 vs 平面应变别搞混:薄板(厚度小、两面自由)用平面应力;很厚的结构(如水坝、长挡土墙,沿长度方向不变形)用平面应变。两者 D 矩阵不同,结果不同。

2.4 弹性力学三大方程 + 边界条件

弹性力学的全部理论围绕三组方程:

方程形式含义
平衡方程∂σij∂xj+bi=0\dfrac{\partial \sigma_{ij}}{\partial x_j} + b_i = 0内应力梯度与体力 bb 平衡
几何方程εij=12(ui,j+uj,i)\varepsilon_{ij}=\dfrac12(u_{i,j}+u_{j,i})应变-位移关系
本构方程σ=D ε\sigma = D\,\varepsilon应力-应变关系

加上边界条件:

  • 力边界 Γt\Gamma_t:给定的面力 σn=tˉ\sigma n = \bar t(某面受力)
  • 位移边界 Γu\Gamma_u:给定的位移 u=uˉu=\bar u(某面固定/约束)

三大方程 + 边界条件 = 强形式(strong form),每一点都要满足。问题是:复杂几何下强形式解不出来。这就引出弱形式。

三、从强形式到弱形式:FEM 的数学根基

这一章是有限元”为什么能算”的核心。新手常卡在这,我尽量讲直觉。

3.1 强形式的”苛刻”

强形式要求微分方程在区域内每一点都成立。但真实结构形状复杂、载荷复杂,精确满足几乎不可能。我们要”放宽要求”。

3.2 弱形式(加权残值法)

思路:找一个近似解 u~\tilde u,它不满足方程时会有残值 R=L(u~)−f≠0R = L(\tilde u) - f \neq 0。我们不让残值处处为零(太难),只要求残值在”加权平均”意义上为零:

∫ΩW R  dΩ=0\int_\Omega W\,R\;\mathrm d\Omega = 0

WW 是权函数(test function)。不同选 WW 的方法不同:Galerkin 法最常用 —— 权函数取成和试探函数(形函数)同族。

直觉:强形式像”每个学生都必须考满分”(做不到);弱形式像”全班平均分达到要求就行”(可近似)。分片近似的有限元,就是在每个单元上用简单函数凑这个”平均分”。

通过分部积分,平衡方程的弱形式可化成(对弹性问题):

∫ΩδεTσ dΩ=∫ΩδuTb dΩ+∫ΓtδuTtˉ dΓ\int_\Omega \delta\varepsilon^T \sigma\,\mathrm d\Omega = \int_\Omega \delta u^T b\,\mathrm d\Omega + \int_{\Gamma_t} \delta u^T \bar t\,\mathrm d\Gamma

左边 = 内力虚功,右边 = 外力虚功。这就是虚功方程。

3.3 虚功原理

引入虚位移 δu\delta u(想象的、满足约束的、无穷小位移),虚功原理表述为:

  内力虚功=外力虚功  \boxed{\;\text{内力虚功} = \text{外力虚功}\;}

即弹性体平衡的充要条件是:对任意虚位移,内力(应力)做的虚功等于外力(体力 + 面力)做的虚功。

虚功原理是弱形式的物理化身。它把”微分方程每点成立”换成”能量平衡”,后者能离散成代数方程。

3.4 最小势能原理

等价的另一表述:真实位移使系统总势能取驻值(最小)。

总势能:

Π(u)=12∫ΩεTD ε dΩ⏟应变能 U−(∫ΩuTb dΩ+∫ΓtuTtˉ dΓ)⏟外力功 W\Pi(u) = \underbrace{\tfrac12\int_\Omega \varepsilon^T D\,\varepsilon\,\mathrm d\Omega}_{\text{应变能 } U} - \underbrace{\left(\int_\Omega u^T b\,\mathrm d\Omega + \int_{\Gamma_t} u^T\bar t\,\mathrm d\Gamma\right)}_{\text{外力功 } W}

真实解满足 δΠ=0\delta\Pi=0。Ritz 法就是用形函数近似 uu,代入 Π\Pi,对节点未知量求驻值,正好得到 Ku=FKu=F。

三条路径殊途同归:加权残值(Galerkin)/ 虚功原理 / 最小势能原理,对线弹性都导出同一个 Ku=FKu=F。这也是为什么有限元有坚实的数学基础。

四、形函数(插值函数)

4.1 形函数的作用

有限元的”魔法”全在形函数:单元内任意一点的位移,用节点位移插值出来。

u(x)=∑i=1nNi(x) ui=N ueu(x) = \sum_{i=1}^{n} N_i(x)\,u_i = N\,u_e
  • uiu_i:节点 ii 的位移(未知量)
  • Ni(x)N_i(x):形函数(shape function),只和位置有关
  • ueu_e:单元节点位移向量

形函数的两条基本性质:

  1. Kronecker 性质:NiN_i 在节点 ii 处 = 1,在其他节点处 = 0
  2. ** partition of unity**:所有形函数之和 = 1(保证常位移能精确表示)

4.2 一维线性形函数(杆单元)

最简单的 2 节点杆单元,节点 1(x1x_1)、节点 2(x2x_2),长 L=x2−x1L=x_2-x_1:

N1(x)=x2−xL,N2(x)=x−x1LN_1(x)=\frac{x_2-x}{L},\qquad N_2(x)=\frac{x-x_1}{L}

验证:x=x1x=x_1 时 N1=1,N2=0N_1=1,N_2=0;x=x2x=x_2 时 N1=0,N2=1N_1=0,N_2=1;N1+N2=1N_1+N_2=1。✓

位移线性插值 u(x)=N1u1+N2u2u(x)=N_1 u_1 + N_2 u_2。这就是为什么杆单元内位移是线性分布、应变是常数。

4.3 自然坐标与等参元

为了统一处理任意形状单元,引入自然坐标 ξ∈[−1,1]\xi\in[-1,1](局部坐标),把物理坐标 xx 也用形函数插值:

x(ξ)=∑Ni(ξ) xix(\xi) = \sum N_i(\xi)\,x_i

如果位移和坐标用同一套形函数插值,就叫等参元(isoparametric element)。这是现代有限元的主流,让任意四边形/六面体/曲边单元都能用高斯积分统一计算。

五、杆单元:最简单的有限元

杆(bar/truss element)是入门最佳对象:一维、2 节点、只受轴力。完整推导一遍,你就理解了所有单元的推导套路。

5.1 单元模型

2 节点杆,长 LL,截面积 AA,弹性模量 EE。节点位移 u1,u2u_1,u_2(沿杆轴线),节点力 f1,f2f_1,f_2。

5.2 从位移到刚度(三步推导)

第一步:位移插值(用形函数)

u(x)=N1u1+N2u2,N=[x2−xLx−x1L]u(x)=N_1 u_1 + N_2 u_2,\qquad N=\begin{bmatrix}\dfrac{x_2-x}{L} & \dfrac{x-x_1}{L}\end{bmatrix}

第二步:应变(应变-位移关系,几何方程)

ε=dudx=dNdx ue=B ue\varepsilon = \frac{\mathrm du}{\mathrm dx} = \frac{\mathrm dN}{\mathrm dx}\,u_e = B\,u_e

其中应变-位移矩阵:

B=dNdx=[−1L1L]B = \frac{\mathrm dN}{\mathrm dx} = \begin{bmatrix} -\dfrac{1}{L} & \dfrac{1}{L} \end{bmatrix}

(形函数是线性的,求导是常数 → 杆单元内应变是常数,符合直觉)

第三步:单元刚度矩阵(由虚功原理 / 应变能)

应变能 U=12∫εTσ dVU=\tfrac12\int \varepsilon^T\sigma\,\mathrm dV,代入 ε=Bue\varepsilon=B u_e、σ=Dε=Eε\sigma=D\varepsilon=E\varepsilon:

U=12 ueT∫0LBTEB A dx⏟ke ueU = \tfrac12\,u_e^T \underbrace{\int_0^L B^T E B\,A\,\mathrm dx}_{k_e}\,u_e

得单元刚度矩阵:

  ke=EAL[1−1−11]  \boxed{\;k_e = \frac{EA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\;}

5.3 物理意义

  • 对角元 EAL\frac{EA}{L}:在节点 1 加单位位移,需要 EAL\frac{EA}{L} 的力 → 这就是”杆的刚度”k=EA/Lk=EA/L
  • 非对角 −EAL-\frac{EA}{L}:节点 1 加力,节点 2 会受反向力(牛顿第三定律)
  • 行和 = 0 → 矩阵奇异(没加约束时刚体位移不定),这是后续必须加边界条件的原因

这根杆,其实就是一根”刚度为 EA/LEA/L 的弹簧”,有限元只是把它写成矩阵形式。

六、桁架结构:坐标变换与组装

6.1 桁架特点

桁架(truss) 由杆铰接而成,杆只受轴向力(拉/压),节点是理想铰(不传递弯矩)。桥梁、屋架、起重机臂都是桁架。

分析桁架要解决:杆是斜的(局部 1D),结构是 2D/3D 的,怎么把单根杆的刚度”装”进整体?

6.2 局部坐标 vs 整体坐标

单根杆的 kek_e(上节)是在杆的局部坐标(沿杆轴)推导的。但结构里杆是斜的,要转换到整体坐标(全局 X-Y)。

设杆与整体 x 轴夹角 θ\theta,方向余弦 c=cos⁡θc=\cos\theta, s=sin⁡θs=\sin\theta。2D 杆每节点 2 个自由度(整体坐标下),变换矩阵:

T=[cs0000cs]T = \begin{bmatrix} c & s & 0 & 0 \\ 0 & 0 & c & s \end{bmatrix}

整体坐标下的单元刚度:

  Ke=TT ke T  \boxed{\;K_e = T^T\,k_e\,T\;}

展开后是 4×44\times4 矩阵(2 节点 × 2 自由度):

Ke=EAL[c2cs−c2−cscss2−cs−s2−c2−csc2cs−cs−s2css2]K_e = \frac{EA}{L}\begin{bmatrix} c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2 \end{bmatrix}

直觉:斜杆在整体坐标里”投影”到各方向,刚度分量按 cos⁡2θ\cos^2\theta、cos⁡θsin⁡θ\cos\theta\sin\theta、sin⁡2θ\sin^2\theta 分配。

6.3 直接刚度法组装

把所有单元的 KeK_e 按自由度编号叠加到全局 KK:

K=∑eKe(扩展到全局尺寸)K = \sum_e K_e^{(\text{扩展到全局尺寸})}

每个单元的 4 个自由度对应全局的某 4 个自由度,把这 4×4 的小块”放”到全局 KK 的对应位置,重叠的相加。

这就是”积零为整”。全局 KK 的大小 = 总自由度数 × 总自由度数。

七、刚度矩阵:性质、组装与边界条件

7.1 刚度矩阵的性质

单元刚度矩阵和全局刚度矩阵都有这些性质:

  1. 对称:Kij=KjiK_{ij}=K_{ji}(由应变能对称性 / Maxwell-Betti 互等定理)
  2. 正定(施加足够约束后):能量 12uTKu≥0\tfrac12 u^T K u \geq 0
  3. 奇异(无约束时):存在刚体位移模式,∣K∣=0|K|=0,不能直接求逆
  4. 稀疏 + 带状:只有相邻节点自由度耦合,大矩阵大部分是零(工程软件靠这个省内存)

7.2 为什么必须加边界条件?

不加约束的 KK 是奇异的(刚体位移不确定),det⁡K=0\det K=0,K−1K^{-1} 不存在,Ku=FKu=F 解不出。必须固定结构的刚体运动(至少约束掉所有刚体平动和转动),才能求解。

7.3 边界条件施加方法

已知某自由度位移 ui=uˉiu_i=\bar u_i(常是固定 uˉi=0\bar u_i=0),两种主流做法:

① 置 1 法(降阶法):把 uiu_i 对应的行列从 K,FK,F 中删掉,降阶求解。概念清晰,但编程要重排。

② 置大数法(惩罚法):把 KiiK_{ii} 改成极大值 β\beta(如 102010^{20}),FiF_i 改成 βuˉi\beta\bar u_i。求解时 uiu_i 被强制 ≈ uˉi\bar u_i。编程简单,商业软件常用。

数学上(惩罚法):

Kii→β,Fi→β uˉiK_{ii}\to \beta,\qquad F_i \to \beta\,\bar u_i

7.4 求解

加完边界条件,Ku=FKu=F 变成可解的对称正定方程组,解出节点位移 uu:

u=K−1Fu = K^{-1}F

(实际软件用 LU、Cholesky 等高效解法,不求显式逆,因为 KK 稀疏且大)

求出 uu 后,回代算每个单元的应变 ε=Bue\varepsilon=B u_e 和应力 σ=Dε\sigma=D\varepsilon,以及约束处的支座反力 R=Ku−FR=Ku-F。

八、完整求解流程(带手算例子)

把前面所有步骤串起来。以一根最简单的两段杆为例,手算一遍。

问题:两根同材料同截面杆串联,EAEA、长度都为 LL,左端固定,右端加拉力 PP。求节点位移。

  固定          节点2           节点3(自由)
   ●──────●──────●  ← P
   1      单元①   2     单元②    3
          长L           长L

步骤 1:节点/单元编号

节点 1(固定）、2（中间）、3（自由）。单元①= 节点 1-2,单元②= 节点 2-3。

步骤 2:单元刚度(都在同一轴向,无需坐标变换)

k(1)=k(2)=EAL[1−1−11]k^{(1)}=k^{(2)}=\frac{EA}{L}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}

步骤 3:组装全局 KK(3 个自由度 u1,u2,u3u_1,u_2,u_3)

把每个 k(e)k^{(e)} 放到对应位置叠加:

K=EAL[1−10−12−10−11]K=\frac{EA}{L}\begin{bmatrix} 1 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{bmatrix}

(中间节点 2 被两个单元共享,所以 K22=1+1=2K_{22}=1+1=2)

步骤 4:加边界条件

节点 1 固定 u1=0u_1=0。用置 1 法删掉第 1 行列:

EAL[2−1−11][u2u3]=[0P]\frac{EA}{L}\begin{bmatrix}2&-1\\-1&1\end{bmatrix}\begin{bmatrix}u_2\\u_3\end{bmatrix}=\begin{bmatrix}0\\P\end{bmatrix}

步骤 5:求解

[u2u3]=LEA[2−1−11]−1[0P]=LEA[1112][0P]=PLEA[12]\begin{bmatrix}u_2\\u_3\end{bmatrix}=\frac{L}{EA}\begin{bmatrix}2&-1\\-1&1\end{bmatrix}^{-1}\begin{bmatrix}0\\P\end{bmatrix}=\frac{L}{EA}\begin{bmatrix}1&1\\1&2\end{bmatrix}\begin{bmatrix}0\\P\end{bmatrix}=\frac{PL}{EA}\begin{bmatrix}1\\2\end{bmatrix}
  • u2=PL/EAu_2=PL/EA(中间节点位移)
  • u3=2PL/EAu_3=2PL/EA(右端位移,是两段伸长之和,符合直觉 ✓)

步骤 6:应力 / 反力

单元②应变 ε(2)=(u3−u2)/L=(2PL/EA−PL/EA)/L=P/EA\varepsilon^{(2)}=(u_3-u_2)/L=(2PL/EA-PL/EA)/L=P/EA,应力 σ=Eε=P/A\sigma=E\varepsilon=P/A。✓

支座反力 R1=−(K u)1=−PR_1 = -(K\,u)_1 = -P(和所加外力 PP 平衡,验证了整体平衡 ✓)。

这就是有限元求解一个问题的完整闭环:建模 → 单元刚度 → 组装 → 边界条件 → 解位移 → 算应力/反力。再复杂的结构、再多的单元,流程完全一样,只是矩阵更大,交给计算机。

九、网格、收敛与误差

9.1 收敛性

有限元结果是近似解,网格越密越接近真解。两种收敛方式:

  • h-收敛:单元尺寸 hh 减小(加密网格)
  • p-收敛:提高形函数阶数 pp(线性→二次→高次)

工程上常做 h-细化:加密关键区域(应力集中处),看结果是否稳定(收敛性研究)。

9.2 网格质量

网格差会引入误差,关注:

  • 长宽比(aspect ratio):单元别太扁
  • 雅可比 / 畸变:单元别严重扭曲(否则 BB 矩阵病态)
  • 过渡:疏密过渡要平缓
  • 应力集中区:加密,最好用规则单元

9.3 误差来源与验证

误差来自:

  1. 离散误差(单元近似)—— 加密网格可减
  2. 数值积分误差(高斯积分点不足)
  3. 模型简化(平面化、边界条件近似)
  4. 舍入误差(大矩阵计算)

验证手段:

  • 反力之和 = 外力之和(整体平衡)
  • 与解析解/手册解对比(简单情形)
  • 网格收敛性研究(加密一倍,结果变化 < 几 % 算收敛)
  • 对称结构结果应对称

十、CAE 工程实践要点

这些是”用软件不踩坑”的实战经验,新手最容易忽略。

  1. 单位制必须自洽:全模型统一(如 mm-N-s-MPa,或 m-N-s-Pa)。最常见错误是混单位(如长度用 mm、弹性模量填 Pa),结果差 10610^6 倍。

  2. 圣维南原理:远离载荷施加处的应力分布,与载荷具体施加方式无关。→ 应力集中只看局部,远处分布很快”平滑”。

  3. 边界条件要”够”:不能少约束(机构/可动,奇异),也别过约束(多约束会引入虚假内力)。对称结构可用对称边界减半模型。

  4. 载荷简化:集中力其实分布在小面积上,别在单节点加极大集中力(会数值奇异),用耦合/分布更稳。

  5. 结果看什么:

    • 先看变形(整体对不对、量级合不合理)
    • 再看应力(注意应力集中、单元边界应力跳跃 = 网格不够)
    • 最后反力平衡验证

  6. 奇异点:尖角、点载荷、点约束处应力会”发散”(理论上无穷),别迷信这些点的数值,要看稍远处的值或圆角处理。

十一、常用软件与学习路径

常用软件

软件特点
ANSYS工业界最广,多物理场齐全
Abaqus非线性(材料/几何/接触)强,学术界常用
COMSOL多物理场耦合友好,GUI 直观
HyperMesh前处理/网格划分强,常配其他求解器
Nastran航空航天、振动分析经典
SolidWorks Simulation集成 CAD,快速初评

学习路径建议

  1. 先懂原理(这份笔记的主线):Ku=FKu=F 怎么来的,为什么加约束,怎么验证结果
  2. 手算几个小例子(杆/桁架):建立对”组装-求解”的肌肉记忆
  3. 上手一个软件(ANSYS 或 Abaqus):照着例题做,理解每一步对应原理
  4. 做收敛性研究:同一个问题不同网格,看结果变化
  5. 进阶:非线性、动力学、传热、流体……

十二、二维单元(平面问题)

前面杆/桁架是一维。真实结构多是二维(平板)或三维。这一章讲平面问题(平面应力/平面应变)的两类基本单元:三角形单元和四边形单元。

回顾:平面问题每节点有 2 个自由度 (u,v)(u,v),应变向量 ε=[εx,εy,γxy]T\varepsilon=[\varepsilon_x,\varepsilon_y,\gamma_{xy}]^T,弹性矩阵 DD 见 §2.3。

12.1 三节点三角形单元(CST,Constant Strain Triangle)

最简单的二维单元。3 个节点 (xi,yi)(x_i,y_i),每节点 2 自由度,单元共 6 自由度。

① 面积坐标(自然坐标)

三角形内任一点 PP,把三角形分成 3 个小三角形(面积 A1,A2,A3A_1,A_2,A_3),定义面积坐标:

Li=AiA,L1+L2+L3=1L_i = \frac{A_i}{A},\quad L_1+L_2+L_3=1

其中 AA 是三角形总面积。面积坐标和直角坐标是线性关系,可互相转换。

② 形函数(线性)

对 3 节点线性单元,Ni=LiN_i=L_i。用节点坐标展开:

Ni=12A(ai+bix+ciy),i=1,2,3N_i = \frac{1}{2A}(a_i + b_i x + c_i y),\quad i=1,2,3

系数(2×2\times 面积归一)由节点坐标定,轮换关系:

a1=x2y3−x3y2,  b1=y2−y3,  c1=x3−x2(下标 1→2→3 轮换)a_1=x_2 y_3-x_3 y_2,\; b_1=y_2-y_3,\; c_1=x_3-x_2\quad(\text{下标 }1\to2\to3\text{ 轮换})

③ 应变-位移矩阵 BB(关键)

位移 u=∑Niuiu=\sum N_i u_i,v=∑Niviv=\sum N_i v_i。应变 ε=L de\varepsilon=L\,d_e(de=[u1,v1,u2,v2,u3,v3]Td_e=[u_1,v_1,u_2,v_2,u_3,v_3]^T)。由于 NiN_i 是线性,∂Ni/∂x=bi/2A\partial N_i/\partial x=b_i/2A、∂Ni/∂y=ci/2A\partial N_i/\partial y=c_i/2A 都是常数 → 应变在单元内是常数 → 故名”常应变三角”:

B=12A[b10b20b300c10c20c3c1b1c2b2c3b3]B=\frac{1}{2A}\begin{bmatrix} b_1&0&b_2&0&b_3&0\\ 0&c_1&0&c_2&0&c_3\\ c_1&b_1&c_2&b_2&c_3&b_3 \end{bmatrix}

④ 单元刚度矩阵

由虚功 / 应变能,ke=∫ABTDB t dAk_e=\int_A B^TDB\,t\,\mathrm dA。因为 B,DB,D 都是常数(厚度 tt 也常数),积分解析:

  ke=BTDB t A  \boxed{\;k_e = B^T D B\,t\,A\;}

(6×66\times6 矩阵,对称正定)

⑤ CST 的特点与缺陷

  • ✅ 简单、形状适应(任意三角形拼合任意区域)
  • ❌ 应变/应力在单元内是常数 → 单元边界应力跳跃,需细分网格才有合理应力分布
  • ❌ 对弯曲问题过刚(over-stiff,位移偏小),要很密的网格才准

经验:CST 只在应力梯度小、或粗算时用;精细分析用 LST 或 Q4。

12.2 六节点三角形单元(LST,Linear Strain Triangle)

在 CST 基础上加 3 个边中节点(共 6 节点),用二次形函数 → 应变线性 → 精度远高于 CST。

二次形函数(用面积坐标表示):

角节点:  N1=L1(2L1−1),  N2=L2(2L2−1),  N3=L3(2L3−1)边中节点:  N4=4L1L2,  N5=4L2L3,  N6=4L3L1\begin{aligned} \text{角节点:}\;&N_1=L_1(2L_1-1),\;N_2=L_2(2L_2-1),\;N_3=L_3(2L_3-1)\\ \text{边中节点:}\;&N_4=4L_1L_2,\;N_5=4L_2L_3,\;N_6=4L_3L_1 \end{aligned}

(12 自由度,可表示完全二次多项式)。LST 应力在单元内线性变化,适合应力梯度大的区域。

12.3 四节点四边形单元(Q4,等参元)—— 工程主力

四边形是 2D 工程最常用单元。关键难点:任意四边形在物理坐标下形函数难写,等参元思路把它统一到自然坐标 ξ,η∈[−1,1]\xi,\eta\in[-1,1]。

① 自然坐标下的形函数(双线性)

标准正方形 4 角节点 (ξi,ηi)=(±1,±1)(\xi_i,\eta_i)=(\pm1,\pm1):

  Ni=14(1+ξiξ)(1+ηiη)  \boxed{\;N_i=\frac14(1+\xi_i\xi)(1+\eta_i\eta)\;}

② 等参映射(物理坐标也用同一套形函数)

x=∑i=14Ni(ξ,η) xi,y=∑i=14Ni(ξ,η) yix=\sum_{i=1}^{4} N_i(\xi,\eta)\,x_i,\qquad y=\sum_{i=1}^{4} N_i(\xi,\eta)\,y_i

“位移和坐标用同样的形函数” → 等参(isoparametric)。这一招让任意四边形(甚至曲边,用更多节点)都能统一处理。

③ 雅可比矩阵(自然→物理坐标变换)

J=[∂x∂ξ∂x∂η∂y∂ξ∂y∂η],∣J∣=det⁡JJ=\begin{bmatrix}\dfrac{\partial x}{\partial\xi}&\dfrac{\partial x}{\partial\eta}\\[4pt]\dfrac{\partial y}{\partial\xi}&\dfrac{\partial y}{\partial\eta}\end{bmatrix},\qquad |J|=\det J

导数变换(链式法则):

[∂∂x∂∂y]=J−1[∂∂ξ∂∂η]\begin{bmatrix}\dfrac{\partial}{\partial x}\\[2pt]\dfrac{\partial}{\partial y}\end{bmatrix}=J^{-1}\begin{bmatrix}\dfrac{\partial}{\partial\xi}\\[2pt]\dfrac{\partial}{\partial\eta}\end{bmatrix}

∣J∣|J| 必须为正(单元不能畸变到翻转),否则 J−1J^{-1} 出问题 → 这就是”网格质量(雅可比)“检查的由来。

④ 应变-位移矩阵 BB

BB 由 ∂Ni/∂x,∂Ni/∂y\partial N_i/\partial x, \partial N_i/\partial y 组成,通过 J−1J^{-1} 从自然坐标导数求得。BB 是 ξ,η\xi,\eta 的函数(不是常数) → 刚度积分不能用解析公式。

⑤ 单元刚度(必须数值积分)

ke=∫−11 ⁣∫−11BTDB ∣J∣ t dξ dη  ≈  ∑i∑jwiwj (BTDB)ij ∣J∣ij tk_e=\int_{-1}^{1}\!\int_{-1}^{1} B^T D B\,|J|\,t\,\mathrm d\xi\,\mathrm d\eta\;\approx\;\sum_i\sum_j w_i w_j\,(B^TDB)_{ij}\,|J|_{ij}\,t

这就必须用高斯积分(下一章)。

Q4 比 CST 准(应变线性,非常数)、比 LST 简单(4 节点),是平面问题首选。

十三、高斯数值积分

13.1 为什么需要数值积分?

等参元的 kek_e 积分在自然坐标域,被积函数 BTDB∣J∣B^TDB|J| 含形函数导数 + 雅可比,无法解析积分(尤其畸变单元)。必须数值积分。

高斯积分是有限元标配:它用最少的点达到最高代数精度。

13.2 一维高斯-勒让德积分

∫−11f(ξ) dξ  ≈  ∑i=1nwi f(ξi)\int_{-1}^{1} f(\xi)\,\mathrm d\xi \;\approx\; \sum_{i=1}^{n} w_i\,f(\xi_i)

选定积分点 ξi\xi_i 和权 wiw_i。nn 个高斯点能精确积分 2n−12n-1 阶多项式。

常用点表(背下来):

点数 nn积分点 ξi\xi_i权 wiw_i精确阶
100221(线性)
2±13≈±0.5774\pm\dfrac{1}{\sqrt3}\approx\pm0.57741,  11,\;13(三次)
30,  ±35≈±0.77460,\;\pm\sqrt{\dfrac35}\approx\pm0.774689,  59,59\dfrac89,\;\dfrac59,\dfrac595

13.3 多维(张量积)

二维 = 两个一维相乘:

∫−11 ⁣∫−11f dξ dη=∑i∑jwiwj f(ξi,ηj)\int_{-1}^{1}\!\int_{-1}^{1} f\,\mathrm d\xi\,\mathrm d\eta=\sum_i\sum_j w_i w_j\,f(\xi_i,\eta_j)
  • Q4 平面:2×2=42\times2=4 个高斯点(精确 3 阶,足够覆盖 Q4 的双线性)
  • H8 三维:2×2×2=82\times2\times2=8 个高斯点

13.4 积分阶的取舍:全积分 vs 减缩积分

  • 全积分(full):用够阶(如 Q4 用 2×22\times2)。精确,但薄板/细长结构的 Q4 会剪切锁死(shear locking,过刚,位移严重偏小)。
  • 减缩积分(reduced):降一阶(如 Q4 用 1×1=11\times1=1 点)。软化剪切、避免锁死、省算力,但可能引入零能模式/沙漏(hourglass) —— 一种不产生应变、刚度为零的虚假变形,需额外控制。

工程实践:选择性减缩积分(剪切项用减缩、弯曲项用全积分),或软件内置的沙漏控制。商业软件(ANSYS/Abaqus)默认处理,但用户要懂原理才能判结果。

十四、三维单元

三维问题每节点 3 自由度 (u,v,w)(u,v,w),应变向量 6 个分量 ε=[εx,εy,εz,γxy,γyz,γzx]T\varepsilon=[\varepsilon_x,\varepsilon_y,\varepsilon_z,\gamma_{xy},\gamma_{yz},\gamma_{zx}]^T,DD 是 6×66\times6。

14.1 四节点四面体(T4)

3D 版 CST。4 节点,12 自由度。

体积坐标 ζ1,ζ2,ζ3,ζ4\zeta_1,\zeta_2,\zeta_3,\zeta_4(类比三角形面积坐标,∑ζi=1\sum\zeta_i=1)。

形函数:Ni=ζiN_i=\zeta_i(线性)→ 常应变四面体。

  • ✅ 网格自动性极好(任意 3D 几何都能自动切四面体)
  • ❌ 精度低(常应变,过刚,像 CST),需很密网格

14.2 十节点二次四面体(T10)

T4 加边中节点,二次形函数 → 精度大幅提升。工程上复杂几何的首选(精度 + 自动性平衡)。ANSYS/Abaqus 自动网格默认生成 T10。

14.3 八节点六面体(H8,brick)

3D 版 Q4。8 节点,自然坐标 ξ,η,ζ∈[−1,1]\xi,\eta,\zeta\in[-1,1]。

形函数(三线性):

Ni=18(1+ξiξ)(1+ηiη)(1+ζiζ)N_i=\frac18(1+\xi_i\xi)(1+\eta_i\eta)(1+\zeta_i\zeta)

等参映射 + 3×33\times3 雅可比 + 2×2×2=82\times2\times2=8 点高斯积分。

  • ✅ 精度高(应变线性)
  • ❌ 网格划分难(需结构化/扫掠,复杂几何难自动生成六面体)

14.4 单元选择对比

单元网格自动性精度典型用途
T4 四面体极易低(过刚)不推荐粗用
T10 二次四面体易高复杂几何主力
H8 六面体难高规则几何、大模型
H20 二次六面体很难很高高精度关键件

经验:复杂几何用 T10,规则几何用 H8/H20。别用 T4 做精细分析(除非网格极密)。

十五、高阶单元与形函数族

15.1 提高精度的两条路(再强调)

  • h-收敛:加密网格(减小单元尺寸 hh)
  • p-收敛:提高单元阶数 pp(线性→二次→三次)

p-收敛收敛速率指数级(比 h 的代数级快),但每单元算力大。自适应 FEM 就是自动选 h 或 p。

15.2 形函数的两大族

① Lagrange 单元:形函数由完整多项式构造(含内部节点)。例:9 节点四边形(4 角 + 4 边中 + 1 中心)。一维形式:

Ni(ξ)=∏j≠iξ−ξjξi−ξjN_i(\xi)=\prod_{j\ne i}\frac{\xi-\xi_j}{\xi_i-\xi_j}

优点:完备性好;缺点:内部节点增加自由度。

② Serendipity 单元:只有角节点 + 边中节点,无内部节点(如 8 节点四边形)。形函数靠构造满足节点性质。优点:省内部节点(更经济),工程常用;缺点:高阶时完备性略差。

15.3 Pascal 三角形(完备性保证)

二维多项式按 Pascal 三角选择,保证完备性(包含所有低阶项)→ 满足收敛条件:

        1              ← 常数(刚体)
      x   y            ← 线性(刚体+常应变)
    x²  xy  y²         ← 二次
  x³ x²y xy² y³        ← 三次
  • CST:取到线性(3 项,x,yx,y + 常数)
  • Q4:双线性(1,x,y,xy1,x,y,xy,4 项)
  • LST:完全二次(6 项)
  • 完备性 → 单元能精确表示常应变 → 收敛保证(patch test 满足)

十六、梁单元(弯曲问题)

杆只受轴力;梁受弯矩 + 剪力(横向弯曲)。梁单元是结构分析另一大基础。

16.1 欧拉-伯努利梁假设

  • 变形前垂直中性轴的截面,变形后仍平面且垂直中性轴(忽略剪切变形)
  • 适合细长梁(L/h>10L/h>10)

挠度 w(x)w(x),转角 θ=dw/dx\theta=\mathrm dw/\mathrm dx。每节点 2 个自由度:横向位移 ww + 转角 θ\theta。

16.2 Hermite 三次形函数(C¹ 连续)

梁要求位移和转角都连续(C¹ 连续),普通 Lagrange 形函数(只保证 C⁰ 位移连续)不够。要用 Hermite 多项式。

2 节点 4 自由度 [w1,θ1,w2,θ2][w_1,\theta_1,w_2,\theta_2],令 s=x/L∈[0,1]s=x/L\in[0,1],三次 Hermite 形函数为(N1,N3N_1,N_3 对应位移、N2,N4N_2,N_4 对应转角):

  N1=1−3s2+2s3N2=L(s3−2s2+s)N3=3s2−2s3N4=L(s3−s2)  \boxed{\; \begin{aligned} N_1&=1-3s^2+2s^3\\ N_2&=L(s^3-2s^2+s)\\ N_3&=3s^2-2s^3\\ N_4&=L(s^3-s^2) \end{aligned}\;}

性质:N1,N3N_1,N_3 控制位移(节点处=1,导数=0);N2,N4N_2,N_4 控制转角(位移=0,导数=1)。两单元交接处 ww 和 θ=dw/dx\theta=\mathrm dw/\mathrm dx 都连续 → C¹ 连续。

16.3 梁单元刚度矩阵

梁的应变能 U=12∫0LEI (w′′)2 dxU=\tfrac12\int_0^L EI\,(w'')^2\,\mathrm dx(曲率 w′′w'' 是变形核心)。代入 w=∑Nidiw=\sum N_i d_i 积分,得 4×44\times4 刚度矩阵(自由度 [w1,θ1,w2,θ2][w_1,\theta_1,w_2,\theta_2]):

  ke=EIL3[126L−126L6L4L2−6L2L2−12−6L12−6L6L2L2−6L4L2]  \boxed{\; k_e=\frac{EI}{L^3}\begin{bmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^2 & -6L & 2L^2\\ -12 & -6L & 12 & -6L\\ 6L & 2L^2 & -6L & 4L^2 \end{bmatrix}\;}

读法:k11=12EI/L3k_{11}=12EI/L^3 是”在节点 1 加单位横向位移(转角固定)所需的力”;k22=4EI?/Lk_{22}=4EI^? /L… 含 L2L^2 的项对应转角自由度。

16.4 Timoshenko 梁(厚梁)

当梁不够细长(L/hL/h 小),剪切变形不可忽略,假设”截面保持平面但不垂直中性轴”,转角 θ≠dw/dx\theta\ne\mathrm dw/\mathrm dx。需剪切修正,适合短粗梁、复合材料。

十七、动力学分析简介

静力是 Ku=FKu=F;动力问题(振动、冲击)加惯性和阻尼:

  Mu¨+Cu˙+Ku=F(t)  \boxed{\;M\ddot u + C\dot u + Ku = F(t)\;}
  • MM:质量矩阵(∫ρNTN dV\int \rho N^TN\,\mathrm dV,一致质量;或对角化的集中质量)
  • CC:阻尼矩阵(常用 Rayleigh 阻尼 C=αM+βKC=\alpha M+\beta K)
  • u¨,u˙\ddot u,\dot u:加速度、速度

17.1 模态分析(自由振动)

无阻尼自由振动 Mu¨+Ku=0M\ddot u+Ku=0,设 u=ϕ eiωtu=\phi\,e^{i\omega t},得特征值问题:

(K−ω2M)ϕ=0(K-\omega^2 M)\phi = 0

解出固有频率 ωi\omega_i(和频率 fi=ωi/2πf_i=\omega_i/2\pi)与振型 ϕi\phi_i。这是结构动力特性的”指纹”(避共振、抗震设计基础)。

17.2 瞬态动力响应

F(t)F(t) 随时间变化,需时间积分:

  • 隐式(Newmark-β\beta):大步长稳定,适合低速/准静
  • 显式(中心差分):步长受 CFL 限制(很小),适合冲击/爆炸/碰撞(ANSYS LS-DYNA、Abaqus/Explicit)

十八、非线性简介

前面都是线性(小变形、线弹性,KK 常数)。真实问题常非线性:KK 随位移/状态变,Ku=FKu=F 不能直接解,需迭代。

K(u) u=F⟹牛顿-拉夫逊迭代K(u)\,u = F \quad\Longrightarrow\quad \text{牛顿-拉夫逊迭代}

18.1 三类非线性

类型来源例子
材料非线性σ=σ(ε)\sigma=\sigma(\varepsilon) 非线性塑性(屈服后)、橡胶超弹性、蠕变
几何非线性大变形,应变-位移非线性薄板后屈曲、缆索、橡胶大变形
接触非线性边界条件随变形变(接触/分离)装配、碰撞、螺栓

18.2 牛顿-拉夫逊迭代

  1. 猜初值 u0u_0
  2. 算残差 R=F−K(un)unR=F-K(u_n)u_n
  3. 算切线刚度 KT=∂(Ku)/∂uK_T=\partial(Ku)/\partial u
  4. 修正 Δu=KT−1R\Delta u=K_T^{-1}R,un+1=un+Δuu_{n+1}=u_n+\Delta u
  5. 重复直到 RR 收敛(足够小)

非线性比线性贵得多(每次迭代解一次线性方程组),且可能不收敛(需调步长/方法)。

完结:CAE/有限元知识全景

至此,CAE 与有限元法的基础知识体系完整覆盖:

物理层:  弹性力学(应力/应变/本构/三大方程)
   ↓ 强形式 → 弱形式(虚功/势能)
数学层:  加权残值(Galerkin)→ Ku=F
   ↓
离散层:  形函数(插值)→ 单元方程
   ↓
单元库:  杆 / 桁架 / 梁(Hermite)
        二维:CST / LST / Q4(等参)
        三维:T4 / T10 / H8 / H20
        高阶:Lagrange / Serendipity(p-收敛)
   ↓
求解层:  组装(直接刚度法)→ 边界条件 → 解 Ku=F
        数值积分(高斯)→ 网格收敛(h/p)
   ↓
拓展层:  动力学(Mü̈+Cu̇+Ku=F)→ 非线性(材料/几何/接触)

后续深入方向(超出基础课范围):显式动力学、断裂力学/疲劳、多物理场耦合(热-力、流-固)、拓扑优化、meshless(无网格)、isogeometric(等几何)、AI 辅助建模……

参考资料

这是一篇持续更新的课程笔记。基础部分(〇~十八)已完整覆盖 CAE/有限元入门;后续若涉及专题(断裂/疲劳/多物理场/优化),另开新篇。

参考资料

持续更新笔记。下次课后追加新章节,更新顶部 modDatetime 即可置顶。

编辑此页面
Previous Post
Emscripten 在 Windows 下的安装