pde分享------一维热传导方程的理论与应用1

前言

       随着科技的进步和工程应用的不断拓展，对热传导方程的理解和应用已经成为现代工程技术中的一个重要方面。在我们日常生活中，热传导现象无处不在。无论是烹饪时的火焰，还是冬季保暖的衣物，热传导都在其中发挥着关键作用。然而，要深入理解并应用热传导，我们需要依赖于数学模型来描述这一物理现象。其中，最基本也最重要的模型之一就是一维热传导方程。
        一维热传导方程是一个描述热量如何在物体中传递的数学方程。通过这个方程，我们可以了解物体内部温度的分布情况，以及温度随时间和空间的变化规律。这一模型的建立，为我们提供了分析和解决实际热传导问题的有力工具。
在接下来的篇幅中，我们将深入探讨一维热传导方程的理论和应用。通过详尽的数学推导来得到这些方程，此外，还将介绍有限差分法（显式格式）来得到一维热传导方程的数值解。

       首先，引入这样一个场景，如下图所示，最左边有一个恒温热源100℃，其通过一根杆将热量传递到小人的皮肤上，在此之前，需要做一些假设：

1. 只考虑棒的长度方向上的热传导，忽略其他方向。这样，就将情景简化成一个维度上的问题。
2. 杆的末端与人体皮肤表面的热交换情况暂时不做考虑，即热源的热量传导杆的末端即结束，杆末端与皮肤的传热方程将在后面的文章中介绍。
   
          人皮肤表面初始温度37℃的初始条件已经确定，随后需要确定其他区域的初始条件，不妨设这根杆上在t=0时刻的初始温度也为37℃。至此，所有的边界条件和初始条件都已经确立，可以进行下面的模型建立。

一维热传导方程模型的建立

       设温度为时间空间的函数：$u=u\left(x,t\right)$，这个式子代表x位置处，t时刻所具有的温度。
热传导方程基于能量守恒定律推导。在微小的时间内，棒的某一部分的温度变化与该部分的热流量梯度成正比。其热量变化满足如下方程：
$$\frac{\partial u\left(x,t\right)}{\partial t}=\alpha \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}$$
其中，$\alpha =\frac{k}{\rho c},$ $k$ 是物体的导热系数，$\rho$ 是物体的密度，c是物体的比热。
$\alpha$ 被称为热扩散系数，常见的热扩散系数如下表所示：

       这是一个偏微分方程（PDE）,为了求出具体物理场景下的数值，还需要确定一些初始条件和边界条件，以及其他参数。
       设杆子长度是1m,并且认为杆子为银材质，那么从上表就可以查出，此银材质杆子的热扩散系数$\alpha =149\times 10^{-6}$ ，并且设置传热仿真的时间为500s
       初始条件：物体内部各点在初始时刻t = 0 的温度分布$u\left(x,0\right)=37,\left(x\ne 0\right)$
       边界条件：边界x=0处有一个恒温热源$u\left(0,t\right)=100$
       中间层满足热传导方程：$\frac{\partial u\left(x,t\right)}{\partial t}=\alpha \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}$ ，
其中，$\alpha =\frac{k}{\rho c}$ ，$k$是物体的导热系数，$\rho$ 是物体的密度，$c$是物体的比热。
综上所述，整个物理过程满足如下方程组：
$$\{\begin{array}{l}\frac{\partial u\left(x,t\right)}{\partial t}=\alpha \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\\ u\left(0,t\right)=100\\ u\left(x,0\right)=37\\ \alpha =\frac{k}{\rho c}=149\times 10^{-6}\end{array}$$
       至此，这个物理过程中，所有的变化形式都可以用上述方程组来描述，只要求出这组偏微分方程的解，即可清晰地知晓整个热传导的过程，在求解的过程中，采用有限差分法，有限差分法包含了显式格式和隐式格式，要注意的是，显式格式是有收敛条件的，如果不满足收敛条件，求出的解是错误的。而有限差分法的隐式格式则是无条件收敛的，但是其表示比较复杂。有限差分法的显式格式的好处是表示简单，清晰易懂，因此，下文求解时采用有限差分法的显式格式。

一维热传导方程模型的求解

采用有限差分法，显式格式：
有限差分法的基本原理即使用近似方法处理微分方程中的微分项。
首先，对求解区域进行网格划分，我们将杆子的长度区间[0,1]划分为等距的$N_x$ 个部分，每一个小网格的长度为$\Delta x=\frac{1}{N_x}$ ，将仿真时间500s划分为等间隔的$N_t$ 个部分，每一个小时间步长的长度为$\Delta t=\frac{500}{N_t}$
接下来，对原来的热传导方程进行差分，用差分来代替微分运算。
热传导方程差分：对二阶差分项，我们使用二阶中心差商，上述方程可以近似写成
$$\frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\Delta t}=\alpha \frac{u\left(x_{i+1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i-1},t_j\right)}{\Delta x^2}$$
化简可得：
$$u\left(x_i,t_{j+1}\right)=\frac{\alpha \cdot \Delta t}{\Delta x^2}\cdot \left(u\left(x_{i+1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i-1},t_j\right)\right)+u\left(x_i,t_j\right)$$
边界条件差分：
$$\{\begin{array}{l}u\left(x_i,0\right)=37,i=1,2,...,N_x\\ u\left(0,t_j\right)=100,u\left(N_x,t_j\right)=37,j=1,2,...,N_t\end{array}$$
综上所述，整个物理过程差分后形式如下所示：
$$\{\begin{array}{l}u\left(x_i,t_{j+1}\right)=\frac{\alpha \cdot \Delta t}{\Delta x^2}\cdot \left(u\left(x_{i+1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i-1},t_j\right)\right)+u\left(x_i,t_j\right)\\ u\left(x_i,0\right)=37,i=1,2,...,N_x\\ u\left(0,t_j\right)=100,u\left(N_x,t_j\right)=37,j=1,2,...,N_t\end{array}$$
       也就是说，用前一个时刻$t_j$ 中，位置$x_{i-1},x_i,x_{i+1}$ 处的温度来得出下一个时刻$t_{j+1}$ 中，位置$x_i$ 处的位置。过程较为简单，用matlab写个循环即可实现，设置好相应的边界条件，初始条件及参数，即可对一维传热情况进行数值模拟。
       但是，上文也提到过，有限差分法的显式格式是有收敛条件的，根据蔡志杰老师在杂志《数学建模及其应用》中一篇名为“高温作业专用服装设计”的文章指出，该情况下需要满足$\frac{\alpha \cdot \Delta t}{\Delta x^2}\le \frac{1}{2}$ ，方程的解才会收敛，本代码中已经调整过相应的参数，解出来的热传导方程是收敛的。
Matlab代码如下：

clc
clear all
close all
Nx=100;
Nt=100;
alpha=149e-6;
total_location=1;
%杆子长度1m
total_time=500;
%仿真500秒
dx=1/Nx;
dt=1/Nt;
delta=dt/(dx*dx);
a=alpha*(dt/(dx*dx))
%a=0.01;
u=zeros(total_location/dx,total_time/dt);%初始化温度矩阵
u(:,1)=37;
u(1,:)=100;
u(end,:)=37;
%有限差分法显式格式计算各层的温度
for j=1:1:total_time/dt-1
     for i=2:1:(total_location/dx)-1
         u(i,j+1)=(1-2*a)*u(i,j)+a*(u(i+1,j)+u(i-1,j));
     end
 end
figure 
mesh(u);%求出不同时间t和不同位置x下杆的温度，并可视化
xlabel('时间/s');
ylabel('杆子长度/m');
zlabel('温度/℃');
grid on

figure 
plot(u(end-1,:));%画出杆末端温度的变化情况
xlabel('时间/s');
ylabel('温度/m');
title('最后一层温度');
grid on

杆子上空间-时间温度分布三维图如下所示：

杆子末端温度随时间变化图如下所示，可以看出，随着时间推移，杆子末端温度逐步上升中……