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

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

《一维对流方程在三种差分格式下的解_精品文档资料下载.pdf》由会员分享，可在线阅读，更多相关《一维对流方程在三种差分格式下的解_精品文档资料下载.pdf（12页珍藏版）》请在冰豆网上搜索。

2014年03月31号温度：

18C学号：

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（）nnnniiiiiitxxPage2of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：

工学院姓名：

刘广学号：

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或2READ（*,*）DX_DTWRITE（*,*）请选择差分格式，A格式输入1，B格式输入2，C格式输入3READ（*,*）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中山大学工学院、理论与应用力学刘广编制实验编号及题目：

2014年03月31号SUBROUTINEFORMAT_A（DX_DT）IMPLICITNONE!

I,J,K为循环变量，REC_T为时间网格总数INTEGER:

I,J,K,REC_T!

输入dx/dt的比值REAL:

声明动态矩阵，用于存储计算结果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中山大学工学院、理论与应用力学刘广编制实验编号及题目：

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）!

打开通道号8，输出OUTPUT_C.CSVOPEN（UNIT=8,FILE=OUTPUT_B.CSV）!

确定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中山大学工学院、理论与应用力学刘广编制实验编号及题目：

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）!

打开通道号8，输出OUTPUT_C.CSVOPEN（UNIT=8,FILE=OUTPUT_C.CSV）!

计算时间T的网格，节点数为网格数加1Page6of12中山大学工学院、理论与应用力学刘广编制实验编号及题目：

2014年03月31号REC_T=160/DX_DT!

计算三角网格部分数值，循环部分从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中山大学工学院、理论与应用力学刘广编制实验编号及题目：

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中山大学工学院、理论与应用力学刘广编制实验编号及题目：

2014年03月31号很明显，我们可以看到，随着时间的推移，A格式明显呈现出发散状态。

我们再选择B格式进行说明。

B格式，根据所的结果，我们绘制不同时刻的波形图像，选取