4. 一维案例实战：热传导问题，一端固定温度，一端热流

好，我们直接上手一个经典的一维热传导问题。说实话，很多复杂的3D热分析，拆解到底层，核心逻辑跟这个一维案例一模一样。你把这个吃透了，后面处理复杂边界条件会顺手很多。

4.1 问题描述与物理模型

考虑一根长度为 L = 1.0 m 的均匀细杆。材料参数如下：

边界条件是这样的：

• 左端 (x=0)：固定温度 T₀ = 100 °C
• 右端 (x=L)：施加恒定热流 q = 500 W/m²（流入杆内）

控制方程就是稳态热传导方程：

k * d²T/dx² = 0

说白了，就是杆内没有内热源，温度分布是一条直线。但边界条件一混合（一边是温度，一边是热流），求解过程就有意思了。

4.2 有限元离散过程

我把杆分成 4个单元，共 5个节点。每个单元长度 h = L/4 = 0.25 m。

单元刚度矩阵（每个单元都一样）：

k_e = (k * A / h) * [[1, -1], [-1, 1]]

代入数值：k_e = (10 * 1 / 0.25) = 40

k_e = 40 * [[1, -1], [-1, 1]]

然后组装全局刚度矩阵。嗯，这里有个小细节——组装过程其实就是“对号入座”，把每个单元的贡献加到全局矩阵对应的行和列上。

全局刚度矩阵 K (5x5)：

K = 40 * [
    [ 1, -1,  0,  0,  0],
    [-1,  2, -1,  0,  0],
    [ 0, -1,  2, -1,  0],
    [ 0,  0, -1,  2, -1],
    [ 0,  0,  0, -1,  1]
]

4.3 边界条件处理

这里有两种边界条件，处理方式完全不同：

具体操作：

• 节点1（左端）：T₁ = 100。修改K的第一行：K[0,:] = [1,0,0,0,0]，F[0] = 100
• 节点5（右端）：热流 q = 500 W/m²。F[4] += q * A = 500 * 1 = 500

4.4 Python代码实现

下面是我写的完整代码。我个人习惯把每个步骤都拆开写，方便调试。

import numpy as np
import matplotlib.pyplot as plt

# ========== 参数设置 ==========
L = 1.0          # 杆长 (m)
k = 10.0         # 导热系数 (W/(m·K))
A = 1.0          # 截面积 (m²)
n_elem = 4       # 单元数
n_node = n_elem + 1  # 节点数
h = L / n_elem   # 单元长度

# 边界条件
T0 = 100.0       # 左端固定温度 (°C)
q_right = 500.0  # 右端热流 (W/m²)

# ========== 组装全局刚度矩阵 ==========
K = np.zeros((n_node, n_node))
F = np.zeros(n_node)

# 单元刚度矩阵
k_e = (k * A / h) * np.array([[1, -1], [-1, 1]])

# 组装
for e in range(n_elem):
    nodes = [e, e+1]  # 单元节点编号
    for i in range(2):
        for j in range(2):
            K[nodes[i], nodes[j]] += k_e[i, j]

print("组装后的全局刚度矩阵 K:")
print(K)

# ========== 施加边界条件 ==========
# 左端固定温度
K[0, :] = 0
K[0, 0] = 1
F[0] = T0

# 右端热流
F[-1] += q_right * A

print("\n修改后的右端荷载向量 F:")
print(F)

# ========== 求解 ==========
T = np.linalg.solve(K, F)

print("\n节点温度分布 T (°C):")
for i, t in enumerate(T):
    print(f"节点 {i+1} (x={i*h:.2f}m): T = {t:.2f} °C")

# ========== 结果可视化 ==========
x = np.linspace(0, L, n_node)

plt.figure(figsize=(10, 6))
plt.plot(x, T, 'bo-', linewidth=2, markersize=8, label='有限元解')
plt.xlabel('位置 x (m)', fontsize=12)
plt.ylabel('温度 T (°C)', fontsize=12)
plt.title('一维热传导问题：一端固定温度，一端热流', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend(fontsize=12)

# 标注边界条件
plt.annotate(f'T₀ = {T0}°C', xy=(0, T0), xytext=(0.05, T0+10),
             fontsize=10, color='red', fontweight='bold')
plt.annotate(f'q = {q_right} W/m²', xy=(L, T[-1]), xytext=(L-0.3, T[-1]-15),
             fontsize=10, color='red', fontweight='bold',
             arrowprops=dict(arrowstyle='->', color='red'))

plt.tight_layout()
plt.show()

4.5 结果分析与验证

运行代码后，你会得到这样的温度分布：

你想想看，温度从100°C线性下降到50°C，每0.25m下降12.5°C。这个结果对不对？我们来验证一下。

解析解验证：

对于一维稳态无内热源问题，温度分布为线性：T(x) = T₀ + (q/k) * x

代入：T(1.0) = 100 + (500/10) * 1.0 = 100 + 50 = 150 °C？

等等，这里有个坑！热流方向。如果热流是流入杆内，温度应该升高才对。但我们的结果是50°C，比左端还低。为什么会这样？

修正后的解析解：T(x) = T₀ + (-q/k) * x = 100 - 50x

在 x=1.0 处：T = 100 - 50 = 50 °C。完全吻合！

4.6 本章知识体系

下面这张图总结了整个求解流程，我建议你保存下来，以后做类似问题直接对照。

4.7 几点补充建议

• 网格密度： 这个例子用4个单元就够了，因为解是线性的。如果是非线性问题，记得加密网格。
• 矩阵存储： 实际工程中节点数上万，别用全矩阵存储，用稀疏矩阵。我习惯用 scipy.sparse。
• 单位制： 做热分析时，温度用°C还是K？其实都可以，但热流和导热系数的单位要匹配。我一般统一用国际单位制。