热传导方程以及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求解
对各点差分化， 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)之间可以相互表示
我们可以用图直观表示
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://javaforall.net/218237.html原文链接：https://javaforall.net
