热传导方程以及Matlab求解

热传导方程以及Matlab求解

简略推导

设$u=u(\vec{x},t)$表示某一个均匀物体$\Omega \subset {\mathbb{R}}^3$在位置$\vec{x}$，时刻$t$的温度（单位：$K$）,$c$表示比热（单位：$J/(kg\cdot K)$）,$\rho$表示密度（单位：$kg/m^3$），令$V$为$\Omega$内的任何光滑区域，$\vec{q}$表示温度场的热流密度（单位：$W/m^2$）且$f_0=f_0(\vec{x},t)$表示在位置$x$，时刻$t$产生的热量密度，由能量守恒定律我们得到

$$\frac{d}{dt}\underset{V}{\iiint }c\rho u(\vec{x},t)dv+\underset{\partial V}{\iint }\vec{q}\cdot \vec{n}dS=\underset{V}{\iiint }\rho f_{0}dv$$
其中$\vec{n}$为$\partial V$的单位外法向量，利用Gauss公式我们知道
$$c\rho \underset{V}{\iiint }\frac{\partial u(\vec{x},t)}{\partial t}dv+\underset{V}{\iiint }\mathrm{div}\vec{q}dv=\underset{V}{\iiint }\rho f_0(\vec{x},t)dv$$
由于区域$V$的任意性，我们有
$$c\rho \frac{\partial u(\vec{x},t)}{\partial t}+\mathrm{div}\vec{q}=\rho f_0(\vec{x},t)$$
由Fourier热传导定律，温度场中的热流密度与温度梯度成正比，而且热量总是从高温处流向低温处
$$\vec{q}=-k(\frac{\partial u(\vec{x},t)}{\partial x_1},\frac{\partial u(\vec{x},t)}{\partial x_2},\dots ,\frac{\partial u(\vec{x},t)}{\partial x_n})$$
$k>0$是物体的导热系数，代入上式之后，我们得到热方程
$$\frac{\partial u(\vec{x},t)}{\partial t}-a^2(\frac{{\partial }^{2}u(\vec{x},t)}{\partial x_1^2}+\frac{{\partial }^{2}u(\vec{x},t)}{\partial x_2^2}+\cdots +\frac{{\partial }^{2}u(\vec{x},t)}{\partial x_n^2})=f(\vec{x},t)$$
这里$a=\sqrt{\frac{k}{c\rho }},f(\vec{x},t)=\frac{f_0(\vec{x},t)}{c}$

初始条件和边界条件的物理意义

初始条件：给出物体内部各点在初始时刻$t=0$的温度分布
$$u(\vec{x},0)=\phi (\vec{x}),x\in \Omega$$
其中$\phi (x)$是已知函数

边值条件：给出物体边界$\partial \Omega$在时间$t>0$的温度分布或受周围介质的影响情况，通常分为以下三类

(1)已知边界$\partial \Omega$的温度分布
$$u(\vec{x},t)=g(\vec{x},t),x\in \partial \Omega ,t\ge 0$$
当$g$为常数时，表示物体的边界保持恒温。

(2)已知通过物体边界$\partial \Omega$流入或流出物体$\Omega$的热量
$$k\frac{\partial }{\partial \vec{n}}u(\vec{x},t)=g(\vec{x},t),x\in \partial \Omega ,t\ge 0$$
其中$\vec{n}$是$\partial \Omega$的单位外法向量，当$g\ge 0$时表示热量流入，$g\le 0$时表示热量流出，$g=0$时表示物体绝热。

(3)已知通过物体边界$\partial \Omega$与周围介质的热交换强度
$$k\frac{\partial }{\partial \vec{n}}u(\vec{x},t)={\alpha }_0(\vec{x},t)(g_0(\vec{x},t)-u(\vec{x},t){)}_{,}x\in \partial \Omega ,t\ge 0$$
或
$$\frac{\partial }{\partial \vec{n}}u(\vec{x},t)+\alpha (\vec{x},t)u(\vec{x},t)=g(\vec{x},t),x\in \partial \Omega ,t\ge 0$$
其中$\vec{n}$是$\partial \Omega$的单位外法向量，$g_0$表示周围介质的温度，${\alpha }_0\ge 0$表示热交换系数，$\alpha =\frac{{\alpha }_0}{k}>0,g=\frac{g_0}{k}$

有限差分法以及Matlab求解

我们给出一个例题
$$\{\begin{array}{l}\frac{\partial u(x,t)}{\partial t}=\alpha \frac{{\partial }^{2}u(x,t)}{\partial x^2}\\ u(0,t)=T_1=65,u(x_d,t)=T_2=37\\ u(x,0)=20(x\ne 0)\\ \alpha =\frac{k}{\rho c},k=0.082,\rho =300,c=1377\end{array}$$
该题由2018年数学建模国赛A题改编而成，$u$是在$x=0$处有一个65度的热源，经过时间$t$后，距离热源$d$处的温度。
题目中的边界条件是Dirchilet边界条件，$u(0,t)=T_1=65$的物理意义是任何时刻$t$距离恒定热源$0$的温度是$T_1=65$，$u(x_d,t)=T_2=37$的物理意义是防护服的厚度是$x_d$，即距离热源$x_d$的温度恒为37度，即人体体温。
我们用一个图表示题目中给出的边界条件。$u(x,0)=20$的物理意义是在$t=0$时刻防护服内部的任何位置都是环境温度20度。

上图中，黄色的点代表已知的点，黑色的点代表未知的点。

对各点差分化，$x$轴设置$N_x$个分点，$t$轴每个单位长度设置$N_t$个分点。
$$\{\begin{array}{l}{x}_i=\frac{i}{N_x},i=1,2,\dots ,N_x\\ t_j=\frac{j}{N_t},j=1,2,\dots ,\infty \\ x_{i+1}-x_i=d_x=\frac{x_d}{N_x}\\ t_{j+1}-t_j=d_t=\frac{1}{N_t}\end{array}$$
对二阶差分项，我们使用二阶中心差商，上述方程可以近似写成
$$\frac{u(x_i,t_{j+1})-u(x_i,t_j)}{d_t}=\alpha \frac{u(x_{i+1},t_j)-2u(x_i,t_j)+u(x_{i-1},t_j)}{d_x^2}$$
这说明，$(i+1,j),(i,j),(i,j+1),(i,j-1)$之间可以相互表示

我们可以用图直观表示

方程可以改写成
$$\{\begin{array}{ll}u(x_i,t_{j+1})=(1-2\frac{\alpha d_t}{d_x^2})u(x_i,t_j)+\frac{\alpha d_t}{d_x^2}(u(x_{i+1},t_j)+u(x_{i-1},t_j))& \\ u(x_i,0)=20,& i=1,2,\dots ,N_x\\ u(0,t_j)=65,u(N_x,t_j)=37,& j=0,1,2,\dots ,\infty \end{array}$$
我们根据下面一行已知点，求解出上面一行未知点，逐层向上求解，可以使用Matlab求解。

Nx=100;Nt=100;%我们对x轴分成Nx段，有Nx+1个点，对t轴分成Nt段，有Nt+1个点
xd=1;t=1000;
dt=1/Nt;dx=xd/Nx;
alpha=1000*0.082/(300*1377);%热传导系数，可以自行调整，alpha越大表示传热越快
u=zeros(Nx+1,floor(t*Nt)+1);%生成(Nx+1)*(Nt+1)的矩阵
u(end,:)=37;u(:,1)=20;
u(1,:)=65;
for j=1:floor(t*Nt)
    for i=2:Nx
        u(i,j+1)=(1-2*alpha*dt/dx^2)*u(i,j)+(alpha*dt/dx^2)*(u(i+1,j)+u(i-1,j));
    end
end
f=@(x,t)u(floor(x*Nx/xd)+1,floor(t*Nt)+1);
figure
X=0:0.05:1;[~,xl]=size(X);
T=0:1:1000;[~,Tl]=size(T);
[xx,tt]=meshgrid(X,T);
Z=zeros(xl,Tl);
for ii=1:xl
    for jj=1:Tl
        Z(ii,jj)=f(X(ii),T(jj));
    end
end
mesh(xx,tt,Z')
xlabel('与热源之间的距离d')
ylabel('经过的时间t')
zlabel('温度')
%下面我们要找出x=d时，温度随时间的变化趋势
figure
for d=0.01:0.01:0.12
    subplot(3,4,100*d)
    t_x_d=0:1:1000;
    [~,len]=size(t_x_d);
    U_d=zeros(1,len);
    for i=1:len
        U_d(i)=f(d,t_x_d(i));
    end
    plot(t_x_d,U_d);
    xlabel('经过的时间t')
    ylabel('温度')
    title(['d=' num2str(d)])
end
figure
for d=0.2:0.05:0.95
    subplot(4,4,20*d-3)
    t_x_d=0:1:1000;
    [~,len]=size(t_x_d);
    U_d=zeros(1,len);
    for i=1:len
        U_d(i)=f(d,t_x_d(i));
    end
    plot(t_x_d,U_d);
    xlabel('经过的时间t')  
    ylabel('温度')
    title(['d=' num2str(d)])
end