二维对流方程在-A、-B、-C-三种差分格式

一、实验目的1、用数值方法计算二维对流方程在A、B、C三种差分格式下的解以及用时间分裂法计算C格式。2、计算范围为在−16≤x≤16及−16≤y≤16，取Δx=Δy=1，计算至少90个的时间步长。二、实验原理及步骤0yuxutu其他当,044,44,10,,yxyxu1、编程用迎风C格式、A格式在ytxt,=0.5,1,2三种情况下计算。2、计算范围为在1616,1616yx，取,0,1,1tyx计算至90个的时间步长。3、选择另外两种格式，讨论在ytxt,=0.5,1,2下的差分格式的稳定性。4、分别绘制出三种格式在不同时间t下的yxu图，并选择若干将其制作成为动画。三、实验程序首先A格式，我们对微分方程进行离散化，得出差分解的表达式为：同样的，B、C格式都可以写出。这里再给出时间分裂法C格式的差分形式：初始条件上面已经给出。由此可以编写如下的Fortran程序：programmainimplicitnonerealdt_dx!定义变量表示dt_dx和dt_dy的值integertype,r_t,t,i,j,k,l!定义变量，type确定格式类型，r_t为时间网格数，其余为循环变量real,allocatable::s(:,:,:),ts(:,:,:)!定义存储矩阵s和ts,ts用于时间分裂法write(*,*)输入dt_dx和dt_dy的值，0.5，1或者2read(*,*)dt_dxwrite(*,*)输入格式类型，A、B、C格式分别输入1、2、3，时间分裂法C格式输入4read(*,*)typeselectcase(type)!根据类型输出相应的文件case(1)open(unit=8,file='out_a.csv')case(2)open(unit=8,file='out_b.csv')case(3)open(unit=8,file='out_c.csv')case(4)open(unit=8,file='out_ch.csv')endselectr_t=90!时间网格总数为90个步长allocate(s(r_t+1,33,33))!分配存储矩阵的空间s(:,:,:)=0!时间分裂法C格式分配ts矩阵空间并赋初值if(type==4)thenallocate(ts(r_t+1,33,33))ts(:,:,:)=0endif!第一层赋初值,只要赋值不为零的部分doi=13,21,1doj=13,21,1s(1,i,j)=1enddoenddo!循环时间层dot=2,r_t,1selectcase(type)!根据类型选择不同的递推公式，边界自动赋值为零case(1)doi=2,32,1doj=2,32,1s(t,i,j)=s(t-1,i,j)-((s(t-1,i+1,j)-s(t-1,i-1,j))+(s(t-1,i,j+1)-s(t-1,i,j-1)))*(dt_dx)/2enddoenddocase(2)doi=1,32,1doj=1,32,1s(t,i,j)=s(t-1,i,j)-((s(t-1,i+1,j)-s(t-1,i,j))+(s(t-1,i,j+1)-s(t-1,i,j)))*dt_dxenddoenddocase(3)doi=2,33,1doj=2,33,1s(t,i,j)=s(t-1,i,j)-((s(t-1,i,j)-s(t-1,i-1,j))+(s(t-1,i,j)-s(t-1,i,j-1)))*dt_dxenddoenddocase(4)!为提高精度，时间分裂法分奇偶讨论if(mod(t,2)==0)thendoi=2,33,1doj=2,33,1ts(t,i,j)=s(t-1,i,j)-(s(t-1,i,j)-s(t-1,i,j-1))*dt_dxs(t,i,j)=ts(t,i,j)-(ts(t,i,j)-ts(t,i-1,j))*dt_dxenddoenddoelsedoi=2,33,1doj=2,33,1ts(t,i,j)=s(t-1,i,j)-(s(t-1,i,j)-s(t-1,i-1,j))*dt_dxs(t,i,j)=ts(t,i,j)-(ts(t,i,j)-ts(t,i,j-1))*dt_dxenddoenddoendifendselectenddo!完成提示并输出相应文件dot=1,r_t,1doi=1,33,1doj=1,33,1write(8,*)t,,,i,,,j,,,s(t,i,j)enddoenddoenddowrite(*,*)'数据已输出至源目录'pausestopendprogram四、实验结果及分析运行程序后，依据选择的格式和ΔtΔx、ΔtΔy的值不同，会在源目录下生成相应的文件，将数据导入Matlab软件后可以画出图形以及做成动画。理论上可以证明，C格式的收敛条件为𝛥𝑡𝛥𝑥+Δ𝑡Δ𝑦≤1，时间分裂法C格式收敛条件为ΔtΔx≤1,ΔtΔy≤1，其余情况均发散。限于篇幅，这里选取部分情形做出图像说明。首先我们作出初始情形的图形，这里跟格式、步长等都无关。（图放在下一页）以A格式，𝛥𝑡𝛥𝑥=𝛥𝑡𝛥𝑦=1的情况为例。除初始时刻外，再选取t=5和t=10三个节点作图对比。（图均在下一页）明显发现，随着时间推移A格式发散，B格式也有类似的结果。再看C格式，作为对比，先来看𝛥𝑡𝛥𝑥=𝛥𝑡𝛥𝑦=0.5的情形，这时C格式是收敛的，仍取t=5和t=10作图。（图片为上一页最后两张）再取𝛥𝑡𝛥𝑥=𝛥𝑡𝛥𝑦=1情形来作图，意料之中，它是发散的。（下面前两张）最后我们考察时间分裂法的C格式，其在𝛥𝑡𝛥𝑥=𝛥𝑡𝛥𝑦=1（上一页中间两张）和𝛥𝑡𝛥𝑥=𝛥𝑡𝛥𝑦=0.5（上一页后两张）时都是收敛的，作图。五、总结计算的结果收敛性和之前的理论分析完全相同。A、B格式完全发散，迎风C格式在二维问题中收敛条件进一步加强为𝛥𝑡𝛥𝑥+Δ𝑡Δ𝑦≤1，而时间分裂法的迎风C格式收敛条件放宽为𝛥𝑡𝛥𝑥,Δ𝑡Δ𝑦≤1。六、附录动画保存为gif格式上交，这里附上Matlab动画代码。fig=figures=0fort=1:90s=s+1xx=-16:1:16;yy=xx;[x,y]=meshgrid(xx,yy);z=c2(:,:,t);%此处更换不同数据文件surf(x,y,z);axis([-20,20,-20,20,0,1])light('color','r')shadinginterppause(0.1);f=getframe(fig);set(gca,'nextplot','replacechildren','visible','off');[im(:,:,1,s),map]=rgb2ind(f.cdata,256,'nodither');holdoffendimwrite(im,map,'c2.gif','DelayTime',0,'LoopCount',inf)%文件名与数据文件一致