Fortran计算实例ppt课件

13 实例之一 1 建模 如流体力学方程组 要点 作些简化后确定维数 注意诸项之取舍 如热传导 辐射 重力等 激波管问题 初条 2 无量纲化 对1D情况 2维流体力学方程组 号已略去 3 确定解域 划分网格及设置边界条件 可划分为上千个格点 模拟区域 网格 均匀网格 最好用自适应网格 边条 在激波传播至边界之前有效 4 选取适当的差分格式 Lax格式或Friedrichs格式 精度 稳定条件 模型方程 对流体力学问题 其稳定条件为 5 编程及调试 programshocktubeparameter im 1001 implicitdoubleprecision a h o z dimensionx im d m u im t im dimensiond1 m u1 im t1 im d0 1 67d 5t0 3 d2r 2 78d 2v0 dsqrt r t0 L0 10gamma 1 4d0 可将3个变量一起放在一个2维数组 设定常数 d 密度 u 速度 t 温度 dx 2 d0 dble im 1 do10i 1 imx i dx dble i 1 1 continuedo20i 1 im 1 2d1 i 1 d0u1 i 0 d0t1 i 1 d0do30i im 1 2 1 imd1 i 1 d 1u1 i 0 d0t1 I 1 d0 划分网格 设定初值 10 20 30 dt 0 9d0 dxcontinuedo40i 1 imd i d1 i u i u1 i t i t1 i do50i 2 im 1d1 i 0 5d0 d i 1 d i 1 dt u i d i 1 d i 1 2 d0 dx d i u i 1 u i 1 2 d0 dx u1 i 0 5d0 u i 1 u i 1 dt u i u i 1 u i 1 2 d0 dx t i 1 t i 1 2 d0 dx t i d i d i 1 d i 1 2 d0 dx t1 i 0 5d0 t i 1 t i 1 dt u i t i 1 t i 1 2 d0 dx gamma 1 d0 t i u i 1 u i 1 2 d0 dx 初始时间步长 赋新值 999 40 主体部分 50 可考虑加入人为耗散 tim 0 d0 初始时间 d1 1 d 1 u1 1 u 1 t1 1 t 1 d1 im d im u1 im u im t1 im t im tim tim dtdt 1 d0do60i 2 im 1dt dmin1 dt dx dsqrt gamma t1 i dabs u1 i Continuewrite tim dt 由柯朗条件定下一时间步长 固定边界条件 主体部分 60 局地无量纲声速 屏幕输出供调试 if tim在某些时刻 calloutput tim d1 u1 t1 if 输出结果的次数等于tend 0 02 goto888If dt le 1 e 9 then输出 dt太小 并goto888goto999close 10 EndSubroutineoutput tim d1 u1 t1 If tim eq 0 open 10 form unformatted write 10 d u tReturnend 主体部分 如0 02 0 04 0 06 0 08 即每0 02时间间隔则将结果输出到数据文件 停止运算的无量纲时间 可采用另一计时变量 如tp 其初值为0 tp tp dt 待tp快接近0 02 dt 0 02 tp 然后输出结果 并重新令tp 0 888 实例之二 磁场扩散 磁场扩散方程 仅1D 1 建模 初条 解析解 2 无量纲化 3 确定解域 划分网格及设置边界条件 可划分为上千个格点 模拟区域 网格 均匀网格 边条 扩散至边界之前有效 4 选取适当的差分格式 模型方程 全隐格式 精度 稳定条件 恒稳 5 编程及调试 programdiffusionparameter im 1001 implicitdoubleprecision a h o z dimensionx im bb im dimensionb1 im Dimensiona im b im c im s im b0 1 67d 5t0 1 d5L0 1 d5 可将3个变量一起放在一个2维数组 设定常数 a b c系数 dx 1 dble im 1 do10i 1 imx i dx dble im 1 continuedo20i 2 imb1 i 1 d0b1 1 0 d0 划分网格 设定初值 10 20 dt x 2 x 1 2continuedo30i 1 imbb i b1 i do50i 2 im 1a i 2 d0 dx dx dtb i 1 d0c i 1 d0s i dx dx dt bb i a 1 1 d0b 1 0 d0c 1 0 d0s 1 0 d0a im 1 d0b im 0 d0c im 0 d0s im 1 d0 初始时间步长 赋新值 999 30 主体部分 50 tim 0 d0 初始时间 计算三对角矩阵系数 此过程由上依次向下的递推过程俗称 追 此过程由下依次向上的回代过程俗称 赶 主体部分 dt 1 d0do60i 2 im 1dt 5 d0 dmin1 dt dx dx 2 d0 Continuetim tim dtwrite tim dt 60 暂取步长为显格式的5倍 if tim在某些时刻 calloutput tim d1 u1 t1 if 输出结果的次数等于tend 0 02 goto888goto999close 10 EndSubroutineoutput tim d1 u1 t1 If tim eq 0 open 10 Returnend 主体