扩散方程的差分解法

1、扩散方程的差分解法在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。1.扩散方程一维扩散方程为：（1）式中，u为因知量，为扩散系数，x为坐标，t为时间。其定解条件如下： 初始条件：（2） 边界条件：（3）一般假定函数，满足连接条件，即，。

2、2.有限差分法有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。其控制方程中的导数用离散点上函数值的差商代替。差分格式可以分为显格式和隐格式。所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。因而一个差分方程常常包括几个相邻结点

3、上的未知数，未知数的个数取决于格式的构成形式。为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。因此，隐格式总是伴随着求解巨大的代数方程组。隐格式的主要缺点是计算工作量大，因而不如显格式计算得快，但这只是就时间步长一样的情况而言的。隐格式的主要优点是时间步长可以比显格式能够采用的最大步长大很多。显格式的时间步长受到稳定性条件的限制，而隐格式则几乎不受限制。3.方程的离散3.1 显格式采用时间前差及第n时间层的空间中心差，得一维扩散方程的显格式解：（4）即（5）式中，3.2 全隐格式采用时间前差及第（n+1）时间层的空间中心差，得一维扩散方程的全隐格式解：（6）

4、4.差分解的基本问题 差分解的基本问题包括：适定性、相容性、收敛性和稳定性四个方面。4.1 适定性 在用差分方程作微分方程数值解时，首先，要求微分方程的问题是适定的。所谓适定性问题是指这一微分方程在一定的初始和边界条件下要有唯一解，并且在初始条件和边界条件稍有改变时，微分方程的解也只是稍有偏离。从数学的角度讲，若微分方程的解存在并且是唯一的，同时连续依赖于数据（初始条件、边界条件），则问题是适定的。对于抛物线方程，这是初值问题，要求给出初始条件，并且在区域边界上给出边界条件。在本文的一维扩散问题中，除了给出t=0时的初值外，应在左、右两端边界，即及，各给定一个边界条件，第一类边界条件u的值，或

5、给定第二类边界条件的值或给定第三类边界条件u与的组合。 由于抛物型方程具有扩散性质，它往往是随时间扩散的，u值将不断因扩散而衰减。 在本文的一维扩散方程中，给定了初始条件，同时在区域的左、右端边界给定了边界条件，满足适定性要求。4.2 相容性 将一个偏微分方程用差分格式化为相应差分方程，当步长和趋近于零时，这个差分方程应当收敛于原微分方程，也就是说，相应的差分方程和微分方程之间的截断误差在任一时刻任一网格点上均应趋近于零，这样的差分方程和微分方程才是相容的。 对一维扩散方程：（7）采用显格式差分格式，令，则（8）用Taylor级数展开代入上述差分方程中，则有（9）当，时，上式化为（10）可见，

6、此差分格式所构成的差分方程与原来的微分方程是相容的，故该显格式为相容格式。采用全隐格式差分格式，令，则（11）用Taylor级数展开代入上述差分方程中，则有（12）当，时，上式化为（13）可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该全隐格式为相容格式。4.3 稳定性 差分法计算中所产生的误差（舍入误差，参数误差等）随时间衰减或不增大，则称离散格式是稳定的，反之，则是不稳定的。分析差分方程稳定性有不同的方法，如矩阵方法，谐波分析法等。下面用谐波分析法对以上两种离散格式的稳定性进行分析：（1）显格式显格式的差分方程为（14）即（15）其误差方程为（16）任取一k次谐波分量（17）

7、则（18）误差放大因子为（19）要满足稳定性条件，则要求对所有的k值均有，须，即。因此，一维扩散方程显格式的稳定性条件为（20）（2）全隐格式全隐格式的差分方程为（21）即（22）其误差方程为（23）任取一k次谐波分量（24）则，（25）则误差方程为（26）误差放大因子为（27）要满足稳定性条件，则要求对所有的k值均有。从（28）式中可以看出，当（即）时，恒成立。因此，全隐格式是无条件稳定的。4.4 收敛性 如果差分方程的解为，微分方程的解为，若当，时，差分方程的解与微分方程的解之差（28）则称差分格式是收敛的。 拉克斯（Lax）等价定理指出，如果问题是适定的，并且差分格式满足相容性条件，那么

8、差分格式的稳定性就是该格式收敛性的充分而必要的条件。 由4.1、4.2和4.3的分析可知，显格式和全隐格式都满足适定性、相容性及稳定性的条件，因而这两种格式满足收敛性要求。5.差分方程的求解5.1 初始条件及边界条件由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。在本文计算中，取，。初始条件（时）（29）边界条件为（30）其初始时刻（）时的u分布如图1所示，x=0m处u随时间变化情况如图2所示，x=10m处u随时间变化情况如图3所示。图1 初始时刻u分布图图2 x=0处u随时间变化图图3 x=10m处u随时间变化图5.2 空间步长及时间步长 通过对差分格式的稳定性分

9、析知，显格式空间步长与时间步长间应满足一定的关系，即式（21）。本文选用的步长如表1所示，分别用显格式和全隐格式进行数值计算，差分网格如图4所示。表1 时间步长及空间步长取值表时间步长/s空间步长/m时间步数计算时长/s空间步数0.0030.150000015001000.3图4 差分网格示意图5.2 求解方法 对于显格式，有，若n时间层上的u已知，则可直接求解出（n+1）时间层上的u值。因此，由于给出了初始条件，显格式无需迭代，直接从t=0时刻开始，逐一计算下一个时间层上的u值，便可求解出各时间层上各空间点的u值。对于全隐格式，由于与，有联系，不能直接求解，必须联解代数方程组。此时方程组为三

10、对角方程组，可采用追赶法进行求解。5.3 计算程序 5.3.1显格式!-显格式求解扩散方程-!-参数设置-！parameter (nt=20001,nx=101) !nt:时间节点数;nx:空间节点数dimension u(nx,nt)Sx=10.0!计算长度（m）dt=0.003!时间步长（s）dx=0.1!空间步长（m）alfa=1.0!扩散系数!-!open(1,file='显格式沿程变化.txt')open(2,file='显格式随时间变化.txt')r=alfa*dt/dx*2!-边界条件-!do n=1,nt u(1,n)=0 u(nx,n)=0en

11、ddo!-!-初始条件-!do j=1,nx x=(j-1)*dx if (x.le.sx/2) then u(j,1)=2*x /sx else u(j,1)=2*(1-x/sx) endifenddo!-!-求解差分方程-!do n=2,nt do j=2,nx-1 u(j,n)=u(j,n-1)+r*(u(j-1,n-1)-2*u(j,n-1)+u(j+1,n-1) enddoenddo!-!-每隔15s输出沿程变化-!do n=1,nt t=(n-1)*dt if (mod(int(t*1000),15000).eq.0) then write(1,*) 't=',t,

12、's' do j=1,nx write(1,*) (j-1)*dx,u(j,n) enddo endifenddo!-!-每隔2.5m输出随时间变化-!do j=1,nx x=dx*(j-1) if (mod(int(x*10),25).eq.0) then write(2,*) 'x=',x,'m' do n=1,nt,200 write(2,*) (n-1)*dt,u(j,n) enddo endifenddo!-!end 5.3.2全隐格式!-全隐格式求解扩散方程-!-参数设置-！!implicit double precision (a-

13、h,o-z)parameter (nt=20001,nx=101) !nt:时间节点数;nx:空间节点数dimension u(nx,nt),a(nx),b(nx),c(nx-1),d(nx-1),xx(nx)Sx=10!计算长度（m）dt=0.003!时间步长（s）dx=0.1!空间步长（m）alfa=1.0!扩散系数!-!open(1,file='全隐格式沿程变化.txt')open(2,file='全隐格式随时间变化.txt')r=alfa*dt/dx*2!-边界条件-!do n=1,nt u(1,n)=0 u(nx,n)=0enddo!-!-初始条件-!

14、do j=1,nx x=(j-1)*dx if (x.le.sx/2) then u(j,1)=2*x/sx else u(j,1)=2*(1-x/sx) endifenddo!-!-三对角方程组系数矩阵-!a(1)=1do i=1,nx-1 d(i)=-r c(i)=-r a(i+1)=1+2*renddoc(1)=0!-!-三对角方程组常数项-!do i=1,nx b(i)=u(i,1)enddo!-!-求解差分方程-!do n=2,nt call systri(nx,d,a,c,b,xx) do j=1,nx-1 b(j)=xx(j) u(j,n)=xx(j) enddoenddo!-!

15、-每隔15s输出沿程变化-!do n=1,nt t=(n-1)*dt if (mod(int(t*1000),15000).eq.0) then write(1,*) 't=',t,'s' do j=1,nx write(1,*) (j-1)*dx,u(j,n) enddo endifenddo!-!-每隔2.5m输出随时间变化-!do j=1,nx x=dx*(j-1) if (mod(int(x*10),25).eq.0) then write(2,*) 'x=',x,'m' do n=1,nt,200 write(2,*)

16、(n-1)*dt,u(j,n) enddo endifenddo!-!end subroutine systri(n,d,a,c,b,x)! d：下对角线! c：上对角线! a：主对角线! b：方程右边常数项! x：方程的解! implicit double precision (a-h,o-z) dimension a(n),d(n-1),c(n-1) dimension p(n),q(n-1),y(n),x(n),b(n) p(1)=a(1) do i=1,n-1 q(i)=c(i)/p(i) p(i+1)=a(i+1)-d(i)*q(i) enddo y(1)=b(1)/p(1) do

17、i=2,n y(i)=(b(i)-d(i-1)*y(i-1)/p(i) enddo x(n)=y(n) do i=n-1,1,-1 x(i)=y(i)-q(i)*x(i+1) enddo return end6 计算结果与分析6.1显格式计算结果与分析 6.1.1随时间变化 经计算，沿程断面u随时间变化如下图5-图7所示，从图中可以看出，各断面u值随着时间逐渐递减，在t=60s时，u值衰减至接近于0，曲线的斜率随着时间逐渐趋缓，可见u的衰减速率逐渐减小。由初始条件可知，初始时刻，u分布关于x=5m对称，对比x=2.5m及x=7.5m的u值随时间变化可知，这两个断面的u随时间变化关系一致，说明对

18、称分布的初始条件随着时间的变化，其分布仍具有对称性。图5 x=2.5m处u随时间变化图6 x=5m处u随时间变化图7 x=7.5m处u随时间变化 6.1.2沿程分布图8-图11给出了t=15s、t=30s、t=45s及t=60s时u值的分布情况，从图中可以看出，u值的分布呈抛物线分布，在x=5m处达到最大值，同时u分布关于x=5m对称。可见，u分布仍保持初始时刻的对称关系，只是分布曲线趋于平缓。对比各时刻u的峰值，可知，u的峰值随着时间逐渐减小，至t=60s时，其峰值为0.002左右，此时最大值与最小值之差在0.002左右，可见，在t=60s时，扩散趋于稳定，即整个区域u值基本相同。图8 t=

19、15s时u沿程分布图9 t=30s时u沿程分布图10 t=45s时u沿程分布图11 t=60s时u沿程分布 6.2全隐格式计算结果与分析 6.2.1随时间变化 经计算，沿程断面u随时间变化如下图12-图14所示，从图中可以看出，各断面u值随着时间逐渐递减，在t=60s时，u值衰减至接近于0，曲线的斜率随着时间逐渐趋缓，可见u的衰减速率逐渐减小。由初始条件可知，初始时刻，u分布关于x=5m对称，对比x=2.5m及x=7.5m的u值随时间变化可知，这两个断面的u随时间变化关系一致，说明对称分布的初始条件随着时间的变化，其分布仍具有对称性。图12 x=2.5m处u随时间变化图13 x=5.0m处u随

20、时间变化图14 x=7.5m处u随时间变化 6.2.2沿程分布图15-图18给出了t=15s、t=30s、t=45s及t=60s时u值的分布情况，从图中可以看出，u值的分布近似呈抛物线分布，在x=5m处达到最大值，同时u分布关于x=5m对称。可见，u分布仍保持初始时刻的对称关系，只是分布曲线趋于平缓。对比各时刻u的峰值，可知，u的峰值随着时间逐渐减小，至t=60s时，其峰值为0.002左右，此时最大值与最小值之差在0.002左右，可见，在t=60s时，扩散趋于稳定，即整个区域u值基本相同。图15 t=15s时u沿程分布图16 t=30s时u沿程分布图17 t=45s时u沿程分布图18 t=60s时u沿程分布6.3计算结果对比 6.3.1随时间变化对比 图19-图21给出了x=2.5m，x=5m，x=7.5m处显格式与全隐格式计算所得u值随时间变化的对比图，从图中可以看出，显格式与全隐格式计算得到的曲线几乎重合基本一致，说明两者计算结