第 4 有限元法
4.1 有限元法理论基础
有限元法（Finite Element Method, FEM）是一种强大的数值分析方法，广泛应用于求解工程和科学中的偏微分方程问题。它由Richard Courant于1943年首次提出，经过几十年的发展，已成为现代数值计算的重要工具之一。
4.1.1 历史发展与应用背景
有限元法的发展经历了几个重要阶段：
- 理论奠基期（1943-1960）：Courant提出了使用分片多项式和变分方法求解偏微分方程的思想
- 工程应用期（1960-1980）：在结构力学和航空航天领域得到广泛应用
- 理论完善期（1980-2000）：数学理论体系日趋完善，收敛性和误差分析取得重大进展
- 多学科融合期（2000至今）：在地球科学、生物医学、材料科学等领域得到广泛应用
在地球科学中，有限元法特别适用于： - 地下水流动模拟 - 地质构造变形分析 - 地震波传播计算 - 地热传输模拟 - 污染物扩散预测
4.1.2 基本思想与哲学
有限元法的核心思想体现了”分而治之”的哲学理念，将复杂的连续问题分解为简单的离散问题。这种思想具有以下特点：
- 局部化原理：复杂的全局问题可以通过求解简单的局部问题来近似
- 分片逼近：使用简单的多项式函数在每个小区域内逼近真解
- 变分基础：基于能量原理或变分原理，确保解的物理合理性
- 系统性方法：提供了一套系统的、标准化的求解流程
4.1.3 详细求解步骤
有限元法的完整求解过程包括以下七个核心步骤：
4.1.3.1 1. 域离散化（Domain Discretization）
将连续的求解域 \(\Omega\) 划分为 \(N_e\) 个不重叠的子域（单元）： \[\Omega = \bigcup_{e=1}^{N_e} \Omega_e, \quad \Omega_i \cap \Omega_j = \emptyset \text{ (when } i \neq j \text{)}\]
单元类型的选择原则： - 一维问题：线段单元（2节点或3节点） - 二维问题：三角形单元或四边形单元 - 三维问题：四面体单元或六面体单元
4.1.3.2 2. 形函数构造（Shape Function Construction）
在每个单元内定义形函数 \(N_i(x)\)，满足： - 插值性质：\(N_i(x_j) = \delta_{ij}\)（Kronecker delta） - 完备性质：能够精确表示多项式函数 - 连续性质：在单元边界上保持连续性
4.1.3.5 5. 全局矩阵组装（Global Matrix Assembly）
将所有单元的贡献组装成全局矩阵： \[\mathbf{K} = \sum_{e=1}^{N_e} \mathbf{A}^e \mathbf{K}^e (\mathbf{A}^e)^T\]
其中 \(\mathbf{A}^e\) 是单元组装矩阵。
4.1.4 变分原理的深入理解
4.1.4.1 变分原理的数学基础
变分原理是有限元法的理论基石，它建立了偏微分方程与泛函极值问题之间的等价关系。这种等价性不仅为数值求解提供了新的途径，更重要的是为解的存在性和唯一性提供了理论保证。
泛函的定义： 泛函是从函数空间到实数的映射。对于函数 \(u(x)\)，泛函 \(J[u]\) 将函数 \(u\) 映射到一个实数值： \[J[u] = \int_\Omega F(x, u, \nabla u) dx\]
其中 \(F\) 是给定的函数，\(\nabla u\) 是 \(u\) 的梯度。
4.1.4.2 从强形式到弱形式的详细推导
考虑一般的二阶椭圆型偏微分方程： \[-\nabla \cdot (k(x) \nabla u) + c(x) u = f(x) \quad \text{in} \quad \Omega\]
配以边界条件： \[u = g \quad \text{on} \quad \Gamma_D \quad \text{(Dirichlet边界)}\] \[k \frac{\partial u}{\partial n} = h \quad \text{on} \quad \Gamma_N \quad \text{(Neumann边界)}\]
其中 \(\Gamma_D \cup \Gamma_N = \partial\Omega\)，\(\partial\Omega\) 是 \(\Omega\) 的边界。
步骤1：构造测试函数空间
定义测试函数空间： \[V_0 = \{v \in H^1(\Omega) : v = 0 \text{ on } \Gamma_D\}\]
其中 \(H^1(\Omega)\) 是Sobolev空间，包含所有在 \(\Omega\) 上平方可积且弱导数也平方可积的函数。
步骤2：建立弱形式
将强形式方程乘以测试函数 \(v \in V_0\) 并在 \(\Omega\) 上积分： \[\int_\Omega [-\nabla \cdot (k \nabla u) + c u] v \, dx = \int_\Omega f v \, dx\]
步骤3：应用Green公式（分部积分）
对第一项应用Green公式： \[-\int_\Omega v \nabla \cdot (k \nabla u) dx = \int_\Omega k \nabla u \cdot \nabla v \, dx - \int_{\partial\Omega} v k \frac{\partial u}{\partial n} ds\]
由于 \(v = 0\) 在 \(\Gamma_D\) 上，边界积分只在 \(\Gamma_N\) 上非零： \[-\int_{\partial\Omega} v k \frac{\partial u}{\partial n} ds = -\int_{\Gamma_N} v h \, ds\]
步骤4：得到弱形式
最终的弱形式为：寻找 \(u \in V_g = \{w \in H^1(\Omega) : w = g \text{ on } \Gamma_D\}\)，使得对所有 \(v \in V_0\)：
\[\int_\Omega k \nabla u \cdot \nabla v \, dx + \int_\Omega c u v \, dx = \int_\Omega f v \, dx + \int_{\Gamma_N} h v \, ds\]
这可以写成抽象形式： \[a(u,v) = L(v)\]
其中： - \(a(u,v) = \int_\Omega k \nabla u \cdot \nabla v \, dx + \int_\Omega c u v \, dx\) （双线性形式） - \(L(v) = \int_\Omega f v \, dx + \int_{\Gamma_N} h v \, ds\) （线性形式）
4.1.4.3 变分原理的物理意义
变分原理具有深刻的物理意义：
- 能量最小原理：对于弹性力学问题，弱形式对应于势能的最小化
- 虚功原理：弱形式表达了虚位移做的虚功等于外力做的虚功
- 守恒定律：弱形式自然地体现了各种守恒定律（质量、动量、能量等）
4.1.5 形函数理论的完整体系
形函数（Shape Functions）是有限元法的核心概念，它们不仅定义了单元内的插值关系，更是连接离散节点值与连续函数空间的桥梁。深入理解形函数的性质和构造方法，对于掌握有限元法的精髓至关重要。
4.1.5.1 形函数的基本性质
形函数必须满足以下四个基本性质：
1. 插值性质（Interpolation Property） \[N_i(x_j) = \delta_{ij} = \begin{cases} 1 & \text{if } i = j \\ 0 & \text{if } i \neq j \end{cases}\]
这确保了节点值的准确插值。
2. 单位分解性质（Partition of Unity） \[\sum_{i=1}^{n} N_i(x) = 1, \quad \forall x \in \Omega_e\]
这保证了常数函数能够被精确表示。
3. 完备性性质（Completeness Property） 形函数集合能够精确表示完全多项式，至少到某个给定阶数 \(p\)： \[\mathcal{P}_p \subset \text{span}\{N_1, N_2, \ldots, N_n\}\]
其中 \(\mathcal{P}_p\) 是阶数不超过 \(p\) 的多项式空间。
4. 连续性性质（Continuity Property） 在单元边界上，形函数保持 \(C^0\) 连续性，即函数值连续但导数可能不连续。
4.1.5.2 一维形函数的详细构造
线性形函数（Linear Shape Functions）
对于标准单元 \(\xi \in [-1, 1]\)： \[N_1(\xi) = \frac{1-\xi}{2}, \quad N_2(\xi) = \frac{1+\xi}{2}\]
其导数为： \[\frac{dN_1}{d\xi} = -\frac{1}{2}, \quad \frac{dN_2}{d\xi} = \frac{1}{2}\]
二次形函数（Quadratic Shape Functions）
对于三节点单元，节点位于 \(\xi = -1, 0, 1\)： \[N_1(\xi) = \frac{\xi(\xi-1)}{2}, \quad N_2(\xi) = 1-\xi^2, \quad N_3(\xi) = \frac{\xi(\xi+1)}{2}\]
其导数为： \[\frac{dN_1}{d\xi} = \xi - \frac{1}{2}, \quad \frac{dN_2}{d\xi} = -2\xi, \quad \frac{dN_3}{d\xi} = \xi + \frac{1}{2}\]
三次形函数（Cubic Shape Functions）
对于四节点单元，可以使用Hermite插值： \[N_1(\xi) = \frac{1}{16}(1-\xi)^2(2+\xi)\] \[N_2(\xi) = \frac{1}{16}(1-\xi)^2(1+\xi)\] \[N_3(\xi) = \frac{1}{16}(1+\xi)^2(2-\xi)\] \[N_4(\xi) = \frac{1}{16}(1+\xi)^2(1-\xi)\]
4.1.5.3 Lagrange插值基础
形函数的构造通常基于Lagrange插值理论。对于 \(n+1\) 个节点 \(\xi_0, \xi_1, \ldots, \xi_n\)，Lagrange插值多项式为：
\[L_i(\xi) = \prod_{j=0, j \neq i}^{n} \frac{\xi - \xi_j}{\xi_i - \xi_j}\]
这自然满足插值性质：\(L_i(\xi_j) = \delta_{ij}\)。
4.1.5.4 等参变换理论
坐标变换
从物理坐标 \(x\) 到参考坐标 \(\xi\) 的变换： \[x = \sum_{i=1}^{n} x_i N_i(\xi)\]
其中 \(x_i\) 是节点的物理坐标。
雅可比矩阵
变换的雅可比矩阵为： \[J = \frac{dx}{d\xi} = \sum_{i=1}^{n} x_i \frac{dN_i}{d\xi}\]
对于一维情况，雅可比行列式 \(|J| = J\)。
导数变换
物理空间中的导数通过链式法则计算： \[\frac{dN_i}{dx} = \frac{dN_i}{d\xi} \frac{d\xi}{dx} = \frac{1}{J} \frac{dN_i}{d\xi}\]
4.1.5.5 二维形函数
三角形单元
对于三角形单元，使用面积坐标（重心坐标）\(L_1, L_2, L_3\)： \[N_1 = L_1, \quad N_2 = L_2, \quad N_3 = L_3\]
其中 \(L_1 + L_2 + L_3 = 1\)。
四边形单元
对于双线性四边形单元： \[N_1(\xi, \eta) = \frac{(1-\xi)(1-\eta)}{4}\] \[N_2(\xi, \eta) = \frac{(1+\xi)(1-\eta)}{4}\] \[N_3(\xi, \eta) = \frac{(1+\xi)(1+\eta)}{4}\] \[N_4(\xi, \eta) = \frac{(1-\xi)(1+\eta)}{4}\]
4.1.5.6 形函数的数值积分
形函数涉及的积分通常通过高斯积分法计算：
一维高斯积分 \[\int_{-1}^{1} f(\xi) d\xi \approx \sum_{i=1}^{n} w_i f(\xi_i)\]
其中 \(w_i\) 是权重，\(\xi_i\) 是积分点。
常用高斯积分点
- 1点积分：\(\xi_1 = 0\), \(w_1 = 2\)
- 2点积分：\(\xi_{1,2} = \pm\frac{1}{\sqrt{3}}\), \(w_{1,2} = 1\)
- 3点积分：\(\xi_1 = 0, \xi_{2,3} = \pm\sqrt{\frac{3}{5}}\), \(w_1 = \frac{8}{9}, w_{2,3} = \frac{5}{9}\)
4.2 数学公式论证
4.2.1 变分形式的推导
考虑一维热传导方程：
\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \quad \text{in} \quad (0,L) \times (0,T)\]
其中 \(u(x,t)\) 是温度，\(\alpha\) 是热扩散系数，\(L\) 是杆长，\(T\) 是总时间。
边界条件： \[u(0,t) = u_0, \quad u(L,t) = u_L\]
初始条件： \[u(x,0) = u_{ic}(x)\]
将方程乘以测试函数 \(v(x)\) 并在空间域上积分：
\[\int_0^L v(x) \frac{\partial u}{\partial t} dx = \alpha \int_0^L v(x) \frac{\partial^2 u}{\partial x^2} dx\]
对右边进行分部积分：
\[\int_0^L v(x) \frac{\partial u}{\partial t} dx = -\alpha \int_0^L \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} dx + \alpha \left[v(x) \frac{\partial u}{\partial x}\right]_0^L\]
由于边界条件已知，边界项可以处理，得到弱形式：
\[\int_0^L v(x) \frac{\partial u}{\partial t} dx + \alpha \int_0^L \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} dx = 0\]
4.2.2 有限元离散化
将空间域 \([0,L]\) 划分为 \(n\) 个单元，每个单元长度为 \(h = L/n\)。在每个单元 \([x_i, x_{i+1}]\) 上，使用线性形函数：
\[u^h(x) = \sum_{j=1}^{n+1} u_j N_j(x)\]
其中 \(u_j\) 是节点 \(j\) 处的温度值，\(N_j(x)\) 是节点 \(j\) 的形函数。
将弱形式离散化，得到：
\[\sum_{i=1}^n \int_{x_i}^{x_{i+1}} N_j(x) \frac{\partial u^h}{\partial t} dx + \alpha \sum_{i=1}^n \int_{x_i}^{x_{i+1}} \frac{\partial N_j}{\partial x} \frac{\partial u^h}{\partial x} dx = 0\]
这可以写成矩阵形式：
\[M \frac{d\mathbf{u}}{dt} + K\mathbf{u} = \mathbf{0}\]
其中 \(M\) 是质量矩阵，\(K\) 是刚度矩阵，\(\mathbf{u}\) 是节点温度向量。
4.2.3 单元矩阵的构造
对于线性单元，质量矩阵和刚度矩阵的元素为：
\[M_{ij} = \int_{x_i}^{x_{i+1}} N_i(x) N_j(x) dx\]
\[K_{ij} = \alpha \int_{x_i}^{x_{i+1}} \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial x} dx\]
对于均匀网格，这些积分可以解析计算：
\[M = \frac{h}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}\]
\[K = \frac{\alpha}{h} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\]
4.3 一维热传导问题的深入分析与FEM求解
4.3.1 工程背景与实际意义
热传导问题在地球科学和工程实践中具有重要意义。本节以铁棍热传导为例，展示有限元法在解决实际工程问题中的应用。这类问题在以下领域具有广泛应用：
地球科学应用： - 地热能开发：地下岩层的热传导分析 - 冻土工程：永久冻土层的热力学行为 - 地质勘探：地下温度场分布预测 - 环境工程：土壤污染热修复技术
工程应用： - 建筑节能：墙体保温材料设计 - 电子散热：芯片热管理系统 - 材料加工：金属热处理工艺 - 核工程：反应堆热传递分析
4.3.2 详细问题描述
考虑一根长度为 \(L = 1.0\) m 的均匀铁棍，其物理参数如下：
材料参数： - 热传导系数：\(k = 50\) W/(m·K) - 密度：\(\rho = 7850\) kg/m³ - 比热容：\(c = 460\) J/(kg·K) - 热扩散系数：\(\alpha = \frac{k}{\rho c} = 1.384 \times 10^{-5}\) m²/s
几何参数： - 长度：\(L = 1.0\) m - 横截面积：\(A = 1.0\) cm²（假设为细棍，忽略横向传热）
边界条件： - 左端（\(x = 0\)）：恒温 \(u(0,t) = 100°C\) - 右端（\(x = L\)）：恒温 \(u(L,t) = 50°C\)
初始条件： - 初始温度：\(u(x,0) = 25°C\)（室温）
4.3.3 CLAMS方法的系统化建模
CLAMS方法为复杂物理问题提供了系统化的建模框架，确保建模过程的完整性和科学性。
4.3.3.1 1. 概念描述 (Conceptual Model)
物理系统描述： 一维热传导系统，包含以下关键要素： - 系统边界：铁棍的几何边界 - 控制体积：整个铁棍内部 - 热源/热汇：两端的恒温边界条件 - 传热机制：纯导热（忽略对流和辐射）
简化假设的物理合理性： 1. 一维假设：当铁棍的长度远大于其直径时，横向温度梯度可忽略 2. 均匀材料：工业铁棍在宏观尺度上可视为均匀介质 3. 恒定物性：在给定温度范围内，铁的热物性变化较小
4.3.3.2 2. 物理定律 (Laws of Physics)
傅里叶热传导定律： \[\mathbf{q} = -k \nabla u\]
对于一维情况： \[q_x = -k \frac{\partial u}{\partial x}\]
其中 \(q_x\) 是单位面积的热流密度（W/m²）。
能量守恒定律： 基于控制体积的能量平衡： \[\frac{\partial}{\partial t}(\rho c u) + \nabla \cdot \mathbf{q} = S\]
其中 \(S\) 是内热源项。对于无内热源的一维情况： \[\rho c \frac{\partial u}{\partial t} + \frac{\partial q_x}{\partial x} = 0\]
热力学第二定律： 热量自发地从高温区域传向低温区域，这确保了问题解的物理合理性。
4.3.3.3 3. 假设条件 (Assumptions)
几何假设： 1. 一维传热：\(\frac{\partial u}{\partial y} = \frac{\partial u}{\partial z} = 0\) 2. 直杆假设：忽略几何形状对传热的影响 3. 边界绝热：侧表面无热交换
物理假设： 4. 恒定物性：\(k\)、\(\rho\)、\(c\) 为常数 5. 无内热源：\(S = 0\) 6. 线性传热：忽略非线性效应（如温度相关的物性）
数值假设： 7. 连续介质：忽略分子尺度的不连续性 8. 准静态过程：忽略惯性效应
假设的有效性分析： - 对于长径比 > 10 的铁棍，一维假设误差 < 5% - 在0-200°C温度范围内，铁的热物性变化 < 10% - 忽略对流和辐射在低温差情况下是合理的
4.3.3.4 4. 数学推导 (Mathematical Derivation)
控制方程的推导：
将傅里叶定律代入能量守恒方程： \[\rho c \frac{\partial u}{\partial t} - \frac{\partial}{\partial x}\left(k \frac{\partial u}{\partial x}\right) = 0\]
对于恒定物性： \[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\]
其中热扩散系数 \(\alpha = \frac{k}{\rho c}\) 具有量纲 [m²/s]。
边界条件的数学表达： - 第一类边界条件（Dirichlet）： \[u(0,t) = u_0 = 100°C\] \[u(L,t) = u_L = 50°C\]
初始条件： \[u(x,0) = u_{ic}(x) = 25°C\]
无量纲化处理： 为了便于数值计算和参数分析，引入无量纲变量： - 无量纲坐标：\(\xi = \frac{x}{L}\) - 无量纲温度：\(\theta = \frac{u - u_L}{u_0 - u_L}\) - 无量纲时间：\(\tau = \frac{\alpha t}{L^2}\)
无量纲控制方程： \[\frac{\partial \theta}{\partial \tau} = \frac{\partial^2 \theta}{\partial \xi^2}\]
边界条件：\(\theta(0,\tau) = 1\)，\(\theta(1,\tau) = 0\) 初始条件：\(\theta(\xi,0) = \frac{25-50}{100-50} = -0.5\)
4.3.3.5 5. 求解策略 (Solver Strategy)
有限元离散化策略：
- 空间离散化：
- 单元类型：一维线性单元
- 网格密度：根据Peclet数确定
- 单元数量：\(N_e = 20\)（经收敛性分析确定）
- 时间离散化：
- 时间积分格式：隐式Euler法（确保无条件稳定）
- 时间步长：\(\Delta t = 100\) s（满足精度要求）
- 总时间：\(T = 10000\) s（达到稳态）
- 线性方程组求解：
- 矩阵类型：对称正定稀疏矩阵
- 求解方法：直接法（LU分解）
- 边界条件处理：罚方法或消元法
数值稳定性分析： - CFL条件：隐式格式无条件稳定 - 质量矩阵：确保质量守恒 - 能量守恒：验证总能量平衡
收敛性检验： - 网格收敛性：比较不同网格密度的结果 - 时间步长收敛性：验证时间积分精度 - 解析解对比：在特殊情况下与解析解比较
4.3.5 完整实现代码
# 有限元法求解一维热传导问题
# 作者：数值地球教科书项目
# 日期：2025年
library(Matrix)
# 加载字体配置
source("Functions/font_setup.R")
# 设置中文字体（在绘图前）
tryCatch({
setup_bookdown_fonts()
}, error = function(e) {
cat("字体设置失败，使用默认字体：", e$message, "\n")
})## 字体设置完成： Arial Unicode MS## [1] "Arial Unicode MS"# 参数设置
L <- 1.0 # 杆长 (m)
n_elements <- 20 # 单元数量
n_nodes <- n_elements + 1 # 节点数量
h <- L / n_elements # 单元长度
alpha <- 1e-4 # 热扩散系数 (m²/s)
dt <- 100 # 时间步长 (s)
T_max <- 10000 # 总时间 (s)
n_steps <- T_max / dt # 时间步数
# 初始条件
u_ic <- rep(25, n_nodes) # 初始温度 (°C)
u_0 <- 100 # 左端温度 (°C)
u_L <- 50 # 右端温度 (°C)
# 构造质量矩阵和刚度矩阵
M <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
K <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
# 组装全局矩阵
for (i in 1:n_elements) {
# 局部节点编号
local_nodes <- c(i, i+1)
# 局部质量矩阵
M_local <- h/6 * matrix(c(2, 1, 1, 2), 2, 2)
# 局部刚度矩阵
K_local <- alpha/h * matrix(c(1, -1, -1, 1), 2, 2)
# 组装到全局矩阵
M[local_nodes, local_nodes] <- M[local_nodes, local_nodes] + M_local
K[local_nodes, local_nodes] <- K[local_nodes, local_nodes] + K_local
}
# 时间推进（隐式欧拉法）
u_current <- u_ic
u_history <- matrix(0, n_nodes, n_steps + 1)
u_history[, 1] <- u_current
# 构造线性方程组系数矩阵
A <- M + dt * K
b <- M %*% u_current
# 施加边界条件
A[1, ] <- 0; A[1, 1] <- 1
A[n_nodes, ] <- 0; A[n_nodes, n_nodes] <- 1
for (step in 1:n_steps) {
# 更新右端项
b <- M %*% u_current
# 施加边界条件
b[1] <- u_0
b[n_nodes] <- u_L
# 求解线性方程组
u_new <- solve(A, b)
# 更新解
u_current <- as.vector(u_new)
u_history[, step + 1] <- u_current
}
# 结果可视化
x_coords <- seq(0, L, length.out = n_nodes)
time_coords <- seq(0, T_max, by = dt)
# 绘制不同时刻的温度分布
par(mfrow = c(1, 1))
# 选择要显示的时刻
time_indices <- c(1, round(n_steps/16), round(n_steps/8), n_steps + 1)
time_labels <- paste('Time =', time_indices, 's')
colors <- c("blue", "red", "green", "purple")
# 绘制第一条曲线（初始时刻）
plot(x_coords, u_history[, time_indices[1]], type = 'l',
col = colors[1], lwd = 2, ylim = c(0, 100),
xlab = '位置 (m)', ylab = '温度 (°C)',
main = '不同时刻的温度分布')
# 添加其他时刻的曲线
for (i in 2:length(time_indices)) {
lines(x_coords, u_history[, time_indices[i]],
col = colors[i], lwd = 2)
}
# 添加图例
legend("topright", legend = time_labels,
col = colors, lwd = 2, lty = 1, cex = 0.8)
grid()
# 绘制温度随时间变化图（选择几个关键位置）
par(mfrow = c(1, 1))
# 选择要显示的位置
position_indices <- c(1, round(n_nodes/4), round(n_nodes/2), round(3*n_nodes/4), n_nodes)
position_labels <- c("左端", "1/4处", "中点", "3/4处", "右端")
colors_time <- c("blue", "red", "green", "orange", "purple")
# 绘制第一条曲线（左端）
plot(time_coords, u_history[position_indices[1], ], type = 'l',
col = colors_time[1], lwd = 2, ylim = c(0, 100),
xlab = '时间 (s)', ylab = '温度 (°C)',
main = '不同位置的温度随时间变化')
# 添加其他位置的曲线
for (i in 2:length(position_indices)) {
lines(time_coords, u_history[position_indices[i], ],
col = colors_time[i], lwd = 2)
}
# 添加图例
legend("topright", legend = position_labels,
col = colors_time, lwd = 2, lty = 1, cex = 0.8)
grid()
# 返回结果
result <- list(
x_coords = x_coords,
time_coords = time_coords,
u_history = u_history,
parameters = list(
L = L,
n_elements = n_elements,
alpha = alpha,
dt = dt,
T_max = T_max
)
)
# 打印结果摘要
cat("有限元法求解完成！\n")## 有限元法求解完成！## 空间网格点数： 21## 时间步数： 100cat("最终温度范围：", round(min(u_history[, n_steps + 1]), 2), "°C 到",
round(max(u_history[, n_steps + 1]), 2), "°C\n")## 最终温度范围： 50 °C 到 100 °C4.3.6 运行结果示例
运行上述代码可以得到一维热传导问题的有限元解。结果显示：
- 温度分布演化：从初始的均匀温度25°C逐渐向两端边界温度过渡
- 收敛性：解在时间推进过程中逐渐收敛到稳态解
- 边界效应：温度梯度在边界附近较大，在中心区域较小
4.3.7 高级有限元技术实现
除了基本的线性有限元法，我们还实现了更高级的技术：
# 高级有限元法实现
# 包括高阶形函数、自适应网格和与有限差分法的对比
# 作者：数值地球教科书项目
# 日期：2025年
library(Matrix)
# 加载字体配置
source("Functions/font_setup.R")
# 设置中文字体（在绘图前）
tryCatch({
setup_bookdown_fonts()
}, error = function(e) {
cat("字体设置失败，使用默认字体：", e$message, "\n")
})## 字体设置完成： Arial Unicode MS## [1] "Arial Unicode MS"# ============================================================================
# 高阶形函数有限元法
# ============================================================================
# 二次形函数
quadratic_shape_functions <- function(xi) {
N1 <- xi * (xi - 1) / 2
N2 <- 1 - xi^2
N3 <- xi * (xi + 1) / 2
return(c(N1, N2, N3))
}
# 二次形函数的导数
quadratic_shape_derivatives <- function(xi) {
dN1 <- (2*xi - 1) / 2
dN2 <- -2*xi
dN3 <- (2*xi + 1) / 2
return(c(dN1, dN2, dN3))
}
# 高阶有限元求解器
FEM_quadratic <- function(L = 1.0, n_elements = 10, alpha = 1e-4,
dt = 100, T_max = 10000, u_0 = 100, u_L = 50, u_ic = 25) {
n_nodes <- 2 * n_elements + 1 # 二次单元：每个单元3个节点
h <- L / n_elements
# 构造质量矩阵和刚度矩阵
M <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
K <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
# 高斯积分点和权重
gauss_points <- c(-sqrt(3/5), 0, sqrt(3/5))
gauss_weights <- c(5/9, 8/9, 5/9)
# 组装全局矩阵
for (i in 1:n_elements) {
# 局部节点编号（二次单元）
local_nodes <- c(2*i-1, 2*i, 2*i+1)
# 初始化局部矩阵
M_local <- matrix(0, 3, 3)
K_local <- matrix(0, 3, 3)
# 高斯积分
for (gp in 1:3) {
xi <- gauss_points[gp]
weight <- gauss_weights[gp]
# 形函数和导数
N <- quadratic_shape_functions(xi)
dN <- quadratic_shape_derivatives(xi)
# 雅可比变换
J <- h/2
# 组装局部矩阵
for (a in 1:3) {
for (b in 1:3) {
M_local[a, b] <- M_local[a, b] + weight * N[a] * N[b] * J
K_local[a, b] <- K_local[a, b] + weight * dN[a] * dN[b] * J
}
}
}
# 组装到全局矩阵
M[local_nodes, local_nodes] <- M[local_nodes, local_nodes] + M_local
K[local_nodes, local_nodes] <- K[local_nodes, local_nodes] + alpha * K_local
}
# 时间推进
n_steps <- T_max / dt
u_current <- rep(u_ic, n_nodes)
u_history <- matrix(0, n_nodes, n_steps + 1)
u_history[, 1] <- u_current
# 构造线性方程组系数矩阵
A <- M + dt * K
# 施加边界条件
A[1, ] <- 0; A[1, 1] <- 1
A[n_nodes, ] <- 0; A[n_nodes, n_nodes] <- 1
for (step in 1:n_steps) {
b <- M %*% u_current
b[1] <- u_0
b[n_nodes] <- u_L
u_new <- solve(A, b)
u_current <- as.vector(u_new)
u_history[, step + 1] <- u_current
}
return(list(
u_history = u_history,
x_coords = seq(0, L, length.out = n_nodes),
parameters = list(n_elements = n_elements, n_nodes = n_nodes, h = h)
))
}
# ============================================================================
# 自适应网格有限元法
# ============================================================================
# 误差估计器
error_estimator <- function(u, h) {
# 基于梯度的简单误差估计
n <- length(u)
error <- numeric(n-1)
for (i in 1:(n-1)) {
gradient <- abs(u[i+1] - u[i]) / h
error[i] <- h * gradient^2
}
return(error)
}
# 自适应网格细化
adaptive_mesh_refinement <- function(u, x, error_threshold = 0.01) {
n <- length(u)
error <- error_estimator(u, diff(x))
# 标记需要细化的单元
refine_elements <- which(error > error_threshold)
if (length(refine_elements) == 0) {
return(list(x = x, u = u, refined = FALSE))
}
# 在误差大的单元中点插入新节点
new_x <- x
new_u <- u
for (i in sort(refine_elements, decreasing = TRUE)) {
# 插入中点
mid_x <- (x[i] + x[i+1]) / 2
mid_u <- (u[i] + u[i+1]) / 2
new_x <- c(new_x[1:i], mid_x, new_x[(i+1):length(new_x)])
new_u <- c(new_u[1:i], mid_u, new_u[(i+1):length(new_u)])
}
return(list(x = new_x, u = new_u, refined = TRUE))
}
# 自适应有限元求解器
FEM_adaptive <- function(L = 1.0, initial_elements = 10, alpha = 1e-4,
dt = 100, T_max = 10000, u_0 = 100, u_L = 50, u_ic = 25,
max_refinements = 3) {
# 初始网格
n_elements <- initial_elements
n_nodes <- n_elements + 1
h <- L / n_elements
# 初始条件
u_current <- rep(u_ic, n_nodes)
x_coords <- seq(0, L, length.out = n_nodes)
# 自适应细化循环
for (refinement in 1:max_refinements) {
cat("第", refinement, "次网格细化，节点数：", length(x_coords), "\n")
# 使用当前网格求解
result <- FEM_linear_adaptive(x_coords, u_current, alpha, dt, T_max, u_0, u_L)
# 误差估计和网格细化
mesh_result <- adaptive_mesh_refinement(result$u_final, x_coords)
if (!mesh_result$refined) {
cat("网格细化完成\n")
break
}
# 更新网格和解
x_coords <- mesh_result$x
u_current <- mesh_result$u
}
return(result)
}
# 线性有限元求解器（用于自适应网格）
FEM_linear_adaptive <- function(x_coords, u_ic, alpha, dt, T_max, u_0, u_L) {
n_nodes <- length(x_coords)
n_elements <- n_nodes - 1
# 构造质量矩阵和刚度矩阵
M <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
K <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
# 组装全局矩阵
for (i in 1:n_elements) {
h <- x_coords[i+1] - x_coords[i]
local_nodes <- c(i, i+1)
M_local <- h/6 * matrix(c(2, 1, 1, 2), 2, 2)
K_local <- alpha/h * matrix(c(1, -1, -1, 1), 2, 2)
M[local_nodes, local_nodes] <- M[local_nodes, local_nodes] + M_local
K[local_nodes, local_nodes] <- K[local_nodes, local_nodes] + K_local
}
# 时间推进
n_steps <- T_max / dt
u_current <- u_ic
u_history <- matrix(0, n_nodes, n_steps + 1)
u_history[, 1] <- u_current
A <- M + dt * K
A[1, ] <- 0; A[1, 1] <- 1
A[n_nodes, ] <- 0; A[n_nodes, n_nodes] <- 1
for (step in 1:n_steps) {
b <- M %*% u_current
b[1] <- u_0
b[n_nodes] <- u_L
u_new <- solve(A, b)
u_current <- as.vector(u_new)
u_history[, step + 1] <- u_current
}
return(list(
u_history = u_history,
u_final = u_current,
x_coords = x_coords
))
}
# ============================================================================
# 与有限差分法的对比
# ============================================================================
# 有限差分法求解器（用于对比）
FDM_explicit <- function(L = 1.0, n_points = 21, alpha = 1e-4,
dt = 100, T_max = 10000, u_0 = 100, u_L = 50, u_ic = 25) {
h <- L / (n_points - 1)
n_steps <- T_max / dt
# 检查CFL条件
cfl <- alpha * dt / (h^2)
if (cfl > 0.5) {
cat("警告：CFL条件不满足，CFL =", cfl, "\n")
}
# 初始条件
u_current <- rep(u_ic, n_points)
u_history <- matrix(0, n_points, n_steps + 1)
u_history[, 1] <- u_current
# 显式时间推进
for (step in 1:n_steps) {
u_new <- u_current
# 内部点
for (i in 2:(n_points-1)) {
u_new[i] <- u_current[i] + cfl * (u_current[i+1] - 2*u_current[i] + u_current[i-1])
}
# 边界条件
u_new[1] <- u_0
u_new[n_points] <- u_L
u_current <- u_new
u_history[, step + 1] <- u_current
}
return(list(
u_history = u_history,
x_coords = seq(0, L, length.out = n_points),
cfl = cfl
))
}
# 性能对比函数
performance_comparison <- function() {
cat("有限元法与有限差分法性能对比\n")
cat("=" * 50, "\n")
# 测试参数
L <- 1.0
alpha <- 1e-4
dt <- 100
T_max <- 10000
u_0 <- 100
u_L <- 50
u_ic <- 25
# 测试不同网格密度
grid_sizes <- c(11, 21, 41, 81)
for (n in grid_sizes) {
cat("\n网格点数：", n, "\n")
# 有限差分法
t0 <- Sys.time()
fdm_result <- FDM_explicit(L, n, alpha, dt, T_max, u_0, u_L, u_ic)
t1 <- Sys.time()
fdm_time <- as.numeric(t1 - t0)
# 有限元法（线性）
t0 <- Sys.time()
fem_result <- FEM_linear_adaptive(
seq(0, L, length.out = n), rep(u_ic, n),
alpha, dt, T_max, u_0, u_L
)
t1 <- Sys.time()
fem_time <- as.numeric(t1 - t0)
cat("FDM时间：", round(fdm_time, 4), "秒\n")
cat("FEM时间：", round(fem_time, 4), "秒\n")
cat("时间比：", round(fem_time/fdm_time, 2), "\n")
# 精度对比（使用解析解作为参考）
# 这里简化处理，实际应该与解析解比较
fdm_error <- mean(abs(diff(fdm_result$u_history[, ncol(fdm_result$u_history)])))
fem_error <- mean(abs(diff(fem_result$u_final)))
cat("FDM误差：", round(fdm_error, 6), "\n")
cat("FEM误差：", round(fem_error, 6), "\n")
}
}
# ============================================================================
# 主函数和示例
# ============================================================================
# 运行示例
if (FALSE) {
# 基本有限元法
cat("运行基本有限元法...\n")
result_basic <- FEM_quadratic()
# 自适应网格有限元法
cat("运行自适应网格有限元法...\n")
result_adaptive <- FEM_adaptive()
# 性能对比
performance_comparison()
# 可视化结果
par(mfrow = c(1, 2))
# 基本FEM结果 - 不同时刻的温度分布
time_indices <- c(1, round(ncol(result_basic$u_history)/16),
round(ncol(result_basic$u_history)/8), ncol(result_basic$u_history))
time_labels <- paste('Time =', time_indices, 's')
colors <- c("blue", "red", "green", "purple")
# 绘制第一条曲线
plot(result_basic$x_coords, result_basic$u_history[, time_indices[1]], type = 'l',
col = colors[1], lwd = 2, ylim = c(0, 100),
xlab = '位置 (m)', ylab = '温度 (°C)',
main = '基本FEM - 不同时刻温度分布')
# 添加其他时刻的曲线
for (i in 2:length(time_indices)) {
lines(result_basic$x_coords, result_basic$u_history[, time_indices[i]],
col = colors[i], lwd = 2)
}
# 添加图例
legend("topright", legend = time_labels,
col = colors, lwd = 2, lty = 1, cex = 0.8)
grid()
# 自适应FEM结果 - 不同时刻的温度分布
plot(result_adaptive$x_coords, result_adaptive$u_final, type = 'l', ylim = c(0, 100),
col = 'green', lwd = 2, xlab = '位置 (m)', ylab = '温度 (°C)',
main = '自适应FEM - 最终温度分布')
grid()
# 网格对比图
par(mfrow = c(1, 1))
plot(result_basic$x_coords, rep(1, length(result_basic$x_coords)),
type = 'p', pch = 16, col = 'blue', xlab = '位置 (m)', ylab = '',
main = '网格对比', yaxt = 'n', ylim = c(0.8, 1.4))
points(result_adaptive$x_coords, rep(1.2, length(result_adaptive$x_coords)),
pch = 16, col = 'red')
legend('top', legend = c('基本网格', '自适应网格'),
col = c('blue', 'red'), pch = 16)
grid()
}这些高级技术展示了有限元法的强大功能：
- 高阶形函数：使用二次形函数提高精度，收敛阶从\(O(h^2)\)提升到\(O(h^3)\)
- 自适应网格：根据误差估计自动调整网格密度，在解变化剧烈的区域加密网格
- 性能对比：与有限差分法进行系统比较，展示有限元法的优势
- 数值积分：使用高斯积分法提高计算精度
4.4 地球科学应用案例：地下水流模拟
4.4.2 数学模型
控制方程： 二维承压含水层中的地下水流动由以下偏微分方程描述：
\[\frac{\partial}{\partial x}\left(T_{xx}\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y}\left(T_{yy}\frac{\partial h}{\partial y}\right) = S\frac{\partial h}{\partial t} + Q\]
其中： - \(h(x,y,t)\) 是水头（m） - \(T_{xx}\)、\(T_{yy}\) 是导水系数（m²/s） - \(S\) 是储水系数（无量纲） - \(Q\) 是源汇项（1/s）
边界条件： - Dirichlet边界：\(h = h_0\) 在 \(\Gamma_D\) 上 - Neumann边界：\(T\frac{\partial h}{\partial n} = q\) 在 \(\Gamma_N\) 上
4.4.3 有限元离散化
弱形式： 寻找 \(h \in V\)，使得对所有 \(v \in V_0\)：
\[\int_\Omega T\nabla h \cdot \nabla v \, d\Omega + \int_\Omega S\frac{\partial h}{\partial t}v \, d\Omega = \int_\Omega Qv \, d\Omega + \int_{\Gamma_N} qv \, d\Gamma\]
有限元近似： \[h^h(x,y,t) = \sum_{i=1}^{N} h_i(t) N_i(x,y)\]
其中 \(N_i(x,y)\) 是形函数，\(h_i(t)\) 是节点水头值。
4.4.4 实际应用案例
案例描述： 考虑一个矩形承压含水层，尺寸为1000m × 500m，厚度为50m。含水层参数如下： - 导水系数：\(T = 10^{-3}\) m²/s - 储水系数：\(S = 10^{-4}\) - 初始水头：\(h_0 = 100\) m - 边界条件：左边界 \(h = 120\) m，右边界 \(h = 80\) m
有限元求解： 使用三角形单元对计算域进行离散化，采用线性形函数进行插值。
# 有限元法求解二维地下水流问题
# 作者：数值地球教科书项目
# 日期：2025年
library(Matrix)
library(ggplot2)
# 加载字体配置
source("Functions/font_setup.R")
# 设置中文字体（在绘图前）
tryCatch({
setup_bookdown_fonts()
}, error = function(e) {
cat("字体设置失败，使用默认字体：", e$message, "\n")
})## 字体设置完成： Arial Unicode MS## [1] "Arial Unicode MS"# ============================================================================
# 二维地下水流有限元法求解器
# ============================================================================
# 三角形单元形函数（线性）
triangle_shape_functions <- function(xi, eta) {
N1 <- 1 - xi - eta
N2 <- xi
N3 <- eta
return(c(N1, N2, N3))
}
# 三角形单元形函数导数
triangle_shape_derivatives <- function() {
# 对于线性三角形单元，形函数导数为常数
dN_dxi <- c(-1, 1, 0)
dN_deta <- c(-1, 0, 1)
return(list(dN_dxi = dN_dxi, dN_deta = dN_deta))
}
# 生成三角形网格
generate_triangular_mesh <- function(Lx, Ly, nx, ny) {
# 生成节点坐标
x <- seq(0, Lx, length.out = nx)
y <- seq(0, Ly, length.out = ny)
nodes <- expand.grid(x = x, y = y)
n_nodes <- nrow(nodes)
# 生成三角形单元
elements <- list()
element_count <- 0
for (i in 1:(nx-1)) {
for (j in 1:(ny-1)) {
# 每个矩形分成两个三角形
# 左下角节点索引
bottom_left <- (j-1) * nx + i
bottom_right <- (j-1) * nx + i + 1
top_left <- j * nx + i
top_right <- j * nx + i + 1
# 第一个三角形（左下-右上-左上）
element_count <- element_count + 1
elements[[element_count]] <- c(bottom_left, top_right, top_left)
# 第二个三角形（左下-右下-右上）
element_count <- element_count + 1
elements[[element_count]] <- c(bottom_left, bottom_right, top_right)
}
}
return(list(
nodes = nodes,
elements = elements,
n_nodes = n_nodes,
n_elements = length(elements)
))
}
# 计算单元矩阵
compute_element_matrices <- function(nodes, element, T_xx, T_yy, S) {
# 获取单元节点坐标
node_coords <- nodes[element, ]
# 计算雅可比矩阵
x1 <- node_coords[1, 1]; y1 <- node_coords[1, 2]
x2 <- node_coords[2, 1]; y2 <- node_coords[2, 2]
x3 <- node_coords[3, 1]; y3 <- node_coords[3, 2]
# 雅可比矩阵
J <- matrix(c(x2-x1, x3-x1, y2-y1, y3-y1), 2, 2)
det_J <- det(J)
# 雅可比矩阵的逆
J_inv <- solve(J)
# 形函数导数
dN <- triangle_shape_derivatives()
# 计算形函数在物理坐标中的导数
dN_dx <- J_inv[1, 1] * dN$dN_dxi + J_inv[1, 2] * dN$dN_deta
dN_dy <- J_inv[2, 1] * dN$dN_dxi + J_inv[2, 2] * dN$dN_deta
# 计算单元刚度矩阵
K_elem <- matrix(0, 3, 3)
for (i in 1:3) {
for (j in 1:3) {
K_elem[i, j] <- T_xx * dN_dx[i] * dN_dx[j] + T_yy * dN_dy[i] * dN_dy[j]
}
}
K_elem <- K_elem * abs(det_J) / 2
# 计算单元质量矩阵
M_elem <- matrix(0, 3, 3)
for (i in 1:3) {
for (j in 1:3) {
M_elem[i, j] <- S * (1 + (i==j)) / 12 # 简化的质量矩阵
}
}
M_elem <- M_elem * abs(det_J) / 2
return(list(K = K_elem, M = M_elem))
}
# 主求解函数
solve_groundwater_fem <- function(Lx = 1000, Ly = 500, nx = 21, ny = 11,
T_xx = 1e-3, T_yy = 1e-3, S = 1e-4,
h_left = 120, h_right = 80, h_initial = 100,
dt = 1000, T_max = 10000) {
cat("开始有限元法求解二维地下水流问题...\n")
# 生成网格
mesh <- generate_triangular_mesh(Lx, Ly, nx, ny)
nodes <- mesh$nodes
elements <- mesh$elements
n_nodes <- mesh$n_nodes
n_elements <- mesh$n_elements
cat("网格生成完成：", n_nodes, "个节点，", n_elements, "个单元\n")
# 初始化全局矩阵
K_global <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
M_global <- Matrix(0, n_nodes, n_nodes, sparse = TRUE)
# 组装全局矩阵
for (e in 1:n_elements) {
element <- elements[[e]]
elem_matrices <- compute_element_matrices(nodes, element, T_xx, T_yy, S)
# 组装到全局矩阵
K_global[element, element] <- K_global[element, element] + elem_matrices$K
M_global[element, element] <- M_global[element, element] + elem_matrices$M
}
# 初始条件
h_current <- rep(h_initial, n_nodes)
# 时间推进
n_steps <- T_max / dt
h_history <- matrix(0, n_nodes, n_steps + 1)
h_history[, 1] <- h_current
# 构造线性方程组系数矩阵
A <- M_global + dt * K_global
# 施加边界条件
# 左边界（x = 0）
left_boundary <- which(nodes$x == 0)
A[left_boundary, ] <- 0
for (i in left_boundary) {
A[i, i] <- 1
}
# 右边界（x = Lx）
right_boundary <- which(nodes$x == Lx)
A[right_boundary, ] <- 0
for (i in right_boundary) {
A[i, i] <- 1
}
cat("开始时间推进...\n")
for (step in 1:n_steps) {
# 更新右端项
b <- M_global %*% h_current
# 施加边界条件
b[left_boundary] <- h_left
b[right_boundary] <- h_right
# 求解线性方程组
h_new <- solve(A, b)
h_current <- as.vector(h_new)
h_history[, step + 1] <- h_current
if (step %% 10 == 0) {
cat("时间步", step, "/", n_steps, "完成\n")
}
}
cat("求解完成！\n")
return(list(
nodes = nodes,
elements = elements,
h_history = h_history,
h_final = h_current,
parameters = list(
Lx = Lx, Ly = Ly, nx = nx, ny = ny,
T_xx = T_xx, T_yy = T_yy, S = S,
dt = dt, T_max = T_max
)
))
}
# 可视化函数
plot_groundwater_results <- function(result, time_step = NULL) {
nodes <- result$nodes
h_data <- result$h_final
if (!is.null(time_step)) {
h_data <- result$h_history[, time_step]
}
# 创建数据框
plot_data <- data.frame(
x = nodes$x,
y = nodes$y,
h = h_data
)
# 绘制水头分布
p1 <- ggplot(plot_data, aes(x = x, y = y, fill = h)) +
geom_tile() +
scale_fill_viridis_c(name = "水头 (m)") +
labs(title = "地下水流有限元解 - 水头分布",
x = "距离 (m)", y = "距离 (m)") +
theme_minimal() +
theme(text = element_text(family = "SimHei"))
print(p1)
# 绘制等水头线
p2 <- ggplot(plot_data, aes(x = x, y = y, z = h)) +
geom_contour(aes(color = ..level..), bins = 20) +
scale_color_viridis_c(name = "水头 (m)") +
labs(title = "地下水流有限元解 - 等水头线",
x = "距离 (m)", y = "距离 (m)") +
theme_minimal() +
theme(text = element_text(family = "SimHei"))
print(p2)
return(list(contour_plot = p1, level_plot = p2))
}
# 主函数
main_groundwater_fem <- function() {
cat("=== 二维地下水流有限元法求解 ===\n")
# 求解问题
result <- solve_groundwater_fem()
# 可视化结果
plots <- plot_groundwater_results(result)
# 打印结果摘要
cat("\n=== 求解结果摘要 ===\n")
cat("最终水头范围：", round(min(result$h_final), 2), "m 到",
round(max(result$h_final), 2), "m\n")
cat("水头梯度：", round((max(result$h_final) - min(result$h_final)) / 1000, 4), "m/m\n")
return(result)
}
# 运行示例
if (FALSE) {
# 运行主函数
result <- main_groundwater_fem()
# 绘制不同时刻的水头分布
par(mfrow = c(2, 2))
time_steps <- c(1, round(ncol(result$h_history)/4),
round(ncol(result$h_history)/2), ncol(result$h_history))
for (i in 1:4) {
plot_groundwater_results(result, time_steps[i])
}
}结果分析： 1. 水头分布：从左侧高水头向右侧低水头流动 2. 流线模式：流线垂直于等水头线，符合达西定律 3. 时间演化：水头分布随时间逐渐达到稳态
4.5 有限元法的收敛性分析与误差估计
4.5.1 收敛性理论基础
有限元法的收敛性是数值分析中的核心问题，它回答了”随着网格加密，数值解是否趋向于精确解”这一基本问题。
4.5.2 先验误差估计
4.5.2.1 Céa引理
Céa引理是有限元误差分析的基石：
\[\|u - u^h\|_V \leq \frac{M}{\alpha} \inf_{v^h \in V^h} \|u - v^h\|_V\]
其中： - \(M\) 是双线性形式 \(a(\cdot,\cdot)\) 的连续性常数 - \(\alpha\) 是强制性常数 - \(V^h\) 是有限元空间
这表明有限元误差受最佳逼近误差控制。
4.5.3 后验误差估计
后验误差估计基于已计算的数值解，提供误差的实时评估。
4.5.3.1 残差型误差估计
定义残差： \[R(u^h) = f + \nabla \cdot (k \nabla u^h) \quad \text{in } \Omega_e\] \[r(u^h) = \llbracket k \nabla u^h \cdot \mathbf{n} \rrbracket \quad \text{on } \partial\Omega_e\]
其中 \(\llbracket \cdot \rrbracket\) 表示跳跃。
误差估计子： \[\eta_e^2 = h_e^2 \|R(u^h)\|_{L^2(\Omega_e)}^2 + h_e \|r(u^h)\|_{L^2(\partial\Omega_e)}^2\]
总误差估计： \[\|u - u^h\|_{H^1} \leq C \left(\sum_{e} \eta_e^2\right)^{1/2}\]
4.5.4 自适应网格细化
基于后验误差估计的自适应策略：
4.5.5 数值实验与验证
4.5.5.1 制造解方法
为验证收敛阶，使用制造解： \[u_{\text{exact}}(x,t) = \sin(\pi x) e^{-\pi^2 \alpha t}\]
对应的源项： \[f = \frac{\partial u_{\text{exact}}}{\partial t} - \alpha \frac{\partial^2 u_{\text{exact}}}{\partial x^2}\]
4.6 有限元法与有限差分法的系统对比
4.6.2 精度和稳定性对比
精度特征： - 有限元法：具有最优收敛阶，\(L^2\) 范数超收敛 - 有限差分法：收敛阶取决于差分格式，通常较低
稳定性特征： - 有限元法：隐式格式无条件稳定，质量矩阵确保稳定性 - 有限差分法：显式格式有CFL限制，隐式格式稳定性较好
4.7 高级有限元技术的深入探讨
4.7.1 高阶形函数技术
高阶形函数是提高有限元精度的重要手段，特别适用于解具有高正则性的问题。
4.7.2 自适应有限元技术
4.7.2.1 \(h\)-型自适应
网格细化策略： 1. 规则细化：将单元等分为子单元 2. 不规则细化：只细化部分边或面 3. 各向异性细化：根据解的各向异性特征细化
挂节点处理： 在不规则网格中，挂节点的约束关系： \[u_h = \frac{1}{2}(u_i + u_j)\]
其中\(u_i\)、\(u_j\)是相邻节点的值。
4.7.3 稳定化技术
4.7.3.1 流线迁移扩散方程
对于对流占优问题： \[-\epsilon \Delta u + \mathbf{b} \cdot \nabla u + cu = f\]
当\(\epsilon \ll |\mathbf{b}|\)时，标准Galerkin法会产生数值振荡。
4.7.3.2 SUPG稳定化
流线迁移Petrov-Galerkin方法： 修改的弱形式： \[a(u^h, v^h) + \sum_e \tau_e \int_{\Omega_e} \mathcal{L}u^h \cdot (\mathbf{b} \cdot \nabla v^h) d\Omega = L(v^h) + \sum_e \tau_e \int_{\Omega_e} f \cdot (\mathbf{b} \cdot \nabla v^h) d\Omega\]
其中稳定化参数： \[\tau_e = \frac{h_e}{2|\mathbf{b}|} \coth(Pe_e) - \frac{1}{2|\mathbf{b}|}\]
\(Pe_e = \frac{|\mathbf{b}|h_e}{2\epsilon}\)是单元Peclet数。
4.7.4 多尺度有限元方法
4.7.5 间断有限元方法
4.7.8 高级技术的应用实例
4.8 总结
有限元法作为一种强大的数值分析方法，在地球科学中具有广泛的应用前景。通过本章的学习，读者应该掌握：
- 理论基础：理解变分原理和弱形式的概念
- 数学推导：掌握形函数和单元矩阵的构造方法
- 编程实现：能够使用R语言实现简单的有限元程序
- 方法对比：了解有限元法与有限差分法的优缺点
- 应用前景：认识有限元法在地球科学中的潜在应用
有限元法的优势在于其灵活性和精度，特别适合处理复杂几何和多物理场耦合问题。随着计算技术的发展，有限元法将在地球科学数值模拟中发挥越来越重要的作用。