本内容针对西安交通大学材料学院凝固实验 中 温度场数值模拟

网格剖分

本实验如下图所示，采用矩形网格来划分求解域，网格步长分别为$\delta x$及$\delta y$

微分方程离散化

网格生成后，就要将控制方程离散化，即将控制偏微分方程转化为各个节点上的代数方程组。直角坐标系内二维导热过程温度场控制微分方程：

$$ \frac{\partial T}{\partial t}=\alpha\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}\right) $$

$T(x,y,t)$为温度，$t$为时间，$\alpha$为材料热扩散系数
把控制方程中的各阶导数用相应的差分表达式来代替而形成离散方程的方法称为泰勒级数展开法。由导数定义有：

$$ \frac{\mathrm{d} T}{\mathrm{dx}}=\lim _{\Delta x \rightarrow 0} \frac{T(x+\Delta x)-T(x)}{\Delta x}=\lim _{\Delta x \rightarrow 0} \frac{\Delta T}{\Delta x} $$

泰勒级数展开：

$$ T(x+\Delta x)=T(x)+\Delta x \cdot \frac{\mathrm{~d} T}{\mathrm{~d} x}+\frac{(\Delta x)^2}{2!} \cdot \frac{\mathrm{d}^2 T}{\mathrm{~d} x^2}+\cdots+\frac{(\Delta x)^n}{n!} \cdot \frac{\mathrm{d}^{\mathrm{n}} T}{\mathrm{~d} x^n}+\cdots $$

$$ \frac{T(x+\Delta x)-T(x)}{\Delta x}=\frac{\mathrm{d} T}{\mathrm{~d} x}+\frac{\Delta x}{2!} \cdot \frac{\mathrm{d}^2 T}{\mathrm{~d} x^2}+\cdots+\frac{(\Delta x)^{n-1}}{n!} \cdot \frac{\mathrm{d}^{\mathrm{n}} T}{\mathrm{~d} x^n}+\cdots $$

$$ \frac{\mathrm{d} T}{\mathrm{~d} x}=\frac{T(x+\Delta x)-T(x)}{\Delta x}, \mathrm{O}(\Delta x) $$

其中，$O(\delta x)$代替了高阶导数纸盒，称为 截断误差 。另外，上式称为$(dT/dx)$的向前差分，此外还有向后差分：

$$ \frac{\mathrm{d} T}{\mathrm{~d} x}=\frac{T(x)-T(x-\Delta x)}{\Delta x}, \mathrm{O}(\Delta x) $$

以及中心差分：

$$ \frac{\mathrm{d} T}{\mathrm{~d} x}=\frac{T(x+\Delta x)-T(x-\Delta x)}{2 \Delta x}, \mathrm{O}\left(\Delta x^2\right) $$

对于二阶或更高阶导数的差分表达式，如二阶：

$$ \frac{\mathrm{d}^2 T}{\mathrm{~d} x^2}=\frac{T(x+\Delta x)-2 T(x)+T(x-\Delta x)}{\Delta x^2}, \mathrm{O}\left(\Delta x^2\right) $$

由此，对于二元情况，参照上面所示的网格系统

$$ \frac{T_{i, j}^{n+1}-T_{i, j}^n}{\Delta t}=\alpha\left[\frac{T_{i+1, j}^n-2 T_{i, j}^n+T_{i-1, j}^n}{\Delta x^2}+\frac{T_{i, j+1}^n-2 T_{i, j}^n+T_{i, j-1}^n}{\Delta y^2}\right] $$

因此，各个单元温度计算格式为：

$$ \begin{aligned} &T_{i, j}^{n+1}=M_1\left[T_{i-1, j}^n+T_{i+1, j}^n\right]+M_2\left[T_{i, j-1}^n+T_{i, j+1}^n\right]+\left(1-2 M_1-2 M_2\right) T_{i, j}^n\\ &\text { 其中 }\\ &M_1=\frac{\alpha \Delta t}{\Delta x^2}, \quad M_2=\frac{\alpha \Delta t}{\Delta y^2} \end{aligned} $$

其稳定性条件为式中各项系数应大于零 ，故$2 M_1+2 M_2 \leq 1$
即：

$$ \Delta t \leq \frac{1}{2 \alpha\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)} $$

初始条件和边界条件

初始条件

初始条件指所求解问题的初始温度场，也就是在零时刻温度场的分布。它可以是均匀的，也可以是不均匀的。

边界条件

$$ \left.\begin{array}{l} \left.T\right|_{\Gamma}=T_W \\ \left.T\right|_{\Gamma}=T_W(x, y, t) \end{array}\right\} $$

实验内容

Python 代码实现

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# 常量定义
L = 8.0  # 长度 (cm)
T0 = 500.0  # 初始温度 (摄氏度)
Tw = 100.0  # 表面温度 (摄氏度)
alpha = 0.05  # 热扩散系数 (cm^2/s)
dx = 0.1  # 空间步长 (cm)
dy = 0.1  # 空间步长 (cm)
dt = 0.01  # 时间步长 (s)

# 检查稳定性条件

# 网格大小
nx, ny = int(L/dx), int(L/dy)

# 创建初始温度场
T = T0 * np.ones((nx, ny))

# 边界条件设置
def apply_boundary_conditions(T):
    T[0, :] = Tw
    T[-1, :] = Tw
    T[:, 0] = Tw
    T[:, -1] = Tw
    return T

# 矢量化计算一步的温度分布
def step_temperature(T, alpha, dx, dy, dt):
    T_new = T.copy()
    # 使用矢量化计算非边界点的新温度
    T_new[1:-1, 1:-1] = T[1:-1, 1:-1] + alpha * dt * (
        (T[2:, 1:-1] - 2*T[1:-1, 1:-1] + T[:-2, 1:-1]) / dx**2 +
        (T[1:-1, 2:] - 2*T[1:-1, 1:-1] + T[1:-1, :-2]) / dy**2
    )
    return apply_boundary_conditions(T_new)

# 计算整个温度场随时间变化
total_time = 15 * 60  # 15分钟
total_steps = int(total_time / dt)
T_all = [T0 * np.ones((nx, ny))]
apply_boundary_conditions(T_all[0])
T = T_all[0].copy()
for _ in range(total_steps):
    T = step_temperature(T, alpha, dx, dy, dt)
    T_all.append(T.copy())

# 创建图形窗口和布局
fig, axes = plt.subplots(3, 2, figsize=(12, 12))
fig.subplots_adjust(hspace=0.4)

# 解包子图
ax_anim, ax_curve, ax_1min, ax_5min, ax_10min = axes.flatten()[:5]

# 动画模拟设置
cax = ax_anim.imshow(T_all[0], cmap='coolwarm', vmin=Tw, vmax=T0)
ax_anim.set_title("Temperature Animation")

# 中心点冷却曲线
center_temperatures = [T[nx//2, ny//2] for T in T_all]
time_array = np.arange(total_steps+1) * dt
ax_curve.plot(time_array, center_temperatures)
ax_curve.set_title("Temperature at Center vs Time")
ax_curve.set_xlabel("Time (s)")
ax_curve.set_ylabel("Temperature (°C)")

# 特定时间点的温度场
def plot_fixed_time(ax, time_point, title):
    time_step = int(time_point / dt)
    T_display = T_all[time_step]
    im = ax.imshow(T_display, cmap='coolwarm', vmin=Tw, vmax=T0)
    ax.set_title(title)
    ax.set_xlabel('x (cm)')
    ax.set_ylabel('y (cm)')
    fig.colorbar(im, ax=ax)

plot_fixed_time(ax_1min, 60, 'Temperature at 1 min')
plot_fixed_time(ax_5min, 300, 'Temperature at 5 min')
plot_fixed_time(ax_10min, 600, 'Temperature at 10 min')

# 动画更新函数
def update(frame):
    cax.set_array(T_all[frame])
    ax_anim.set_title(f'Time = {frame*dt/60:.2f} min')
    return cax,

ani = animation.FuncAnimation(fig, update, frames=range(0, total_steps, int(1/dt)), interval=50, blit=False)

plt.show()

其中，可在常量定义部分修改时间步长，如设置为0.1，即可获得非平衡模拟过程

模拟结果