热传导方程差分格式
热传导方程的求解及其应用
热传导方程的求解及其应用热传导是指物质内部由高温区向低温区传递热量的过程，是自然界中十分普遍的现象。
为了更好地理解和研究这一过程，我们需要借助数学模型来描述和求解热传导过程，其中最常用的数学模型就是热传导方程。
一、热传导方程的数学模型热传导方程是描述物质内部温度变化随时间和空间的变化而变化的偏微分方程。
它可以描述均质物质内部的热量传递，以及介质中的温度变化。
热传导方程的数学表示式如下：$$ \frac{\partial u}{\partial t}=\alpha \nabla^2 u $$其中，$u$表示物质内部温度的分布，$t$表示时间，$\alpha$表示热扩散系数，$\nabla^2$表示拉普拉斯算子，表示温度分布的曲率。
二、热传导方程的求解方法热传导方程是一个偏微分方程，需要借助一定的数学方法才能求解。
下面简要介绍两种常见的求解方法：1.分离变量法分离变量法是求解偏微分方程的常见方法之一。
对于热传导方程，我们通常采用分离变量法将其转化为两个方程：$$ \frac{1}{\alpha}\frac{\partial u}{\partial t}= \nabla^2 u $$设$u(x,t)=f(x)g(t)$，代入上式得：$$ \frac{1}{\alpha}\frac{g'(t)}{g(t)}= \frac{f''(x)}{f(x)}=\lambda $$其中，$\lambda$为待定常数，$f(x)$和$g(t)$分别为$x$和$t$的函数。
将上述两个方程分别求解，可以得到形如下面的解：$$ u(x,t)=\sum_{n=1}^{\infty}c_nexp(-\lambda_n\alphat)sin(\frac{n\pi x}{L}) $$其中，$\lambda_n$为常数，$L$为问题的区间长度。
2.有限差分法有限差分法是一种常见的数值求解方法，可以用来求解各种偏微分方程，包括热传导方程。
6第六讲 典型模型方程-热传导方的差分格式程
ADE Methods
同一时间步，同时左右扫描
6、Hopscotch Method
Comments
2、Richardaon’s Method: CTCS
3、Simple Implicit(Laasonen) Method
算子表示：
其中：
未知，每一个时间步长需要求解三对角方程组
放大因子：
4、Crank-Nicolson Method: famous
1 N 1 1 uN uN j 1 2u j j 1
where
G 1
3-D: ADI Methods
假设
∆t左移，右端合并即为C-N格式
其中 rx x 2 , ry y 2
at
at
3-D: ADI Methods
4、Splitting or Fractional-step Methods
5、ADE Methods
1-D:
或
Simple Explicit Method
差分格式的放大因子与精确解的放大因子比较
精度高
Simple Explicit Method
Simple Explicit Method: Example
Simple Explicit Method: Example
没有相位误差，幅值误差：1.88%
2-D: Crank-Nicolson Scheme
通常采用迭代方 法，需要比三对 角更多的计算资 源
3、2-D: ADI Methods
time
1/ 2 n 1 / 2 1/ 2 uin , uin 1, j , ui , j 1, j
1 n 1 n 1 uin, j 1 , ui , j , ui , j 1
热传导方程的一种高精度O(τ 2+h 4)阶差分格式
“ 一 （ ＋ 一 ２ ｉ ｕ一 ） ， ｕ ＋ ｊ。 ，
验估 计和 数值 例子 。
１ 网格 剖 分
取空间步长 ｈ ／ 和时间步长 Ｚ＝ ／ 其 —Ｚ Ｍ －＝ Ｎ， ＝Ｔ
有限差分法及热传导数值计算
有限差分法及热传导数值计算有限差分法（finite difference method）是一种常用的数值计算方法，可以用于求解热传导问题。
它基于热传导方程，通过将连续的热传导问题离散化成离散网格上的代数方程组，然后利用数值迭代方法求解方程组，得到热传导问题的数值解。
热传导方程描述了热量在物体内部传导的过程，它可以写成以下形式：∂T/∂t=α∇²T其中，T是温度场的分布，α是热扩散系数，∇²是拉普拉斯算子。
为了使用有限差分法求解热传导问题，我们需要将时间和空间进行离散化。
时间上，我们将连续的时间区间[0,T]分成N个子区间，每个子区间的长度为Δt，表示为t_i=iΔt，其中i=0,1,2,...,N。
空间上，我们将研究区域Ω划分为M个离散节点，每个节点的坐标为x_j，表示为x_j=jΔx，其中j=0,1,2,...,M。
在离散化后，我们可以用差分近似的方式来近似热传导方程。
对于时间上的导数，我们可以使用前向差分，即∂T(x_j,t_i)/∂t≈(T(x_j,t_{i+1})-T(x_j,t_i))/Δt对于空间上的二阶导数，我们可以使用中心差分，即∇²T(x_j,t_i)≈(T(x_{j-1},t_i)-2T(x_j,t_i)+T(x_{j+1},t_i))/Δx²将上述差分近似带入热传导方程中，我们可以得到如下的差分方程：(T(x_j,t_{i+1})-T(x_j,t_i))/Δt=α*(T(x_{j-1},t_i)-2T(x_j,t_i)+T(x_{j+1},t_i))/Δx²重新整理得到：T(x_j,t_{i+1})=T(x_j,t_i)+α*Δt*(T(x_{j-1},t_i)-2T(x_j,t_i)+T(x_{j+1},t_i))/Δx²这个差分方程可以用于迭代求解热传导问题。
我们可以根据初始条件和边界条件，从t=0的初始时刻开始，按照时间步长Δt进行迭代计算。
热传导方程的差分格式汇总
热传导方程的差分格式汇总1.显式差分格式：显式差分格式是最简单的一种方法，通过将导热方程时间和空间上的导数进行近似，引入差分算子，将方程转化为差分格式。
其中最常见的差分格式有：a. 前向差分法（Forward Difference Method）：利用当前节点和其相邻节点的温度值进行计算。
例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+α(u(i+1,j)-2u(i,j)+u(i-1,j))b. 后向差分法（Backward Difference Method）：利用当前节点和其相邻节点的温度值进行计算。
例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+α(u(i+1,j+1)-2u(i,j+1)+u(i-1,j+1))c. 中心差分法（Central Difference Method）：利用当前和其相邻节点的温度值进行计算。
例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+α(u(i+1,j)-2u(i,j)+u(i-1,j))+β(u(i+1,j)-u(i-1,j))其中α和β是时间和空间步长的比例因子。
2.隐式差分格式：显式差分格式具有较大的稳定性限制。
为了克服这个问题，可以使用隐式差分格式，其中使用下一个时间步长的温度值来求解当前时间步长。
常见的隐式差分格式有：a. C-N差分法（Crank-Nicolson Method）：利用前后两个时间步长的温度值进行计算。
例如，在一维离散情况下，可以使用公式：u(i,j+1)=u(i,j)+0.5α(u(i+1,j+1)-2u(i,j+1)+u(i-1,j+1))+0.5α(u(i+1,j)-2u(i,j)+u(i-1,j))b. 力学模拟法（Finite Element Method）：将空间离散化后，通过引入有限元方法，将热传导问题转化为线性方程组，再通过求解线性方程组得到温度分布。
热传导与有限差分
(6)
研究生课程《高等工程热力学与传热学》电子教案
三类边界条件
第 1 类边界条件： 给定边界上的温度。 第 2 类边界条件： 给定边界上的法向热流密度。 第 3 类边界条件： 给定外部介质的温度和给定边界上的对流换热面的对流 换热系数。
3
研究生课程《高等工程热力学与传热学》电子教案
偏微分方程的求解
) (
)]
(22)
该差分格式称之为 Crank-Nicholson 格式，它的精度要高于显 式格式和隐式格式，其截断误差为 O[( Δt ) 2 + ( Δx ) 2 ] 。
9
研究生课程《高等工程热力学与传热学》电子教案
常见差分格式的稳定性分析
所谓稳定性问题，是指：若定义第 n 步第 j 个节点的数值求解 误差为 ε n j ，如果定义误差放大倍数 ω 满足
13
研究生课程《高等工程热力学与传热学》电子教案
常见差分方程求解
1) 显式差分方程： 直接求解。 2) 隐式差分方程或者 Crank-Nicholson 差分方程： 一维问题：追赶法； 二维问题：高斯消元法或者迭代法（如高斯－赛德尔法） 。 快速算法－交替方向法（俞昌铭，1981）
14
研究生课程《高等工程热力学与传热学》电子教案
+1 εn j ω = n ≤1 εj
(23)
则称该数值解法是稳定的，反之则为不稳定。常见的稳定性分 析方法是傅里叶级数法。一般来说，稳态问题的差分方法不存在稳 定性问题。 1) 显式格式 假设误差函数 ε ( Δx , Δt ) 可以展开成傅里叶级数的形式： ∞ ε ( Δx , Δt ) = ∑ Am cos β m Δx (24)
(7) (8) (9)
一维热传导方程的差分格式
u
x j ,tk1
t
2
2u
x j ,tk1 t 2
o( ).
(2.10)
再将 u xj1,tk1 , u xj1,tk1 分别以 x j , tk1 为中心关于 x 运用泰勒级数展开, 有
u
x j1, tk1
=u
x j , tk1
u
x j , tk1
(h） u
x xj , 0 j M ,
t tk , 0 k N
将 分割成矩形网格.记 h xj | 0 j M , tk |0 k N , h h .
称 x j , tk 为结点[1].
定义 h 上的网格函数
U
k j
|0
j
M,0 k
N
,
其中U
k j
u
xj ,tk
u
x j1, tk
=u
xj ,tk
u
xj ,tk
(h） u
xj ,tk h2 u(xj ,tk（) -h）3
2!
3!
u(4)
xj ,tk
h4
o(h4 ) ,
4!
u
x j1, tk
=u
xj ,tk
u
xj ,tk
u h
xj ,tk 2!
h2 u
xj ,tk 3!
舍去截断误差,
用
u
k j
代替
u
xj ,tk
,得到如下差分方程
u k 1 j
u
k j
a
u k 1 j 1
2u
k j
1
h2
u k 1 j 1
f
k j
1
,
1 j M 1, 1 k N.
热传导方程的差分格式讲解
热传导方程的左分格式—上机卖习报告二零一gg年五月一维抛物方程的初边值问题分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题:du d2u”(兀0) = sin兀X、0 <x <1w(0,O = z/(l，O = 0, r >0在f = 0.05,0.1和0.2时刻的数值解，并与解析解u^t) = e-7：l sm(^x)进行比较。
1差分格式形式设空间步长h = l/N，时间步长r>0, T=M T,网比r = r/h2.(1)向前差分格式向前差分格式，即Z = /C\) ‘“；=0 =心),必=吆=0其中，丿= 1,2,…,N —1,R = 1,2,…,M—l. ^r^at/h2表示网比。
(1)式可改写成如下:M*+1 = + (i-2r)Uj + rw*_! + tfj此格式为显格式。
其矩阵表达式如下：Q-2r r)r l-2r(j、r 1一2广rl吐7、厂1一2、用丿加(2)向后差分格式(1)向后差分格式，即=0=久形)上：=WN =a其中j = 12・・\N_l,k = H,M_L (2)式可改写成- rw ：［： + (l+2r )叶' -中；；=0 + 叭此种差分格式被称为隐格式。
其矩阵表达式如下：rl + 2r -r( j \ I”-r l + 2r-r l + 2r -rW.V-1-r 1 + 2广丿MJ< UN >(3) 六点对称格式六点差分格式：喟-0 _ a加:-2喟+唸;唏- 2”； +吃,—T2L戸 戸 J眄=0产久XJM=H ；=O.将(3)式改写成-g 唸;+ (1 + 时-1 昭=g 略 + (1 - 诃 * * 咯 + /其矩阵表达式如下：(1 + r -r/2<l-r r/2 ) ( j\ -r/2 l + rr/2 1-rui-r/2 l + r -r/2r/2 1-r r/2X-r 1+2匚M丿r/2 l-2r ；E >2利用MATLAB 求解问题的过程对每种差分格式依次取N = 40., r=l/1600, r=l/3200, el/6400,用 MATLAB 求解并图形比较数值解与精确解，用表格列出不同剖分时的Z?误差。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
u(: , i+1) =R * u(: , i);
end
u=[0 ; u(: , i+1) ; 0];
fori=1 : k
forj=1 : N+1
up(j , i)=exp(-pi*pi*t) * sin(pi*(j-1)*h);
end
end
x=0 : h : 1;
plot(x , u ,'b.', x , up(: , k) ,'r--');
t=0.1
0.0443229427820120
0.0442555796634600
0.0442226561119014
t=0.2
0.0265375944429049
0.0264980178171388
0.ቤተ መጻሕፍቲ ባይዱ264788945621501
3方法总结及分析
本文向前差分格式，向后差分格式及六点差分格式都是使用三对角系数矩阵，计算简单。根据matlab作，特别明显的是， 时，图像解析解波动特别大，与真解差距很大。这是因为 差分格式不稳定。根据误差比较，解本题时向后差分格式误差最小。衡量一个差分格式是否经济实用，有点多方面因素决定，应考虑不同的情况决定。
t=0.2
1.16572934157246e+136
0.00126149089942459
0.000315223617970853
向后差分格式
t=0.05
0.00483090554361983
0.00276517069125328
0.00172971295195573
t=0.1
0.00590373504606959
此格式为显格式。
其矩阵表达式如下：
（2）向后差分格式
向后差分格式，即
（2）
其中 （2）式可改写成
此种差分格式被称为隐格式。
其矩阵表达式如下：
（3）六点对称格式
六点差分格式：
（3）
将（3）式改写成
其矩阵表达式如下:
2利用MATLAB求解问题的过程
对每种差分格式依次取 ， ， ， ，用MATLAB求解并图形比较数值解与精确解，用表格列出不同剖分时的 误差。
k = t / t1;
fori = 1 : k
u(:,i+1) =inv(R) * u(:,i);
end
u=[0;u(:,i+1);0];
fori=1:k
forj=1:N+1
up(j,i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);
end
end
x=0:h:1;
plot(x,u,'b.',x,up(:,k),'r--');
legend('数值解','精确解');
err=norm(u-up(:,k));
end
六点差分格式：
function[ u , err ] =ld( t , t1 )
N = 40;
h = 1 / N;
x=[h:h:(1-h)]';
r = t1 / (h^2);
u(:,1)=sin(pi * x);%t=0时刻
R = zeros(N-1 , N-1);
fori = 2 : N-2
R(i ,i) = 1 + 2 * r;
R(i , i+1) = -r;
R(i, i-1) = -r;
end
R(1 , 1) = 1 + 2 * r;
R(1 , 2) = -r;
R(N-1 , N-1) = 1 + 2 * r;
R(N-1, N-2) = -r;
向前差分格式：
：
:
:
:
:
:
:
向后差分格式：
六点差分格式：
误差：
t=1/1600 t=1/3200 t=1/6400
向前差分格式
t=0.05
8.888396579076e+21
0.0014
0.00034640856587426
t=0.1
8.83785296707723e+59
0.0017
0.000422936658240724
fori = 2 : N-2
R(i , i) = 1 - 2 * r;
R(i , i+1) = r;
R(i, i-1) = r;
end
R(1 , 1) = 1 -2 * r;
R(1 , 2) = r;
R(N-1 , N-1) = 1 - 2 * r;
R(N-1, N-2) = r;
k = t / t1;
fori = 1 : k
u(:,i+1) =inv(R) *R1* u(:,i);
end
u=[0;u(:,i+1);0];
fori=1:k
forj=1:N+1
up(j,i)=exp(-pi*pi*t)*sin(pi*(j-1)*h);
end
end
x=0:h:1;
plot(x,u,'b.',x,up(:,k),'r--');
一维抛物方程的初边值问题
分别用向前差分格式、向后差分格式、六点对称格式，求解下列问题：
在 时刻的数值解，并与解析解 。
1差分格式形式
设空间步长 , 时间步长 ， ，网比 .
（1）向前差分格式
该问题是第二类初边值问题（混合问题），我们要求出所需次数的偏微商的函数 ，满足方程 和初始条件 ，及边值条件 。
fori = 2 : N-2
R1(i ,i) = 1 - r;
R1(i , i+1) = r/2;
R1(i, i-1) = r/2;
end
R1(1 , 1) = 1 - r;
R1(1 , 2) = r/2;
R1(N-1 , N-1) = 1 -2 * r;
R1(N-1, N-2) = r/2;
k = t / t1;
legend('数值解','精确解');
err=norm(u - up(: , k));
end
向后差分格式：
function[ u , err ] = xh( t , t1 )
N = 40;
h = 1 / N;
x=[h:h:(1-h)]';
r = t1 / (h^2);
u(:,1)=sin(pi * x);%t=0时刻
R = zeros(N-1 , N-1);
fori = 2 : N-2
R(i ,i) = 1 + r;
R(i , i+1) = -r/2;
R(i, i-1) = -r/2;
end
R(1 , 1) = 1 + r;
R(1 , 2) = -r/2;
R(N-1 , N-1) = 1 + 2 * r;
R(N-1, N-2) = -r;
legend('数值解','精确解');
err=norm(u-up(:,k));
end
已知 在相应区域光滑，并且在 与边值相容，使问题有唯一充分光滑的解。
取空间步长 ，和时间步长 ，其中 都是正整数。用两族平行直线 和 将矩形域 分割成矩形网络，网络格节点为 。以 表示网格内点集合，即位于矩形 的网点集合； 表示闭矩形 的网格集合； 是网格界点的集合。
向前差分格式，即
（1）
，
其中， 以 表示网比。（1）式可改写成如下：
附件程序
向前差分格式：
function[ u , err ] = xq( t , t1 )
% t 是时刻值；
% t1 是时间步长；
N = 40;
h = 1 / N;
x=[h : h : (1-h)]';
r = t1 / (h^2);
u(:,1)=sin(pi * x);%t=0时刻
R = zeros(N-1 , N-1);
0.00337797222842041
0.00211264169361075
t=0.2
0.00440853027126062
0.00252054496720081
0.00157579425393478
六点差分格式
t=0.05
0.0526000555873324
0.0525169033137740
0.0524758656498487