【数学建模】导热微分方程-Crank-Nicolson 方法

本篇文章提供由 Crank-Nicolson 方法求解一维非稳态导热微分方程的代码示例，对原理不会展开讲解。

模型建立

reyuanchuandit0818_

假设有个各向同性，导热系数固定，在热源一侧均匀导热的墙壁。即可建立一维非稳态无内热源导热微分方程： $$ \frac{\partial t}{\partial \tau} = \alpha \frac{\partial^2 t}{\partial x^2} $$

$\alpha$是热扩散系数，单位是$m^2/s$
$t$是温度，单位是℃（或者K）
$\tau$是时间，单位是$s$

初始条件： $$ t_{\tau=0} = 25 $$ 即整个墙壁的温度分布均为25℃。

边界条件： $$ t_{x=0} = \ln{(1+100\tau)} + 25 $$ 即在热源一侧的壁面温度随时间变化已知。

为了数值求解上述偏微分方程，使用Crank-Nicolson 方法。这是一种常用于数值求解偏微分方程的方法，它结合了向前欧拉方法和向后欧拉方法的特点，具有较好的数值稳定性和精度。

Crank-Nicolson 方法的离散形式如下： $$ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{\alpha}{2 (\Delta t)^2} \left[(u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}) + (u_{i+1}^n - 2u_i^n + u_{i-1}^n)\right] $$ 然而，实现 Crank-Nicolson 方法需要解一个线性方程组，通常涉及矩阵求逆等操作，这里不做讲解。

MATLAB代码示例

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49


clear;clc
% 参数和网格设置
L = 1;          % 空间范围长度
T = 1;          % 时间总长
Nx = 50;        % 空间网格点数
Nt = 100;       % 时间步数
alpha = 0.3;    % 热扩散系数

dx = L / Nx;
dt = T / Nt;

x = linspace(0, L, Nx+1);
t = linspace(0, T, Nt+1);

% 初始化温度场
u = zeros(Nx+1, Nt+1);
u(:,:) = 25;

% 设置定解条件
u(1,:) = log(1+100.*linspace(0,T,length(u)))+25;

% Crank-Nicolson 方法求解
A = zeros(Nx-1, Nx-1);
for i = 1:Nx-1
    A(i, i) = 1 + alpha*dt/(2*dx^2);
    if i > 1
        A(i, i-1) = -alpha*dt/(4*dx^2);
    end
    if i < Nx-1
        A(i, i+1) = -alpha*dt/(4*dx^2);
    end
end

for n = 1:Nt
    b = u(2:Nx, n);
    b(1) = b(1) + alpha*dt/(4*dx^2) * u(1, n+1);
    b(end) = b(end) + alpha*dt/(4*dx^2) * u(Nx+1, n+1);
    u(2:Nx, n+1) = A \ b;
end

% 可视化
[X, T] = meshgrid(x, t);
surf(X, T, u');
xlabel('x');
ylabel('t');
zlabel('Temperature');
title('Temperature Distribution');
colormap(flipud(autumn))
view(45,45)