2. 稳态热传导与热应力顺序耦合

任务目标

实现一个二维稳态热传导与热应力的顺序耦合求解器：首先求解给定边界条件下的温度场，然后将温度场作为体积载荷，求解结构在约束下的热应力与热变形。本案例采用有限差分法离散温度场，采用有限元法（常应变三角形单元）求解热应力，形成完整的顺序耦合闭环。

物理模型与耦合流程

考虑一个边长为 L=1.0 m 的方形薄板，材料为各向同性线弹性体。温度场满足稳态热传导方程 ∇²T = 0，边界条件为：左边界 T=100°C，右边界 T=0°C，上下边界绝热。求得温度场后，将温度变化 ΔT = T(x,y) - T_ref（参考温度 T_ref=0°C）作为热载荷，在结构左边界施加固定约束（u_x = u_y = 0），求解热应力场。

顺序耦合的关键在于：温度场求解独立于应力场，但应力场依赖温度场的空间分布。本实现中温度场采用有限差分法在规则网格上求解，应力场采用同一套网格生成的三角形单元进行有限元计算，保证数据传递的一致性。

环境依赖

仅需 Python 3.8+ 及 NumPy，无需其他第三方库。代码自包含全部求解逻辑。

关键完整代码


import numpy as np

# ========== 参数定义 ==========
L = 1.0              # 板边长 (m)
nx, ny = 20, 20      # 网格数（温度场与应力场共用节点）
E = 200e9            # 弹性模量 (Pa)
nu = 0.3             # 泊松比
alpha = 1.2e-5       # 热膨胀系数 (1/°C)
T_ref = 0.0          # 参考温度 (°C)
T_left, T_right = 100.0, 0.0   # 左右边界温度

# ========== 1. 稳态温度场求解（有限差分 + 高斯-赛德尔迭代） ==========
dx = L / nx
dy = L / ny
T = np.zeros((ny + 1, nx + 1))
T[:, 0] = T_left      # 左边界
T[:, -1] = T_right    # 右边界

# 迭代求解拉普拉斯方程（绝热上下边界：Neumann条件）
for _ in range(20000):
    T_old = T.copy()
    for j in range(1, ny):
        for i in range(1, nx):
            T[j, i] = 0.25 * (T_old[j, i+1] + T_old[j, i-1] +
                              T_old[j+1, i] + T_old[j-1, i])
    # 上下边界绝热（∂T/∂y = 0）
    T[0, :] = T[1, :]
    T[-1, :] = T[-2, :]
    if np.max(np.abs(T - T_old)) < 1e-6:
        break

# ========== 2. 热应力有限元求解（常应变三角形单元） ==========
# 生成节点坐标（与温度场网格一致）
n_nodes = (nx + 1) * (ny + 1)
nodes = np.zeros((n_nodes, 2))
for j in range(ny + 1):
    for i in range(nx + 1):
        nodes[j * (nx + 1) + i] = [i * dx, j * dy]

# 生成三角形单元（每个网格分成2个三角形）
elements = []
for j in range(ny):
    for i in range(nx):
        n0 = j * (nx + 1) + i
        n1 = n0 + 1
        n2 = n0 + (nx + 1)
        n3 = n2 + 1
        elements.append([n0, n1, n2])   # 左下三角形
        elements.append([n1, n3, n2])   # 右上三角形
elements = np.array(elements)
n_elem = len(elements)

# 组装刚度矩阵与热载荷向量
K = np.zeros((2 * n_nodes, 2 * n_nodes))
F = np.zeros(2 * n_nodes)

# 材料矩阵（平面应力）
D = E / (1 - nu**2) * np.array([
    [1, nu, 0],
    [nu, 1, 0],
    [0, 0, (1 - nu) / 2]
])

for e in range(n_elem):
    idx = elements[e]
    coords = nodes[idx]   # 3x2
    xe, ye = coords[:, 0], coords[:, 1]
    
    # 形函数导数（常应变三角形）
    B = np.zeros((3, 6))
    area = 0.5 * ((xe[1] - xe[0]) * (ye[2] - ye[0]) -
                  (xe[2] - xe[0]) * (ye[1] - ye[0]))
    if abs(area) < 1e-12:
        continue
    dN_dx = np.array([ye[1] - ye[2], ye[2] - ye[0], ye[0] - ye[1]]) / (2 * area)
    dN_dy = np.array([xe[2] - xe[1], xe[0] - xe[2], xe[1] - xe[0]]) / (2 * area)
    B[0, 0::2] = dN_dx
    B[1, 1::2] = dN_dy
    B[2, 0::2] = dN_dy
    B[2, 1::2] = dN_dx
    
    # 单元刚度矩阵
    Ke = B.T @ D @ B * area
    dofs = np.zeros(6, dtype=int)
    for k in range(3):
        dofs[2*k] = 2 * idx[k]
        dofs[2*k+1] = 2 * idx[k] + 1
    K[np.ix_(dofs, dofs)] += Ke
    
    # 热载荷（等效节点力）
    T_elem = np.mean(T[idx[:, 1], idx[:, 0]])   # 单元平均温度
    dT = T_elem - T_ref
    if abs(dT) > 1e-12:
        # 热应变 ε_th = alpha * dT * [1, 1, 0]^T
        eth = alpha * dT * np.array([1.0, 1.0, 0.0])
        sigma_th = D @ eth
        Fe = B.T @ sigma_th * area
        F[dofs] -= Fe

# 施加约束：左边界节点 (x=0) 固定 ux=uy=0
constrained_dofs = []
for j in range(ny + 1):
    node_id = j * (nx + 1)   # x=0 的节点
    constrained_dofs.extend([2 * node_id, 2 * node_id + 1])
free_dofs = np.setdiff1d(np.arange(2 * n_nodes), constrained_dofs)

# 求解位移
K_ff = K[np.ix_(free_dofs, free_dofs)]
F_f = F[free_dofs]
u_free = np.linalg.solve(K_ff, F_f)

u = np.zeros(2 * n_nodes)
u[free_dofs] = u_free

# ========== 3. 后处理：计算应力与最大位移 ==========
# 提取每个单元的平均应力
elem_stress = np.zeros((n_elem, 3))
for e in range(n_elem):
    idx = elements[e]
    coords = nodes[idx]
    xe, ye = coords[:, 0], coords[:, 1]
    area = 0.5 * ((xe[1] - xe[0]) * (ye[2] - ye[0]) -
                  (xe[2] - xe[0]) * (ye[1] - ye[0]))
    dN_dx = np.array([ye[1] - ye[2], ye[2] - ye[0], ye[0] - ye[1]]) / (2 * area)
    dN_dy = np.array([xe[2] - xe[1], xe[0] - xe[2], xe[1] - xe[0]]) / (2 * area)
    B = np.zeros((3, 6))
    B[0, 0::2] = dN_dx
    B[1, 1::2] = dN_dy
    B[2, 0::2] = dN_dy
    B[2, 1::2] = dN_dx
    dofs = np.zeros(6, dtype=int)
    for k in range(3):
        dofs[2*k] = 2 * idx[k]
        dofs[2*k+1] = 2 * idx[k] + 1
    u_e = u[dofs]
    strain = B @ u_e
    T_elem = np.mean(T[idx[:, 1], idx[:, 0]])
    dT = T_elem - T_ref
    eth = alpha * dT * np.array([1.0, 1.0, 0.0])
    elem_stress[e] = D @ (strain - eth)

# 输出关键结果
ux = u[0::2].reshape(ny + 1, nx + 1)
uy = u[1::2].reshape(ny + 1, nx + 1)
max_disp = np.sqrt(ux**2 + uy**2).max()
max_von_mises = np.sqrt(elem_stress[:, 0]**2 - elem_stress[:, 0]*elem_stress[:, 1] +
                        elem_stress[:, 1]**2 + 3*elem_stress[:, 2]**2).max()

print(f"温度场迭代收敛步数: {_+1}")
print(f"中心点温度 T({L/2},{L/2}) = {T[ny//2, nx//2]:.4f} °C")
print(f"最大热变形: {max_disp:.6e} m")
print(f"最大von Mises应力: {max_von_mises:.6e} Pa")

运行与验证

运行方式：将上述代码保存为 thermal_stress_coupling.py，在命令行执行 python thermal_stress_coupling.py。程序将输出温度场迭代步数、中心点温度、最大变形与最大von Mises应力。

验证方法：

  1. 温度场正确性：由于左右边界温度差为100°C且上下绝热，温度场应为线性分布 T(x) = 100 - 100*x/L。中心点 (0.5, 0.5) 的理论温度为 50°C。程序输出应接近50°C（误差<0.5°C）。
  2. 应力场合理性：在自由膨胀但左边界固定的条件下，板整体受热膨胀。由于温度沿x方向线性递减，板内产生压缩热应力（左侧温度高，膨胀趋势大但被约束）。最大von Mises应力应出现在左边界附近，量级约为 E*alpha*ΔT_max ≈ 200e9 * 1.2e-5 * 100 = 2.4e8 Pa（考虑约束效应后略低）。
  3. 位移趋势：最大位移应出现在右边界（温度最低、自由膨胀最大处），方向为x正方向，量级约为 alpha * ΔT * L ≈ 1.2e-5 * 100 * 1.0 = 1.2e-3 m。

若输出结果与上述理论估算在同一量级且趋势一致，则顺序耦合实现正确。

常见问题

  • 温度场不收敛：检查迭代次数上限是否足够（当前20000次），或收敛容差是否过严（1e-6已足够）。
  • 应力求解奇异：确认约束已正确施加（左边界所有节点两个自由度均固定），且刚度矩阵未出现零主元。
  • 网格依赖性：增大 nx, ny 可提高精度，但计算时间增加。当前20×20网格在普通PC上运行时间约数秒。
← 上一章 📖 返回目录 下一章 →