一维对流方程在三种差分格式下的解_精品文档资料下载.pdf

1、2014 年 03 月 31 号 温度：18 C 学号：11309018 实验二实验二：对流方程差分格式下的解对流方程差分格式下的解 一、一、实验目的实验目的 1、用数值方法计算一维对流方程在 A、B、C 三种差分格式下的解。2、取为 0.05.取值为 0.5，1，2。并作相关讨论。二：实验仪器与设备二：实验仪器与设备:装有 Fortran 的计算机 1 台 三三、实验原理、实验原理 1 x 1 x00）（x,1x0 10）（x,0 x1-x 10）（x,0 t8x8-,0或xxt 学会对 MS-FORTRAN 的基本操作。用 Fortran 编写程序计算一维对流方程在 A、B、C 三种格式下

2、的解。讨论各种格式下不同的tx值的差分格式解的特点。三：实验步骤三：实验步骤:首先 A 格式，我们对微分方程进行离散化，得出 A 格式的差分解的表达式，为 1110（）2（）nnnniiiiiitxx Page 2 of 12 中山大学工学院、理论与应用力学刘广编制 实验编号及题目：工学院 姓名：刘广 学号：11309018 日期：2014 年 03 月 31 号 这里初始时刻的方程为（x,0）1x -1x0（x,0）1 0 x1（x,0）0 x1 x1x 或，即为三角波，其中传播速度为 1，我们可以写如下 Fortran90 程序：!主程序，根据输入选择差分格式，输入1为A格式，2为B格式，

3、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（*,

4、*）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 3 of 12 中山大学工学院、理论与应用力学刘广编制 实验编号及题目：2014 年 03 月

5、 31 号 SUBROUTINE FORMAT_A（DX_DT）IMPLICIT NONE !I,J,K为循环变量，REC_T为时间网格总数 INTEGER:I,J,K,REC_T !输入dx/dt的比值 REAL:声明动态矩阵，用于存储计算结果 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

6、=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

7、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）Page 4 of 12 中山大学工学院、理论与应用力学刘广编制 实验编号及题目：2014 年 03 月 31 号 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

8、SUBROUTINE FORMAT_B（DX_DT）!打开通道号8，输出OUTPUT_C.CSV OPEN（UNIT=8,FILE=OUTPUT_B.CSV）!确定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 5 of 12 中山大学工学院、理论与应用力学刘广编制 实验编号及题目：2014 年 0

9、3 月 31 号 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

10、,K）=0 WRITE（8,*）I,K,M（I,K）END DO END DO RETURN END SUBROUTINE SUBROUTINE FORMAT_C（DX_DT）!打开通道号8，输出OUTPUT_C.CSV OPEN（UNIT=8,FILE=OUTPUT_C.CSV）!计算时间T的网格，节点数为网格数加1 Page 6 of 12 中山大学工学院、理论与应用力学刘广编制 实验编号及题目：2014 年 03 月 31 号 REC_T=160/DX_DT !计算三角网格部分数值，循环部分从I到321 DO J=I,321,1 M（I,J）=M（I-1,J）-（M（I-1,J）-M（I

11、-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 Page 7 of 12 中山大学工学院、理论与应用力学刘广编制 实验编号及题目：2014 年 03 月 31 号 其中对于程序说明在程序注释中已经说明的很清楚，这里不再赘述。四：实验结果四：实验结果 最终程序执行完毕之后会在目录下生成对应的输出文档，根据用户输入的=1xt数值的大小生成对应的网格，然后计算出所有网格节点的数值。五：实验

12、结果分析五：实验结果分析 这里我们选择=1xt为例进行 A，B，C 三种格式的计算结果说明。A 格式，根据所的结果，我们绘制不同时刻的波形图像，选取时刻为初始时刻，第 35 个时间区间时刻（3.5s），和最后的 80 区间时刻（8s），图像如下：Page 8 of 12 0.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中山大学工学院、理论与应用力学刘广编制 实验编号及题目：2014 年 03 月 31 号 很明显，我们可以看到，随着时间的推移，A 格式明显呈现出发散状态。我们再选择 B 格式进行说明。B 格式，根据所的结果，我们绘制不同时刻的波形图像，选取