3. 一维有限元实现：线性单元形函数，单元刚度矩阵组装，Dirichlet边界处理

好，咱们今天来点真格的。前面聊完了加权余量法和变分原理，算是把理论基础铺好了。现在，我要带你亲手走一遍一维有限元的完整流程。从形函数开始，到组装总刚，再到处理边界条件——每一步我都会把项目里踩过的坑指给你看。

3.1 线性单元形函数：为什么是它？

形函数，说白了就是单元内部的“插值规则”。它决定了我们怎么用节点位移算出单元内任意一点的位移。

对于一维线性单元，只有两个节点。形函数必须满足一个基本条件：在本节点处值为1，在另一个节点处值为0。你想想看，这样才能保证节点位移就是它自己的值。

我习惯用自然坐标系来推导。假设单元长度是 L，局部坐标 ξ 从 -1 到 1：

N₁(ξ) = (1 - ξ) / 2
N₂(ξ) = (1 + ξ) / 2

就这么简单。但简单归简单，我见过不少新手在这里翻车——把形函数和插值函数搞混。形函数是定义在单元上的，插值函数是定义在整个域上的。这是两个概念。

3.2 单元刚度矩阵推导

有了形函数，单元刚度矩阵就是水到渠成的事。以一维热传导问题为例，控制方程是：

-k * d²u/dx² = f

单元刚度矩阵的通用表达式是：

K_e = ∫ Bᵀ * D * B dΩ

其中 B 是形函数导数矩阵，D 是材料本构矩阵。对于一维问题：

B = [dN₁/dx, dN₂/dx]
D = k (导热系数)

用链式法则把局部坐标变换到物理坐标：

dNᵢ/dx = dNᵢ/dξ * dξ/dx
dξ/dx = 2/L

所以：

B = [-1/L, 1/L]

代入积分：

K_e = ∫₋₁¹ Bᵀ * k * B * (L/2) dξ
     = (k/L) * [[1, -1], [-1, 1]]

你看，结果就是一个2×2的对称矩阵。我在项目中做过上千次这样的推导，每次看到这个矩阵都觉得特别踏实——它完美体现了“近刚远柔”的物理直觉。

参数	符号	说明
单元长度	L	两个节点间的距离
导热系数	k	材料属性，各向同性
形函数导数矩阵	B	[-1/L, 1/L]
单元刚度矩阵	K_e	(k/L) * [[1, -1], [-1, 1]]

3.3 全局组装：从单元到系统

组装，就是把所有单元的贡献“对号入座”到全局矩阵里。每个单元有两个节点，每个节点对应一个全局自由度编号。

我习惯用一个“自由度映射表”来管理这个过程。比如单元连接两个全局节点 i 和 j：

K_global[i][i] += K_e[0][0]
K_global[i][j] += K_e[0][1]
K_global[j][i] += K_e[1][0]
K_global[j][j] += K_e[1][1]

就这么简单？嗯，原理上确实简单。但实际写代码时，我踩过一个坑：忘记初始化全局矩阵。Python里如果你用列表嵌套列表，不初始化就直接累加，会报错。我后来统一用 numpy.zeros 来初始化，再也没出过问题。

3.4 Dirichlet边界处理：两种主流方法

边界条件处理，是有限元里最容易出bug的地方。Dirichlet边界（强制边界）要求节点位移等于给定值。常用的方法有两种：置大数法和消行法。

3.4.1 置大数法

这个方法很暴力，但很有效。假设第 i 个自由度的位移被指定为 u₀：

K_global[i][i] = 1e30  # 一个大数
F_global[i] = 1e30 * u₀

为什么这样能行？因为方程 K * u = F 中，第 i 行被一个大数主导，其他项可以忽略不计，解出来 u[i] ≈ u₀。

我个人不太喜欢这个方法。虽然实现简单，但大数的选择很玄学——太小了约束不够强，太大了可能引起数值问题。我见过有人用 1e15，有人用 1e30，没有一个统一标准。

3.4.2 消行法（我推荐）

这个方法更干净。直接把已知自由度从方程组里“消掉”：

1. 把 K_global 中第 i 行和第 i 列全部置零
2. 把 K_global[i][i] 设为 1
3. 把 F_global[i] 设为 u₀
4. 把右端项中其他行减去 K_global[j][i] * u₀（如果之前没置零）

我习惯用第二种方式：先备份原始矩阵，然后直接修改。这样代码可读性好，调试也方便。

def apply_dirichlet(K, F, dof_idx, value):
    """消行法处理Dirichlet边界"""
    n = K.shape[0]
    # 备份右端项
    F_copy = F.copy()
    # 修改矩阵和右端项
    for j in range(n):
        if j != dof_idx:
            F_copy[j] -= K[j, dof_idx] * value
        K[dof_idx, j] = 0
        K[j, dof_idx] = 0
    K[dof_idx, dof_idx] = 1
    F_copy[dof_idx] = value
    return K, F_copy

3.5 完整流程与SVG框架图

下面这张图，是我做一维有限元求解时脑子里始终挂着的流程。每一步都有坑，每一步都有技巧。

3.6 完整代码示例

最后，我给你一个完整的Python实现。这段代码我用了好几年，结构清晰，适合教学也适合工程使用。

import numpy as np

def linear_shape_function(xi):
    """一维线性单元形函数"""
    N1 = (1 - xi) / 2
    N2 = (1 + xi) / 2
    return np.array([N1, N2])

def element_stiffness_1d(k, L):
    """一维线性单元刚度矩阵"""
    return (k / L) * np.array([[1, -1], [-1, 1]])

def assemble_global(n_nodes, elements, k_values):
    """全局刚度矩阵组装"""
    K = np.zeros((n_nodes, n_nodes))
    for elem, k in zip(elements, k_values):
        i, j = elem
        Ke = element_stiffness_1d(k, 1.0)  # 假设单元长度归一化
        K[i, i] += Ke[0, 0]
        K[i, j] += Ke[0, 1]
        K[j, i] += Ke[1, 0]
        K[j, j] += Ke[1, 1]
    return K

def apply_dirichlet_elimination(K, F, dof_idx, value):
    """消行法处理Dirichlet边界"""
    n = K.shape[0]
    F_copy = F.copy()
    for j in range(n):
        if j != dof_idx:
            F_copy[j] -= K[j, dof_idx] * value
        K[dof_idx, j] = 0
        K[j, dof_idx] = 0
    K[dof_idx, dof_idx] = 1
    F_copy[dof_idx] = value
    return K, F_copy

# 示例：3节点2单元问题
n_nodes = 3
elements = [(0, 1), (1, 2)]
k_values = [1.0, 1.0]

K = assemble_global(n_nodes, elements, k_values)
F = np.array([0.0, 0.0, 10.0])  # 右端项

# 处理边界：节点0位移为0，节点2位移为5
K, F = apply_dirichlet_elimination(K, F, 0, 0.0)
K, F = apply_dirichlet_elimination(K, F, 2, 5.0)

# 求解
u = np.linalg.solve(K, F)
print("节点位移:", u)