热传导方程以及Matlab求解
简略推导
设 u = u ( x ⃗ , t ) u=u(\vec{x},t) u=u(x,t)表示某一个均匀物体 Ω ⊂ R 3 \Omega\subset\mathbb{R}^3 Ω⊂R3在位置 x ⃗ \vec{x} x，时刻 t t t的温度（单位： K K K）, c c c表示比热（单位： J / ( k g ⋅ K ) J/(kg\cdot K) J/(kg⋅K)）, ρ \rho ρ表示密度（单位： k g / m 3 kg/m^3 kg/m3），令 V V V为 Ω \Omega Ω内的任何光滑区域， q ⃗ \vec{q} q表示温度场的热流密度（单位： W / m 2 W/m^2 W/m2）且 f 0 = f 0 ( x ⃗ , t ) f_0=f_0(\vec{x},t) f0=f0(x,t)表示在位置 x x x，时刻 t t t产生的热量密度，由能量守恒定律我们得到
d
d
t
∭
V
c
ρ
u
(
x
⃗
,
t
)
d
v
+
∬
∂
V
q
⃗
⋅
n
⃗
d
S
=
∭
V
ρ
f
0
d
v
\frac{d}{dt}\iiint\limits_{V}c\rho u(\vec{x},t)dv+\iint\limits_{\partial V}\vec{q}\cdot\vec{n}dS=\iiint\limits_{V}\rho f_0dv
dtdV∭cρu(x,t)dv+∂V∬q⋅ndS=V∭ρf0dv
其中
n
⃗
\vec{n}
n为
∂
V
\partial V
∂V的单位外法向量，利用Gauss公式我们知道
c
ρ
∭
V
∂
u
(
x
⃗
,
t
)
∂
t
d
v
+
∭
V
d
i
v
q
⃗
d
v
=
∭
V
ρ
f
0
(
x
⃗
,
t
)
d
v
c\rho\iiint\limits_{V}\frac{\partial u(\vec{x},t)}{\partial t}dv+\iiint\limits_{V}\mathrm{div}\vec{q}dv=\iiint\limits_{V}\rho f_0(\vec{x},t)dv
cρV∭∂t∂u(x,t)dv+V∭divqdv=V∭ρf0(x,t)dv
由于区域
V
V
V的任意性，我们有
c
ρ
∂
u
(
x
⃗
,
t
)
∂
t
+
d
i
v
q
⃗
=
ρ
f
0
(
x
⃗
,
t
)
c\rho \frac{\partial u(\vec{x},t)}{\partial t}+\mathrm{div}\vec{q}=\rho f_0(\vec{x},t)
cρ∂t∂u(x,t)+divq=ρf0(x,t)
由Fourier热传导定律，温度场中的热流密度与温度梯度成正比，而且热量总是从高温处流向低温处
q
⃗
=
−
k
(
∂
u
(
x
⃗
,
t
)
∂
x
1
,
∂
u
(
x
⃗
,
t
)
∂
x
2
,
…
,
∂
u
(
x
⃗
,
t
)
∂
x
n
)
\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})
q=−k(∂x1∂u(x,t),∂x2∂u(x,t),…,∂xn∂u(x,t))
k
>
0
k>0
k>0是物体的导热系数，代入上式之后，我们得到热方程
∂
u
(
x
⃗
,
t
)
∂
t
−
a
2
(
∂
2
u
(
x
⃗
,
t
)
∂
x
1
2
+
∂
2
u
(
x
⃗
,
t
)
∂
x
2
2
+
⋯
+
∂
2
u
(
x
⃗
,
t
)
∂
x
n
2
)
=
f
(
x
⃗
,
t
)
\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}+\dots+\frac{\partial^2 u(\vec{x},t)}{\partial x_n^2})= f(\vec{x},t)
∂t∂u(x,t)−a2(∂x12∂2u(x,t)+∂x22∂2u(x,t)+⋯+∂xn2∂2u(x,t))=f(x,t)
这里
a
=
k
c
ρ
,
f
(
x
⃗
,
t
)
=
f
0
(
x
⃗
,
t
)
c
a=\sqrt{\frac{k}{c\rho}},f(\vec{x},t)=\frac{f_0(\vec{x},t)}{c}
a=cρk,f(x,t)=cf0(x,t)
初始条件和边界条件的物理意义
初始条件：给出物体内部各点在初始时刻
t
=
0
t=0
t=0的温度分布
u
(
x
⃗
,
0
)
=
φ
(
x
⃗
)
,
x
∈
Ω
u(\vec{x},0)=\varphi(\vec{x}),x\in\Omega
u(x,0)=φ(x),x∈Ω
其中
φ
(
x
)
\varphi(x)
φ(x)是已知函数
边值条件：给出物体边界 ∂ Ω \partial \Omega ∂Ω在时间 t > 0 t>0 t>0的温度分布或受周围介质的影响情况，通常分为以下三类
(1)已知边界
∂
Ω
\partial \Omega
∂Ω的温度分布
u
(
x
⃗
,
t
)
=
g
(
x
⃗
,
t
)
,
x
∈
∂
Ω
,
t
≥
0
u(\vec{x},t)=g(\vec{x},t),x\in\partial\Omega,t\ge0
u(x,t)=g(x,t),x∈∂Ω,t≥0
当
g
g
g为常数时，表示物体的边界保持恒温。
(2)已知通过物体边界
∂
Ω
\partial \Omega
∂Ω流入或流出物体
Ω
\Omega
Ω的热量
k
∂
∂
n
⃗
u
(
x
⃗
,
t
)
=
g
(
x
⃗
,
t
)
,
x
∈
∂
Ω
,
t
≥
0
k\frac{\partial}{\partial\vec{n}}u(\vec{x},t)=g(\vec{x},t),x\in\partial\Omega,t\ge0
k∂n∂u(x,t)=g(x,t),x∈∂Ω,t≥0
其中
n
⃗
\vec{n}
n是
∂
Ω
\partial\Omega
∂Ω的单位外法向量，当
g
≥
0
g\ge0
g≥0时表示热量流入，
g
≤
0
g\le0
g≤0时表示热量流出，
g
=
0
g=0
g=0时表示物体绝热。
(3)已知通过物体边界
∂
Ω
\partial\Omega
∂Ω与周围介质的热交换强度
k
∂
∂
n
⃗
u
(
x
⃗
,
t
)
=
α
0
(
x
⃗
,
t
)
(
g
0
(
x
⃗
,
t
)
−
u
(
x
⃗
,
t
)
)
,
x
∈
∂
Ω
,
t
≥
0
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\ge0
k∂n∂u(x,t)=α0(x,t)(g0(x,t)−u(x,t)),x∈∂Ω,t≥0
或
∂
∂
n
⃗
u
(
x
⃗
,
t
)
+
α
(
x
⃗
,
t
)
u
(
x
⃗
,
t
)
=
g
(
x
⃗
,
t
)
,
x
∈
∂
Ω
,
t
≥
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\ge0
∂n∂u(x,t)+α(x,t)u(x,t)=g(x,t),x∈∂Ω,t≥0
其中
n
⃗
\vec{n}
n是
∂
Ω
\partial\Omega
∂Ω的单位外法向量，
g
0
g_0
g0表示周围介质的温度，
α
0
≥
0
\alpha_0\ge0
α0≥0表示热交换系数，
α
=
α
0
k
>
0
,
g
=
g
0
k
\alpha=\frac{\alpha_0}{k}>0,g=\frac{g_0}{k}
α=kα0>0,g=kg0
有限差分法以及Matlab求解
我们给出一个例题
{
∂
u
(
x
,
t
)
∂
t
=
α
∂
2
u
(
x
,
t
)
∂
x
2
u
(
0
,
t
)
=
T
1
=
65
,
u
(
x
d
,
t
)
=
T
2
=
37
u
(
x
,
0
)
=
20
(
x
≠
0
)
α
=
k
ρ
c
,
k
=
0.082
,
ρ
=
300
,
c
=
1377
\begin{cases} \frac{\partial u(x,t)}{\partial t}=\alpha\frac{\partial^2 u(x,t)}{\partial x^2} \\u(0,t)=T_1=65,\quad u(x_d,t)=T_2=37 \\u(x,0)=20(x\ne0) \\\alpha=\frac{k}{\rho c},k=0.082,\rho=300,c=1377 \end{cases}
⎩
⎨
⎧∂t∂u(x,t)=α∂x2∂2u(x,t)u(0,t)=T1=65,u(xd,t)=T2=37u(x,0)=20(x=0)α=ρck,k=0.082,ρ=300,c=1377
该题由2018年数学建模国赛A题改编而成，
u
u
u是在
x
=
0
x=0
x=0处有一个65度的热源，经过时间
t
t
t后，距离热源
d
d
d处的温度。
题目中的边界条件是Dirchilet边界条件，
u
(
0
,
t
)
=
T
1
=
65
u(0,t)=T_1=65
u(0,t)=T1=65的物理意义是任何时刻
t
t
t距离恒定热源
0
0
0的温度是
T
1
=
65
T_1=65
T1=65，
u
(
x
d
,
t
)
=
T
2
=
37
u(x_d,t)=T_2=37
u(xd,t)=T2=37的物理意义是防护服的厚度是
x
d
x_d
xd，即距离热源
x
d
x_d
xd的温度恒为37度，即人体体温。
我们用一个图表示题目中给出的边界条件。
u
(
x
,
0
)
=
20
u(x,0)=20
u(x,0)=20的物理意义是在
t
=
0
t=0
t=0时刻防护服内部的任何位置都是环境温度20度。
上图中，黄色的点代表已知的点，黑色的点代表未知的点。
对各点差分化，
x
x
x轴设置
N
x
N_x
Nx个分点，
t
t
t轴每个单位长度设置
N
t
N_t
Nt个分点。
{
x
i
=
i
N
x
,
i
=
1
,
2
,
…
,
N
x
t
j
=
j
N
t
,
j
=
1
,
2
,
…
,
∞
x
i
+
1
−
x
i
=
d
x
=
x
d
N
x
t
j
+
1
−
t
j
=
d
t
=
1
N
t
\begin{cases} x_i=\frac{i}{N_x},\quad i=1,2,\dots,N_x \\t_j=\frac{j}{N_t},\quad 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{cases}
⎩
⎨
⎧xi=Nxi,i=1,2,…,Nxtj=Ntj,j=1,2,…,∞xi+1−xi=dx=Nxxdtj+1−tj=dt=Nt1
对二阶差分项，我们使用二阶中心差商，上述方程可以近似写成
u
(
x
i
,
t
j
+
1
)
−
u
(
x
i
,
t
j
)
d
t
=
α
u
(
x
i
+
1
,
t
j
)
−
2
u
(
x
i
,
t
j
)
+
u
(
x
i
−
1
,
t
j
)
d
x
2
\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}
dtu(xi,tj+1)−u(xi,tj)=αdx2u(xi+1,tj)−2u(xi,tj)+u(xi−1,tj)
这说明，
(
i
+
1
,
j
)
,
(
i
,
j
)
,
(
i
,
j
+
1
)
,
(
i
,
j
−
1
)
(i+1,j),(i,j),(i,j+1),(i,j-1)
(i+1,j),(i,j),(i,j+1),(i,j−1)之间可以相互表示
我们可以用图直观表示
方程可以改写成
{
u
(
x
i
,
t
j
+
1
)
=
(
1
−
2
α
d
t
d
x
2
)
u
(
x
i
,
t
j
)
+
α
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
,
…
,
N
x
u
(
0
,
t
j
)
=
65
,
u
(
N
x
,
t
j
)
=
37
,
j
=
0
,
1
,
2
,
…
,
∞
\begin{cases} 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,\quad u(N_x,t_j)=37,& j=0,1,2,\dots,\infty \end{cases}
⎩
⎨
⎧u(xi,tj+1)=(1−2dx2αdt)u(xi,tj)+dx2αdt(u(xi+1,tj)+u(xi−1,tj))u(xi,0)=20,u(0,tj)=65,u(Nx,tj)=37,i=1,2,…,Nxj=0,1,2,…,∞
我们根据下面一行已知点，求解出上面一行未知点，逐层向上求解，可以使用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
文章来源：https://www.toymoban.com/news/detail-461226.html
文章来源地址https://www.toymoban.com/news/detail-461226.html
到了这里，关于热传导方程以及Matlab求解的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
