Ch1 · 连续介质力学与线性有限元分析

📖 1.1 引言：为什么要学习非线性 FEM？

有限元方法 (FEM) 求解工程中的偏微分方程（结构力学 → 传热 → 流体力学 → 电磁）。线性结构系统已经有成熟的建模和求解方法，但非线性系统（弹塑性、大变形、接触、动力响应）需要基于问题特征的不同建模和不同求解方法。

本书目的：直接、清晰地介绍非线性 FEM 分析过程，重点 4 件事：

✅ 代表性非线性问题（不穷举所有理论）
✅ 详细理论推导（含每一步的物理解释）
✅ 求解过程（Newton-Raphson、弧长法、显式 / 隐式）
✅ MATLAB 实现（让初学者能上手）

读者对象：机械 / 土木 / 航天 / 生物医学 / 工业工程的研究生、研究人员、设计工程师。

📐 1.2 向量与张量分析

1.2.1 基底与 Einstein 求和约定

三维向量 = 三个基向量的线性组合。Einstein 求和约定：

(1.1)

u = u_je_j = u₁e₁ + u₂e₂ + u₃e₃

重复指标自动求和（j = 1,2,3 自动求和）。一个指标在单项中只能出现两次。

Kronecker δ 与内积

(1.2)

u · v = u_iv_i = u₁v₁ + u₂v₂ + u₃v₃

(1.3)

δ_ij = 1 if i = j, 0 if i ≠ j （Kronecker δ）

💡 例 1.1：向量运算

u = (3, 4)，v = (1, 2)：

u · v = 3·1 + 4·2 = 11
|u| = √(3² + 4²) = 5
由 δ_ij：v_j δ_ij = v_i（重复指标替换）

置换符号 e_ijk（Levi-Civita）

3D 叉积需要：

(1.20)

e_ijk = +1 (偶置换), −1 (奇置换), 0 (有重复)

(1.22)

Wu = w × u （反对称张量 W 关联向量 w）

1.2.2 向量与张量微积分

梯度算子：

(1.23)

∇ = ∂/∂x_i · e_i

标量场的梯度：∇φ（向量）
向量场的梯度：∇u（2 阶张量）
散度：∇ · u = ∂u_i/∂x_i（标量）
Laplace：∇² = ∂²/∂x_i²

1.2.3 散度定理（重要积分定理）

(1.30)

∫_Ω ∇ · F dΩ = ∫_Γ F · n dΓ

体积分 → 面积分。是 FEM 弱形式推导（虚功原理）的核心数学基础。

图 1.x 散度定理：体积分等于边界通量

图 1.x 向量内积与 Kronecker δ

图 1.x 梯度算子与向量微积分

💪 1.3 应力与应变

1.3.1 应力张量

应力 = 单位面积上的内力。3D 中需要2 阶应力张量 σ_ij 描述：

(1.36)

σ = ⎡σ₁₁ σ₁₂ σ₁₃⎤
⎣σ₁₂ σ₂₂ σ₂₃⎦
⎡σ₁₃ σ₂₃ σ₃₃⎤

对称 σ_ij = σ_ji（由动量矩守恒，无体力偶时）。共 6 个独立分量：3 正应力 + 3 剪应力。

⚠ 应力的"正方向"约定

应力分量 σ_ij 的正负号 = 面的正方向 × 力的正方向。
正法线 + 负力 → 负分量。剪应力在正面 + 沿正轴 → 正。

Voigt 标记

把 6 个独立分量写为 6×1 伪向量：

(1.37)

{σ} = ⎡σ₁₁, σ₂₂, σ₃₃, σ₁₂, σ₂₃, σ₁₃⎤^T

主应力与不变量

解 det(σ - λI) = 0 → 3 个主应力 σ₁, σ₂, σ₃（对应特征向量方向无剪应力）。

I₁ = tr(σ) = σ_kk

I₂ = (1/2)[(tr σ)² − tr(σ²)]

I₃ = det(σ)

3 个不变量在坐标旋转下保持不变。

图 1.5 应力张量正方向约定

1.3.2 应变张量

小变形假设下的无穷小应变张量：

(1.40)

ε_ij = (1/2)(∂u_i/∂x_j + ∂u_j/∂x_i)

对角元素 ε_ii：正应变（拉伸 / 压缩）
非对角元素 ε_ij：剪应变（角变形的一半）

1.3.3 本构关系（Hooke 定律）

各向同性线弹性：

(1.46)

σ_ij = λ ε_kk δ_ij + 2μ ε_ij

其中 λ 和 μ 是 Lamé 常数，与工程常数 E（杨氏模量）、ν（泊松比）的关系：

λ = Eν / [(1+ν)(1−2ν)]

μ = G = E / [2(1+ν)]

用 E、ν 表达：

(1.48)

σ_ij = [E/(1+ν)] · ε_ij + [Eν/((1+ν)(1−2ν))] · ε_kk · δ_ij

图 1.x 各向同性线弹性本构

💡 例 1.2：单轴拉伸应力

钢杆 E = 200 GPa, ν = 0.3，单轴应力 σ₁₁ = 100 MPa（其他方向无应力）：

横向应变：ε₂₂ = ε₃₃ = −ν · σ₁₁/E = −0.3 × 100/200000 = −1.5e−4
轴向应变：ε₁₁ = σ₁₁/E = 100/200000 = 5e−4
体积应变：ε_ii = ε₁₁ + ε₂₂ + ε₃₃ = 5e−4 − 2·1.5e−4 = 2e−4（很小，因 ν < 0.5）

图 1.x 单轴拉伸应力-应变

⚖️ 1.4 连续体力学

1.4.1 边值问题 (BVP)

固体力学 BVP 包含 3 部分：

平衡方程（动量守恒）：σ_ij,j + f_i = 0
几何关系（应变-位移）：ε_ij = (1/2)(u_i,j + u_j,i)
本构关系：σ = f(ε)（如 Hooke 定律）

边界条件：

位移边界（Dirichlet）：u_i = ū_i 在 Γ_u
力边界（Neumann）：σ_ij n_j = t̄_i 在 Γ_t
Γ_u ∪ Γ_t = ∂Ω（互补覆盖）

1.4.2 最小势能原理

对保守系统，真实位移场使总势能 Π 取驻值：

(1.66)

δΠ = 0

总势能 = 应变能 − 外力功：

(1.65)

Π = ∫_Ω (1/2) σ : ε dΩ − ∫_Ω f · u dΩ − ∫_{Γ_t} t̄ · u dΓ

优点：只要求一阶导，比平衡方程（要求二阶导）更宽松 → FEM 直接基于此。

1.4.3 虚功原理

对任意虚位移 δu（在位移边界上为零）：

(1.70)

∫_Ω σ : δε dΩ = ∫_Ω f · δu dΩ + ∫_{Γ_t} t̄ · δu dΓ

虚功原理是变分形式，不要求势能存在 → 适合任何材料（不仅保守系统）。
对保守弹性系统，虚功 ≡ 最小势能。

图 1.x 虚功原理的几何解释

🔧 1.5 线性系统有限元分析

1.5.1 有限元近似

把连续体划分为有限个单元，单元内位移用形函数插值：

(1.78)

u(x) ≈ N_a(x) u_a

其中 N_a 是形函数，u_a 是节点位移（对单元节点 a 求和）。

1.5.2 一维杆单元有限元方程

杆单元（长度 L，面积 A，弹性模量 E）：

(1.85)

K^e = (EA/L) · ⎡ 1 −1 ⎤
⎣−1 1 ⎦

组装：K U = F。

1.5.3 3D 实体单元

八节点六面体（hex8）单元刚度：

(1.137)

K^e = ∫_{V_e} B^T C B dV

其中 B = LN（应变-位移矩阵），C = 材料弹性矩阵，LN 是形函数梯度。

1.5.4 MATLAB 实现

第 1.5.5 节给出 3D 8 节点六面体的 MATLAB 代码（见 PDF p50-58）。

图 1.x 有限元离散与组装

图 1.x 2D 平面应变四边形单元

💡 例 1.5：3D 单元应力计算

8 节点六面体单元，节点位移已知，形函数 N 已知：

① 计算每个高斯点的 Jacobian J = dx/dξ
② 计算 B = LN · J^-1
③ 应变 ε = B · U
④ 应力 σ = C · ε
⑤ 单元刚度 K^e = ∫ B^T C B |J| dξdηdζ

图 1.x 3D 六面体单元网格

📌 1.6 本章总结

✅ 向量与张量基础（Einstein 求和、Kronecker δ、Levi-Civita 置换符号）
✅ 应力（对称 2 阶张量、6 个独立分量、Voigt 标记、主应力）
✅ 应变（无穷小应变张量、正应变 / 剪应变）
✅ 本构关系（各向同性线弹性 Hooke 定律、λ/μ 与 E/ν）
✅ 边值问题（平衡 + 几何 + 本构 + 边界）
✅ 最小势能原理（变分形式）
✅ 虚功原理（不要求势能存在）
✅ 线性 FEM（形函数、单元刚度、组装、3D hex8 实现）

第 2 章将基于这些基础，扩展到非线性系统（Newton-Raphson 求解）。

图 1.x 本章关键概念

📝 1.7 习题（PDF p78-86）

建议先做以下题目：

1.1 向量求和约定（1 维、2 维、3 维）
1.2 应力张量对称性证明（从动量矩）
1.3 主应力与不变量计算
1.4 单轴拉伸的应力-应变关系
1.5 置换符号 e_ijk 性质
1.10 单元刚度矩阵推导（1D 杆 → 2D 桁架 → 3D 梁）
1.20 3D 实体单元的有限元方程推导
1.30 用 MATLAB 实现 8 节点六面体单元

📚 本章 PDF 关键页面

点击放大看原图，对照学习。

第 1 章 · 连续介质力学与线性有限元分析