一维对流方程在ABC三种差分格式.docx
《一维对流方程在ABC三种差分格式.docx》由会员分享,可在线阅读,更多相关《一维对流方程在ABC三种差分格式.docx(11页珍藏版)》请在冰豆网上搜索。
![一维对流方程在ABC三种差分格式.docx](https://file1.bdocx.com/fileroot1/2023-1/11/c6624cf6-2d15-48f4-964c-a23d028e3564/c6624cf6-2d15-48f4-964c-a23d028e35641.gif)
一维对流方程在ABC三种差分格式
一、上机目的
用数值方法计算一维对流方程在A、B、C三种差分格式下的解。
取为0.05.
取值为0.5,1,2。
并作相关讨论。
二、实验原理
三、上机要求:
1.学会对MS-FORTRAN的基本操作。
2.用Fortran编写程序计算一维对流方程在A、B、C三种格式下的解。
3.讨论各种格式下不同的
值的差分格式解的特点。
四、实验程序
以A格式为例,对微分方程进行离散化,得出A格式的差分解的表达式:
B、C格式同理可以写出。
由此编写如下的Fortran程序。
注:
除了循环时间层的计算公式略有不同外,三个程序没有区别,因此这里只用一个主程序,并根据格式选择的不同,采用条件语句判断执行哪一部分循环体。
!
空间节点321,dx=0.05输出依次为(时间,空间,数值)
programmain
implicitnone
realdx_dt!
定义
的值
integerabc,r_t,i,j,k!
定义变量,abc为格式类型,r_t为时间网格数,其余为循环变量
real,allocatable:
:
s(:
:
)!
定义存储矩阵s
write(*,*)"输入dx_dt=0.5,1,2"
read(*,*)dx_dt
write(*,*)"选择格式,A,B,C分别输入1,2或3"
read(*,*)abc
!
根据格式选择生成相应的文件
if(abc==1)then
open(unit=8,file='out_a.csv')
elseif(abc==2)then
open(unit=8,file='out_b.csv')
elseif(abc==3)then
open(unit=8,file='out_c.csv')
endif
r_t=160/dx_dt!
计算时间网格总数
allocate(s(r_t+1,321))!
分配存储矩阵的空间
!
第一层赋初值
doi=1,140,1
s(1,i)=0
write(8,*)"1",",",i,",",s(1,i)
enddo
doi=141,161,1
s(1,i)=1+0.05*(i-161)
write(8,*)"1",",",i,",",s(1,i)
enddo
doi=162,181,1
s(1,i)=1-0.05*(i-161)
write(8,*)"1",",",i,",",s(1,i)
enddo
doi=182,321,1
s(1,i)=0
write(8,*)"1",",",i,",",s(1,i)
enddo
!
循环时间层,根据格式的选择来判断执行哪一部分
if(abc==1)then
doi=2,r_t,1
doj=i,322-i,1
s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j-1))/(dx_dt*2)
write(8,*)i,',',j,',',s(i,j)
enddo
dok=1,i-1,1!
余下部分赋值0,下同
s(i,k)=0
write(8,*)i,',',k,',',s(i,k)
enddo
dok=322-i,321,1
s(i,k)=0
write(8,*)i,',',k,',',s(i,k)
enddo
enddo
elseif(abc==2)then
doi=2,r_t+1,1
doj=1,322-i,1
s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j))/dx_dt
write(8,*)i,',',j,',',s(i,j)
enddo
dok=322-i,321,1
s(i,k)=0
write(8,*)i,',',k,',',s(i,k)
enddo
enddo
elseif(abc==3)then
doi=2,r_t+1,1
doj=i,321,1
s(i,j)=s(i-1,j)-(s(i-1,j)-s(i-1,j-1))/dx_dt
write(8,*)i,',',j,',',s(i,j)
enddo
dok=1,i-1,1
s(i,k)=0
write(8,*)i,',',k,',',s(i,k)
enddo
enddo
endif
!
完成提示
write(*,*)'数据已输出至源目录'
pause
stop
endprogram
五、实验结果及分析
程序运行后在对应目录下生成csv表格文件,根据输入的
的值不同生成对应的网格并计算各节点数值。
这里以
为例,就A、B、C三种格式的结果进行分析。
首先我们用MATLAB软件画出初始波形,以便于对比(这里三个格式没有区别)。
A格式,根据输出结果,选取不同时间节点用MATLAB软件绘制波形图像,这里从初始时刻开始每隔10单位取一点(即取t=1,11,21,…画在同一图里,可明显看出随时间变化的趋势)。
波形画出如下:
(下一页第一)
从图中看出,随着时间推移A格式呈现明显的发散状态。
B格式图像如下(下一页第二):
从图中看出,随着时间推移B格式也呈现明显的发散状态。
C格式图像如下:
对比之下明显可以看出,随着时间的推移,C格式能较好的保证收敛性。
同样的,附上
以及
时的作图结果,从图中可以看出
时三种格式都呈发散状态,而
时C格式能保持收敛性,A、B格式都呈发散状态。
首先给出
时三种格式的波形图(依次为A、B、C格式,直接以点序数0~321位横坐标):
下面是
时的波形图:
六、实验总结
经过分析,可以发现在
时三种格式都不收敛,而
和
时C格式收敛,A、B格式发散,这与理论上的结论也相一致。