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

1.40E+00 1.20E+00 1.00E+00 8.00E-01 6.00E-01 4.00E-01 2.00E-01 0.00E+00 -2.00E-01 0 -4.00E-01 -6.00E-01
1 15 29 43 57 71 85 99 113 127 141 155 169 183 197 211 225 239 253 267 281 295 309
DO I=162,181,1 M(1,I)=1-0.05*(I-161) WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=182,321,1 M(1,I)=0 WRITE(8,*)'1',',',I,',',M(1,I) END DO !时间上从第二层开始循环，逐步积分 DO I=2,REC_T+1,1 !计算三角网格部分数值，循环部分从1到322-I DO J=1,322-I,1 M(I,J)=M(I-1,J)-(M(I-1,J+1)-M(I-1,J))/DX_DT WRITE(8,*)I,',',J,',',M(I,J) END DO !右上三角强行赋值为0 DO K=322-I,321,1 M(I,K)=0 WRITE(8,*)I,',',K,',',M(I,K) END DO END DO RETURN END SUBROUTINE SUBROUTINE FORMAT_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.CSV OPEN(UNIT=8,FILE='OUTPUT_C.CSV') !计算时间T的网格，节点数为网格数加1
t n ( in 1 i 1 ) 2x
Page 2 of 12
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 03 月 31 号
x =1 数值的 t
x =1 为例进行 A，B，C 三种格式的计算结果说明。 t A 格式，根据所的结果，我们绘制不同时刻的波形图像，选取时刻为初始时刻，第 35 个 时间区间时刻（3.5s） ，和最后的 80 区间时刻（8s） ，图像如下：
这里我们选择
初始时刻
1.20E+00 1.00E+00 8.00E-01 6.00E-01 4.00E-01 2.00E-01 0.00E+00
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院 实验人姓名：刘广 日期：2014 年 03 月 31 号 专业：理论与应用力学 参加人姓名：刘广 温度：18° C 学号：11309018 年级：2011 级
实验二：对流方程差分格式下的解 一、 实验目的 1、用数值方法计算一维对流方程在 A、B、C 三种差分格式下的解。 2、取为 0.05. 取值为 0.5，1，2。并作相关讨论。
(x, 0) 1 x -1 x 0 这里初始时刻的方程为 (x, 0) x 1 0 x 1 ，即为三角波，其中传播速度为 1，我们 (x, 0) 0 x 1 或 x 1
可以写如下 Fortran90 程序：
!主程序，根据输入选择差分格式，输入1为A格式，2为B格式，3为C格式 !其中输入DX_DT为差分比值，根据输入的比值划分时间网格，其中空间网格为0.05，区间为-8到8 !在A格式中，根据依赖区域，不能计算的左右上三角强行赋值为0 !在B格式中，根据依赖区域右上三角赋值为0 !在C格式中，根据依赖区域左上三角赋值为0 !最后输出为csv格式表格，其中第一列为时间，第二列为空间，第三列为波形数值 !画图时，需要选定同一时间点的二三两列，即可画图 PROGRAM MAIN IMPLICIT NONE !声明差分比值 REAL:: DX_DT !声明差分类型 INTEGER::TYPE_ABC WRITE(*,*)'请输入Dx/Dt，0.5，1或2' READ(*,*)DX_DT WRITE(*,*)'请选择差分格式，A格式输入1，B格式输入2，C格式输入3' READ(*,*)TYPE_ABC !判断类型调用不同的子程序 IF(TYPE_ABC==1) THEN CALL FORMAT_A(DX_DT) ELSE IF(TYPE_ABC==2) THEN CALL FORMAT_B(DX_DT) ELSE IF(TYPE_ABC==3) THEN CALL FORMAT_C(DX_DT) END IF WRITE(*,*)'请在目录下查看结果' PAUSE STOP END
Page 6 of 12
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 03 月 31 号
REC_T=160/DX_DT !确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移 ALLOCATE(M(REC_T+1,321)) DO I=1,140,1 M(1,I)=0 WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=141,161,1 M(1,I)=1+0.05*(I-161) WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=162,181,1 M(1,I)=1-0.05*(I-161) WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=182,321,1 M(1,I)=0 WRITE(8,*)'1',',',I,',',M(1,I) END DO !时间上从第二层开始循环，逐步积分 DO I=2,REC_T+1,1 !计算三角网格部分数值，循环部分从I到321 DO J=I,321,1 M(I,J)=M(I-1,J)-(M(I-1,J)-M(I-1,J-1))/DX_DT WRITE(8,*)I,',',J,',',M(I,J) END DO !左上三角强行赋值为0 DO K=1,I-1,1 M(I,K)=0 WRITE(8,*)I,',',K,',',M(I,K) END DO END DO RETURN END SUBROUTINE
学会对 MS-FORTRAN 的基本操作。 用 Fortran 编写程序计算一维对流方程在 A、B、C 三 种格式下的解。讨论各种格式下不同的 三：实验步骤: 首先 A 格式，我们对微分方程进行离散化，得出 A 格式的差分解的表达式，为
t 值的差分格式解的特点。 x
in1 in i0 ( xi )
二：实验仪器与设备: ① 装有 Fortran 的计算机 三、实验原理 1台
0, - 8 x 8 t 0 t x (x, 0) 1 x -1 x 0 (x, 0) x 1 0 x 1
(x, 0) 0
x 1 或 x 1
Page 3 of 12
中山大学工学院、理论与应用力学刘广编制
wenku.baidu.com
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院
SUBROUTINE FORMAT_A(DX_DT) IMPLICIT NONE !I,J,K为循环变量，REC_T为时间网格总数 INTEGER::I,J,K,REC_T !输入dx/dt的比值 REAL::DX_DT !声明动态矩阵，用于存储计算结果 REAL,ALLOCATABLE::M(:,:) !打开通道号8，输出OUTPUT_C.CSV OPEN(UNIT=8,FILE='OUTPUT_A.CSV') !计算时间T的网格，节点数为网格数加1 REC_T=160/DX_DT !确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移 ALLOCATE(M(REC_T+1,321)) DO I=1,140,1 M(1,I)=0 WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=141,161,1 M(1,I)=1+0.05*(I-161) WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=162,181,1 M(1,I)=1-0.05*(I-161) WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=182,321,1 M(1,I)=0 WRITE(8,*)'1',',',I,',',M(1,I) END DO !时间上从第二层开始循环，逐步积分，注意，时间层循环上限为Rec_T DO I=2,REC_T,1 !计算三角网格部分数值，循环部分从I到322-I DO J=I,322-I,1 M(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)
姓名：刘广
学号：11309018
日期：2014 年 03 月 31 号
Page 5 of 12
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 03 月 31 号
姓名：刘广
学号：11309018
日期：2014 年 03 月 31 号
Page 4 of 12
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院
END DO !左上三角强行赋值为0 DO K=1,I-1,1 M(I,K)=0 WRITE(8,*)I,',',K,',',M(I,K) END DO !右上三角强行赋值为0 DO K=322-I,321,1 M(I,K)=0 WRITE(8,*)I,',',K,',',M(I,K) END DO END DO RETURN END SUBROUTINE SUBROUTINE FORMAT_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.CSV OPEN(UNIT=8,FILE='OUTPUT_B.CSV') !计算时间T的网格，节点数为网格数加1 REC_T=160/DX_DT !确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移 ALLOCATE(M(REC_T+1,321)) DO I=1,140,1 M(1,I)=0 WRITE(8,*)'1',',',I,',',M(1,I) END DO DO I=141,161,1 M(1,I)=1+0.05*(I-161) WRITE(8,*)'1',',',I,',',M(1,I) END DO
Page 7 of 12
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验二
对流方程差分格式下的解
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 03 月 31 号
其中对于程序说明在程序注释中已经说明的很清楚，这里不再赘述。 四：实验结果 最终程序执行完毕之后会在目录下生成对应的输出文档，根据用户输入的 大小生成对应的网格，然后计算出所有网格节点的数值。 五：实验结果分析