一维对流方程在三种差分格式下的解

中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院专业：理论与应用力学年级：2011级实验人姓名：刘广参加人姓名：刘广日期：2014年03月31号温度：18°C学号：11309018实验二：对流方程差分格式下的解一、实验目的1、用数值方法计算一维对流方程在A、B、C三种差分格式下的解。2、取为0.05.取值为0.5，1，2。并作相关讨论。二：实验仪器与设备:①装有Fortran的计算机1台三、实验原理1x1x00)(x,1x010)(x,0x1-x10)(x,0t8x8-,0或xxt学会对MS-FORTRAN的基本操作。用Fortran编写程序计算一维对流方程在A、B、C三种格式下的解。讨论各种格式下不同的tx值的差分格式解的特点。三：实验步骤:首先A格式，我们对微分方程进行离散化，得出A格式的差分解的表达式，为1110()2()nnnniiiiiitxxPage2of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号这里初始时刻的方程为(x,0)1x-1x0(x,0)10x1(x,0)0x1x1x或，即为三角波，其中传播速度为1，我们可以写如下Fortran90程序：!主程序，根据输入选择差分格式，输入1为A格式，2为B格式，3为C格式!其中输入DX_DT为差分比值，根据输入的比值划分时间网格，其中空间网格为0.05，区间为-8到8!在A格式中，根据依赖区域，不能计算的左右上三角强行赋值为0!在B格式中，根据依赖区域右上三角赋值为0!在C格式中，根据依赖区域左上三角赋值为0!最后输出为csv格式表格，其中第一列为时间，第二列为空间，第三列为波形数值!画图时，需要选定同一时间点的二三两列，即可画图PROGRAMMAINIMPLICITNONE!声明差分比值REAL::DX_DT!声明差分类型INTEGER::TYPE_ABCWRITE(*,*)'请输入Dx/Dt，0.5，1或2'READ(*,*)DX_DTWRITE(*,*)'请选择差分格式，A格式输入1，B格式输入2，C格式输入3'READ(*,*)TYPE_ABC!判断类型调用不同的子程序IF(TYPE_ABC==1)THENCALLFORMAT_A(DX_DT)ELSEIF(TYPE_ABC==2)THENCALLFORMAT_B(DX_DT)ELSEIF(TYPE_ABC==3)THENCALLFORMAT_C(DX_DT)ENDIFWRITE(*,*)'请在目录下查看结果'PAUSESTOPENDPage3of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号SUBROUTINEFORMAT_A(DX_DT)IMPLICITNONE!I,J,K为循环变量，REC_T为时间网格总数INTEGER::I,J,K,REC_T!输入dx/dt的比值REAL::DX_DT!声明动态矩阵，用于存储计算结果REAL,ALLOCATABLE::M(:,:)!打开通道号8，输出OUTPUT_C.CSVOPEN(UNIT=8,FILE='OUTPUT_A.CSV')!计算时间T的网格，节点数为网格数加1REC_T=160/DX_DT!确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移ALLOCATE(M(REC_T+1,321))DOI=1,140,1M(1,I)=0WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=141,161,1M(1,I)=1+0.05*(I-161)WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=162,181,1M(1,I)=1-0.05*(I-161)WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=182,321,1M(1,I)=0WRITE(8,*)'1',',',I,',',M(1,I)ENDDO!时间上从第二层开始循环，逐步积分，注意，时间层循环上限为Rec_TDOI=2,REC_T,1!计算三角网格部分数值，循环部分从I到322-IDOJ=I,322-I,1M(I,J)=M(I-1,J)-(M(I-1,J+1)-M(I-1,J-1))/(2*DX_DT)WRITE(8,*)I,',',J,',',M(I,J)Page4of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号ENDDO!左上三角强行赋值为0DOK=1,I-1,1M(I,K)=0WRITE(8,*)I,',',K,',',M(I,K)ENDDO!右上三角强行赋值为0DOK=322-I,321,1M(I,K)=0WRITE(8,*)I,',',K,',',M(I,K)ENDDOENDDORETURNENDSUBROUTINESUBROUTINEFORMAT_B(DX_DT)!I,J,K为循环变量，REC_T为时间网格总数INTEGER::I,J,K,REC_T!输入dx/dt的比值REAL::DX_DT!声明动态矩阵，用于存储计算结果REAL,ALLOCATABLE::M(:,:)!打开通道号8，输出OUTPUT_C.CSVOPEN(UNIT=8,FILE='OUTPUT_B.CSV')!计算时间T的网格，节点数为网格数加1REC_T=160/DX_DT!确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移ALLOCATE(M(REC_T+1,321))DOI=1,140,1M(1,I)=0WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=141,161,1M(1,I)=1+0.05*(I-161)WRITE(8,*)'1',',',I,',',M(1,I)ENDDOPage5of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号DOI=162,181,1M(1,I)=1-0.05*(I-161)WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=182,321,1M(1,I)=0WRITE(8,*)'1',',',I,',',M(1,I)ENDDO!时间上从第二层开始循环，逐步积分DOI=2,REC_T+1,1!计算三角网格部分数值，循环部分从1到322-IDOJ=1,322-I,1M(I,J)=M(I-1,J)-(M(I-1,J+1)-M(I-1,J))/DX_DTWRITE(8,*)I,',',J,',',M(I,J)ENDDO!右上三角强行赋值为0DOK=322-I,321,1M(I,K)=0WRITE(8,*)I,',',K,',',M(I,K)ENDDOENDDORETURNENDSUBROUTINESUBROUTINEFORMAT_C(DX_DT)!I,J,K为循环变量，REC_T为时间网格总数INTEGER::I,J,K,REC_T!输入dx/dt的比值REAL::DX_DT!声明动态矩阵，用于存储计算结果REAL,ALLOCATABLE::M(:,:)!打开通道号8，输出OUTPUT_C.CSVOPEN(UNIT=8,FILE='OUTPUT_C.CSV')!计算时间T的网格，节点数为网格数加1Page6of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号REC_T=160/DX_DT!确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移ALLOCATE(M(REC_T+1,321))DOI=1,140,1M(1,I)=0WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=141,161,1M(1,I)=1+0.05*(I-161)WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=162,181,1M(1,I)=1-0.05*(I-161)WRITE(8,*)'1',',',I,',',M(1,I)ENDDODOI=182,321,1M(1,I)=0WRITE(8,*)'1',',',I,',',M(1,I)ENDDO!时间上从第二层开始循环，逐步积分DOI=2,REC_T+1,1!计算三角网格部分数值，循环部分从I到321DOJ=I,321,1M(I,J)=M(I-1,J)-(M(I-1,J)-M(I-1,J-1))/DX_DTWRITE(8,*)I,',',J,',',M(I,J)ENDDO!左上三角强行赋值为0DOK=1,I-1,1M(I,K)=0WRITE(8,*)I,',',K,',',M(I,K)ENDDOENDDORETURNENDSUBROUTINEPage7of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号其中对于程序说明在程序注释中已经说明的很清楚，这里不再赘述。四：实验结果最终程序执行完毕之后会在目录下生成对应的输出文档，根据用户输入的=1xt数值的大小生成对应的网格，然后计算出所有网格节点的数值。五：实验结果分析这里我们选择=1xt为例进行A，B，C三种格式的计算结果说明。A格式，根据所的结果，我们绘制不同时刻的波形图像，选取时刻为初始时刻，第35个时间区间时刻（3.5s），和最后的80区间时刻（8s），图像如下：Page8of120.00E+002.00E-014.00E-016.00E-018.00E-011.00E+001.20E+00115294357718599113127141155169183197211225239253267281295309初始时刻-6.00E-01-4.00E-01-2.00E-010.00E+002.00E-014.00E-016.00E-018.00E-011.00E+001.20E+001.40E+00050100150200250300350时间区间35中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验二对流方程差分格式下的解《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年03月31号很明显，我们可以看到，随着时间的推移，A格式明显呈现出发散状态。我们再选择B格式进行说明。B格式，根据所的结果，我们绘制不同时刻的波形图像，选取时刻为初始时刻，第35个时间区间时刻（3.5s），和最后的80区间时刻（8s），图像如下：Page9of12-6.00E+01-4.00E+01-2.00E+010.00E+002.00E+014.00E+016.00E+01050100150200250300350时间区间800.00E+002.00E-014.00E-016.00E-018.00E-011.00E+001.20E+001132537496173859710912113314515716918119320521722924125326527