扩散方程的差分解法

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

2、件：u(0,t) fdt) , u(L,t) f2(t) 一般假定函数f(x)，ft)， f2(t)满足连接条件，即f(0) (0)， f(L) f2(0)。2. 有限差分法有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。其控制方程中的导数用离散点上函数值的差商代替。差分格式可以分为显格式和隐格式。所谓显格式是指在任一结点上因变量在新是时间层 上的值可以通过之前的时间层上相邻结点变量的值显式解出来。由于这些层的变量值是已知的，当时间向前推进时， 空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。隐格式则是指任一结点上变量在新的时间层

3、的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。因此，隐格式总是伴随着求解巨大的代数方程组。隐格式的主要缺点是计算工作量大，因而不如显格式计算得快，但这只是就时间步长一样的情况而言的。 隐格式的主要优点是时间步长可以比显格式能够采用的最大步长大很多。显格式的时间步长受到稳定性条件的限制，而隐格式则几乎不受限制。3. 方程的离散3.1显格式采用时间前差及

4、第 n时间层的空间中心差，得-维扩散方程的显格式解：n 1nnujujuj 12u； u；1(4)t2(x)即n 1nnujujr(uj 12u； un 1)(5)式中，r3.2全隐格式采用时间前差及第（n+1）时间层的空间中心差，得一维扩散方程的全隐格式解:n 1nUj Ujn 1n 1n 1Uj 1 2UjUj 1t(x)24. 差分解的基本问题差分解的基本问题包括：适定性、相容性、收敛性和稳定性四个方面。4.1适定性在用差分方程作微分方程数值解时，首先，要求微分方程的问题是适定的。所谓适定性问题是指这一微分方程在一定的初始和边界条件下要有唯一解，并且在初始条件和边界条件稍有改变时，微分方

5、程的解也只是稍有偏离。从数学的角度讲，若微分方程的解存在并且是 唯一的，同时连续依赖于数据（初始条件、边界条件），则问题是适定的。对于抛物线方程，这是初值问题，要求给出初始条件，并且在区域边界上给出边界条件。在 本文的一维扩散问题中，除了给出t=0时的初值u（x,O）夕卜，应在左、右两端边界，即 x 0及x L,各给定一个边界条件，第一类边界条件U的值，或给定第二类边界条件 u n的值或给定第三类边界条件 u与u n的组合。由于抛物型方程具有扩散性质，它往往是随时间扩散的，u值将不断因扩散而衰减。在本文的一维扩散方程中，给定了初始条件，同时在区域的左、右端边界给定了边界条 件，满足适定性要求。

6、4.2相容性将一个偏微分方程用差分格式化为相应差分方程，当步长t和x趋近于零时，这个差分方程应当收敛于原微分方程，也就是说，相应的差分方程和微分方程之间的截断误差在任时刻任一网格点上均应趋近于零，这样的差分方程和微分方程才是相容的。对一维扩散方程:Lu采用显格式差分格式，令t ， xnUjLhu；ut则nUj2uc20xnuj 12u； u；!用Taylor级数展开代入上述差分方程中，则有nuLhUj仃h2124u4x24O ,h 0（7）（8）（9）2u)n2) jx0， h0时，上式化为nLhW可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该显格式为相容格式。 采用全隐格式差

7、分格式，令t ， x h，则n 1 nn 1n 1 n 1,n uj uj uj 1 2uj uj 1 cLhW20h用Taylor级数展开代入上述差分方程中，则有2驰； 与 O ,h20(12)t x当 0, h 0时，上式化为2nUUnLhUj (2) j 0(13)t x可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该全隐格式为相容格式。4.3稳定性差分法计算中所产生的误差(舍入误差，参数误差等)随时间衰减或不增大，则称离散 格式是稳定的，反之，则是不稳定的。分析差分方程稳定性有不同的方法，如矩阵方法，谐 波分析法等。下面用谐波分析法对以上两种离散格式的稳定性进行分析：(1

8、)显格式显格式的差分方程为LhU；n 1UjtnUjnc nnUj 1 2Uj Uj 102 0(x)2即U；1 (12r)U；r(U； 1 u； 1)其误差方程为n j1 (1n2r) jnnr( j 1j 1)任取一 k次谐波分量nikxj-ikjhjAkeAk e则n 1 jA eikjh(12r)ikh r ee ikh；(1 2r 2r coskh)误差放大因子为(14)(15)(16)(17)(18)2r 2rcoskh 12r(1coskh)(19)要满足稳定性条件，则要求对所有的1，须|14r | 1，即 0因此，一维扩散方程显格式的稳定性条件为(20)严且(2 )全隐格式全隐

9、格式的差分方程为LhU；n 1UjnUjtn 1Uj 1nUjn 1nn 1 n 1、(1 2r)UjUjr(Uj 1 Uj 1)其误差方程为(1n2r) jn 1r( j i任取一 k次谐波分量则则误差方程为ikx jikjhn 1 ikhn 1j ej 11 ikhep /c 、/ ikhikh、匚 n 1(1 2r) r(e e ) j12r(1coskh)(23)(24)(25)(26)误差放大因子为1丄（27）1 2r（1 coskh）要满足稳定性条件，则要求对所有的k值均有| | 1。从（28）式中可以看出，当r 0 （即0 ）时，| 1恒成立。因此，全隐格式是无条件稳定的。4.4

10、收敛性如果差分方程的解为 Uj,微分方程的解为U，若当t 0， x 0时，差分方程的解与微分方程的解之差Uj u（Xj）0（ 28）则称差分格式是收敛的。拉克斯（Lax）等价定理指出，如果问题是适定的，并且差分格式满足相容性条件，那 么差分格式的稳定性就是该格式收敛性的充分而必要的条件。由4.1、4.2和4.3的分析可知，显格式和全隐格式都满足适定性、相容性及稳定性的条 件，因而这两种格式满足收敛性要求。5. 差分方程的求解5.1初始条件及边界条件由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。在本文计算屮，取L 10m ,1。初始条件（t 0时）2x/L0 x L

11、/2u(x,0) f(x)(29)2(1x/ L)L /2 x L边界条件为u(O,t)t(t) 0t 0(30)u(L,t)f2(t)0其初始时刻（t 0）时的u分布如图1所示，x=0m处u随时间变化情况如图 2所示，x=10m 处u随时间变化情况如图 3所示。1u0.5-0.5-1 050100150200250300t /s图2 x=0处u随时间变化图1 -u0.5-0.5050100150200250300图3 x=10m处u随时间变化图t /s5.2空间步长及时间步长通过对差分格式的稳定性分析知，显格式空间步长与时间步长间应满足一定的关系，即式（21）。本文选用的步长如表 1所示，分

12、别用显格式和全隐格式进行数值计算，差分网格 如图4所示。表1时间步长及空间步长取值表时间步长t/s空间步长 x/m时间步数计算时长/s空间步数r2t/( x)0.0030.150000015001000.3丿 tAtj-1,n+1j,n+1j+1,n+1j-1,nj,nj+1,nAxX图4差分网格示意图5.2求解方法对于显格式，有u； 1 u； r（u；i 2u； u；J，若n时间层上的u已知，则可直接求解出（n +1）时间层上的u值。因此，由于给出了初始条件，显格式无需迭代，直接从t=0时刻开始，逐一计算下一个时间层上的u值，便可求解出各时间层上各空间点的u值。对于全隐格式，由于u；1与u；

13、11, u； 1有联系，不能直接求解，必须联解代数方程组。此时方程组为三对角方程组，可采用追赶法进行求解。5.3计算程序5.3.1显格式!显格式求解扩散方程!参数设置!parameter （nt=20001,nx=101） !nt:时间节点数；nx:空间节点数dime nsion 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,ntu(1,n)=0 u(

14、nx,n)=0 enddo! 初始条件 do j=1,nxx=(j-1)*dx if (x.le.sx/2) thenu(j,1)=2*x /sx elseu(j,1)=2*(1-x/sx)endif enddo! 求解差分方程 do n=2,ntdo 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,ntt=(n-1)*dt if (mod(int(t*1000),15000).eq.0) thenwrite(1,*) t=,t,sdo j=1,nxwri

15、te(1,*) (j-1)*dx,u(j,n) enddo endifenddo! 每隔 2.5m 输出随时间变化 do j=1,nxx=dx*(j-1)if (mod(int(x*10),25).eq.0) then write(2,*) x=,x,m do n=1,nt,200write(2,*) (n-1)*dt,u(j,n) enddoendifenddo end5.3.2 全隐格式! 全隐格式求解扩散方程 ! 参数设置 ！!implicit double precision (a-h,o-z)parameter (nt=20001,nx=101) !nt:时间节点数；nx:空间节点数

16、 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,ntu(1,n)=0 u(nx,n)=0 enddo! 初始条件 do j=1,nxx=(j-1)*dxif (x.le.sx/2) thenu(j,1)=2*x/sxelseu(j,1)=2*(1-x/sx)

17、endifenddo! 三对角方程组系数矩阵 a(1)=1 do 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,ntcall systri(nx,d,a,c,b,xx)do j=1,nx-1b(j)=xx(j) u(j,n)=xx(j) enddo enddo! 每隔 15s 输出沿程变化 do n=1,ntt=(n-1)*dt if (mod(int(t*1000),15000).eq.0) thenwrite(1,*) t=,t,sd

18、o j=1,nxwrite(1,*) (j-1)*dx,u(j,n) enddo endifenddo!每隔 2.5m 输出随时间变化 do j=1,nxx=dx*(j-1) if (mod(int(x*10),25).eq.0) then write(2,*) x=,x,m do n=1,nt,200write(2,*) (n-1)*dt,u(j,n) enddo endifenddoendsubroutine systri(n,d,a,c,b,x)!d：下对角线!c：上对角线!a：主对角线!b:方程右边常数项!x：方程的解! implicit double precision (a-h,o

19、-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) enddoy(1)=b(1)/p(1)do i=2,ny(i)=(b(i)-d(i-1)*y(i-1)/p(i)enddox(n)=y(n)do i=n-1,1,-1x(i)=y(i)-q(i)*x(i+1)enddoreturnend6 计算结果与分析6.1 显格式计算结果与分析6.1.1 随时间变化经计算，沿程断面 u 随时间变化如下

20、图 5-图 7 所示， 从图中可以看出， 各断面 u 值随着 时间逐渐递减，在t=60s时，u值衰减至接近于0，曲线的斜率随着时间逐渐趋缓，可见u的衰减速率逐渐减小。由初始条件可知， 初始时刻， u 分布关于 x=5m 对称， 对比 x=2.5m 及 x=7.5m 的 u 值随时间变化可知，这两个断面的 u 随时间变化关系一致，说明对称分布的初 始条件随着时间的变化，其分布仍具有对称性。u0.50x=2.5m0.300.200.100.000.4050t /s601020304050t /s600.100.400.300.200.000.60u0.5050t /s60图7x=7.5m处u随时间

21、变化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值基本相同。图8t=15s时u沿程分布图9t=30s时u沿程分布u0.0080.0060.0040.0020.000图11t=60s时u沿程分布6

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

23、s0.400.300.200.100.000.60u0.50图13 x=5.0m处u随时间变化t /s图14 x=7.5m处u随时间变化622沿程分布图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值基本相同。图15t=15s时u沿程分布图16t=30s时u沿程分布u0.0080.0060.0040.0020.000图18t=60s时u沿程分布6.3计算结果对比6.3.1随时间变化对比图19-图21给出了 x=2.5m , x=5m, x=7.5m处显格式与全隐格式计算所得u值随时间变化的对比图，从图中可以看出，显格式与全隐格式计算得到的曲线几乎重合基