上海交通大学硕士研究生课程《传热流动的数值分析大作业》
上海交通大学硕士研究生课程《传热流动的数值分析大作业》
“传热流动的数值分析”2015年大作业1. 2维条件下的无粘、不可压缩流体通过出口和入口流过箱体，具体情况如图所示，求该箱体内的流线情况，所有单位为厘米。
（1）流线方程为：22220x yψψ∂∂+=∂∂使用Gauss-Seidel 线迭代，0.1x y ∆=∆=，误差0.005ξ=，结果输出中，包括在y=0,1,2,3,4,5 处的所有X 处对应的流函数值。
（2）设出口处纵向速度V =0，试采用PSOR 方法，0.1x y ∆=∆=，计算在x=0,1,2,3,4,5 处的所有Y 处对应的流函数值，以及不同的松弛系数和迭代次数的关系曲线（至少三个系数）。
答：（1）该问题为稳态问题，流线方程为椭圆型方程，在求解方程时，首先对方程在计算域内进行离散。
计算域为：{}(,)05,05x y x y Ω=≤≤≤≤，在离散时，0.1x y h ∆=∆=∆=，因此可以得到流线方程的差分方程为：1,,1,,1,,122220i j i j i j i j i j i j h h ψψψψψψ+-+--+-++=∆∆ （1）整理后可得：1,1,,1,1,4i j i j i j i j i j ψψψψψ+-+-+++=（2）在本题中，采用Gauss-Seidel 线迭代方法进行求解，扫描方向选为自左向右，此时有111,1,11,1,1,4n n n n i j i j i j i jn i jψψψψψ++++--+++++=（3）由于是线推进，因此在当前线方向求解时，之前的扫描线上的参数已经得到更新，所以方程可改写为：11,1,111,1,,111444n ni j i j n n n i j i j i j ψψψψψ+-++++-++-+-=（4） 其中11,n i j ψ+-看着当前迭代层中的已知变量。
从方程（4）中可以明显看出，在扫描线方向上，离散方程组的系数矩阵为三对角矩阵，因此该方程组可以采用追赶法进行求解。
《数值分析》课程教学大纲
《数值分析》课程教学大纲课程编号：07054352课程名称：数值分析英文名称：Numerical Analysis课程类型：学科基础课程要求：必修学时/学分：48/3 （讲课学时：40 上机学时：8）适用专业：计算机科学与技术；软件工程一、课程性质与任务“数值分析”是计算机科学与技术、软件工程等相关专业学生的学科基础课，也是其它理、工科专业本科生及研究生的必修或选修课。
数值分析是研究各种数学问题在计算机上通过数值运算，得到数值解答的方法和理论。
随着计算机系统能力的提高和新型数值软件的不断开发，无论在高科技领域还是在传统学科领域，数值分析的理论和方法的作用和影响巨大，是科学工作者和工程技术人员必备的基础知识和工具。
课程的任务是使学生能了解数值分析的基本概念，熟悉常用数值方法的构造原理，了解数值算法复杂性、误差与收敛性分析的基本方法，了解重要数值算法的软件实现过程，使学生系统掌握数值分析的基本概念和分析问题、解决问题的基本方法，为掌握更复杂的现代计算方法打好基础。
内容包括数值计算的基本方法、线性和非线性方程组解法、插值法、数值积分法及微分方程的数值解法。
二、课程与其他课程的联系先修课程：高等数学，线性代数，C语言程序设计，计算基础。
后续课程：人工智能，数字图像处理技术，大数据分析及应用。
三、课程教学目标1．学习使用计算机进行数值计算的基础知识和基本理论知识，能够分辨、选用合适的数值方法解决工程问题。
（支撑毕业能力要求1和2）2. 能掌握常用数值计算方法的构造原理，根据问题设计和综合运用算法设计问题解决方案。
（支撑毕业能力要求1和2）3. 能运用数值算法复杂性、误差与收敛性分析的基本方法初步进行算法分析。
4. 能用计算机语言实现典型的数值计算算法，得到实验技能的基本训练，并具有利用计算机解决常见数学问题的能力；（支撑毕业能力要求4）5．能通过查询阅读文献资料，了解数值分析的前沿和新发展动向，了解数值分析算法原理应用的典型工程领域。
数值传热_数值传热学大作业3gg
数值传热学 2009-2010 学年第一学期大作业 3
阶 梯 形 标 量 场 （ 长 宽 均 为 1 ） 的 纯 对 流 传 递 （ θ = 340 ）， 控 制 方 程 为 u ∂φ + v ∂φ = 0 ，上游的边界条件都是第一类的，即给定了φ 的分布，下游按开
∂x ∂y 口边界的方式处理。
图 1 示意图
（2）从图中可以得知：数值计算在剧烈变化区域（y=0.5 处），采用 CD、 SUD、QUICK 格式时，产生越界现象。
（3）计算过程中采用 STOIC 格式，计算效果最好。 （4）采用不同的格式时，其表达式在 CBC 线内，则会出现稳定解，否则会 出现解得越界现象，越界现象与对流稳定性不同。 （5）编程过程中，采用 SOR 低松弛迭代方法，所得结果比较理想，计算速 度远远 Gauss—Seidel 方法。在本次作业中，松弛系数取为 0.8。 （6）编程过程中，要将边值点带入循环进行计算，否则会导致计算结果出 错。进而也证明了，对流现象是具有方向性，其扰动只能沿下游传播。
要求：
1、分别利用 FUD,CD,SUD,QUICK,CLAM,EULER,MINMOD,MUSCL,OSHER,SMART, STOIC 来离散对流项，观察它们的计算结果有何不同。 2、用延迟修正进行求解。 3、写出详细的离散过程和求解方法。 4、用 TECPLOT 软件画出整个流场中φ 的分布，并画出 y = 0.5 时，φ 随 x 的 分布。 5、编程采用 C/C++或 FORTRAN 语言。 6、将源程序附于作业之后，程序中要有详细的注释，以反映出思路。 7、源程序电子版和打印版上交时间：截止 2009 年 12 月 28 日。 8、每组交一份作业，给出每位同学的贡献度。(总和为 100%)
哈工程传热大作业
传热学大作业班级：20121515学号：2012151531 姓名：张永宽第一题：如图所示，一个无限长矩形柱体，其横截面的边长分别为L 1和L 2，常物性。
该问题可视为二维稳态导热问题，边界条件如图中所示，其中L 1=0.6m ，L 2=0.4m ， T w1=60℃，T w2=20℃，λ=200W/(m ·K)。
(1) 编写程序求解二维导热方程。
(2) 绘制x =L 1/2和y =L 2/2处的温度场，并与解析解进行比较。
已知矩形内的温度场的解析解为()()()()1211w2w1sh sh sin ,L L L y L x t t y x t πππ+=。
（1）根据课本164页公式（b ）Tm ，n=（Tm+1,n+Tm-1,n+Tm,n+1+Tm,n-1）/4; 取步长为1cm 。
编出以下程序迭代求解内部个点温度。
a=zeros(41,61); %生成41*60的矩阵。
k=0:60;a(41,:)=20*sin(pi.*k/60); %矩形上边温度满足Tw2=sin(pi*x/L1).a=a+60; %使四周都为给定的边界条件。
for x=1:10000 %迭代10000次（估计能满足要求精度）。
for i=2:40for j=2:60a(i,j)=(a(i-1,j)+a(i,j-1)+a(i+1,j)+a(i,j+1))/4; %内部每一个点都为周围四个点温度和的四分之一。
end end end mesh(a)title('第一题（张永宽作请勿抄袭）','Fontsize',18) xlabel('x 轴张永宽作请勿抄袭,单位cm','Fontsize',14) ylabel('y 轴,单位cm','Fontsize',14) zlabel('t 轴,单位℃','Fontsize',14) 迭代一万次后个点温度数据：迭代法温度分布图：x 轴张永宽作请勿抄袭,单位cmy 轴,单位cmt 轴,单位℃（2）Y=L2/2时的温度曲线即把第一问中第21行数据画出图即可。
传热大作业-数值解法-清华-传热学
一维非稳态导热的数值解法一、导热问题数值解法的认识（一）背景所谓求解导热问题，就是对导热微分方程在规定的定解条件下的积分求解。
这样获得的解称为分析解。
近100年来，对大量几何形状及边界条件比较简单的问题获得了分析解。
但是，对于工程技术中遇到的许多几何形状或边界条件复杂的导热问题，由于数学上的困难目前还无法得出其分析解。
另一方面，在近几十年中，随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展十分迅速，并得到日益广泛的应用。
这些数值方法包括有限差分法、有限元法及边界元法等。
其中，有限差分法物理概念明确，实施方法简便，本次大作业即采用有限差分法。
（二）基本思想把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场，用有限个离散点上的值的集合来代替，将连续物理量场的求解问题转化为各离散点物理量的求解问题，将微分方程的求解问题转化为离散点被求物理量的代数方程的求解问题。
（三）基本步骤（1）建立控制方程及定解条件。
根据具体的物理模型，建立符合条件的导热微分方程和边界条件。
（2）区域离散化。
用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。
每一个节点都可以看成是以它为中心的一个小区域的代表，将小区域称之为元体。
（3）建立节点物理量的代数方程。
建立方法主要包括泰勒级数展开法和热平衡法。
（4）设立迭代初场。
（5）求解代数方程组。
（6）解的分析。
对于数值计算所获得的温度场及所需的一些其他物理量应作仔细分析，以获得定性或定量上的一些结论。
对于不符合实际情况的应作修正。
二、问题及求解（一）题目一厚度为0.1m 的无限大平壁，两侧均为对流换热边界条件，初始时两侧流体温度与壁内温度一致，1205f f t t t ===℃；已知两侧对流换热系数分别为h 1=11 W/m 2K 、h 2=23W/m 2K ，壁的导热系数λ=0.43W/mK ，导温系数a=0.3437×10-6 m 2/s 。
上海交通大学2010年硕士研究生入学考试传热学考研试题
上海交通大学2010年硕士研究生入学考试传热学考研试题一、简答题1、传热有哪几种方式，各有什么特点。
2、强制对流和自然对流各有哪几种流态，判别的准则数各是什么，给出具体表达式。
3、冬天，一手拿木棒，一手拿铁棒，另一端都放火上烤，铁棒的感觉到热；夏天，一手摸木板，一手摸铁板，铁板的感觉到凉，用传热原理说明原因。
4、北方窗户采用双层窗户隔热，其原因是什么，隔热条件是什么？5、画出水被加热沸腾的曲线，说明为什么1滴水滴在130度的钢板上比420度的钢板上先烧干。
6、有人说d2t/dx2=0的一维稳态导热温度分布与导热系数无关，因为里面没有λ，请说明这种看法对不对，为什么？7、沸腾传热与冷凝传热的强化原则是什么，请在改善表面结构上给出具体措施。
8、推导有效辐射与投入辐射的关系，并说明有效辐射有哪几部分组成。
二、分析题1、一厂的热电偶上标有时间常数1S，给出其具体表达市，用传热知识说明标注是否准确，为什么？2、自己划分网格，给出如图所示的点的二维，稳态，有内热源，第三类边界条件的离散方程3、冬天菜棚里的空气温度大于0度，但潮湿的地面仍然会结冰，为什么？并给出改善防冻的措施4、液-气换热，管内是液，管外是气，现在要强化，可以增加管内流速，减少管壁导热系数，在管外加肋片，请问哪个最有效，为什么？三、计算题1、图略，给了一个电熨斗，其实就是一维稳态平壁导热，左边是定加热功率Ф，右边是与环境对流（第三类边界条件），给出数学描述和边界条件，并求出温度分布2、（1）求出辐射强度（2）求出在与平面法线成0度和60度时的定向辐射强度3、一个水平管道，与外界发生自然对流，求与外界的换热量。
考查辐射传热和水平管的自然对流，已经给出了2个自然对流的NU数的经验式，自己计算进行选择（1）写出X12，X13，X23，X3，（1+2）；（2）画出辐射网络图；（3）求1与2的换热量。
其中所有条件已给，1为漫射面，2为黑体，3绝热4、一个换热器的计算问题，53根管子的液-液换热器，所有条件都给出。
上海交通大学研究生入学考试试题
d21m 第 3 题 附 图N χ χ =0 χ =N上海交通大学1997年硕士研究生入学考试试题 试题名称 传热学（含流体力学）答案必须写在答题纸上 传热学（含流体力学）1、输气管道内的空气温度t f =100℃，流速u=1/s, 用一支插入套管中的水银温度计测量空气温度 （见附图），温度计的读数是铁管底部的温度t h , 已知铁套管与输气管道连接处的温度t 0=50℃, 套管长度h=140mm,外径d=12mm ，材料的导热 系数λ=58.2w/(m 2·℃)，试问测温误差为多少度？ 已知温度计套管的过余温度分布式为)()]([0mh ch h x m ch -=θθ式中，综合参数 第1f u m λα/=，铁管与空气间的对流换热的准则式为参数为λ=3.21×10-2w/(m ·℃)，ν=23.13×10-6m 2/s. 2、 如附图所示，厚δ初始温度为t o 的大平板 一侧被突然置于∞t的流体中冷却，另一侧保持绝热，已知大平板材料的导热系数，密度和比热 分别为 λρ、c ，试导出大平板内节点 n=1,2,…N-1及边界节点n=0,N 的显式差分方程。
这里，N 表示平板的等分刻度数。
3、一辐射换热系统的加热面布置于顶部，底部为受热表面，顶部表面1和底部表面2间隔为1m ，面积均为1×1 m 2。
已知顶面的黑度ε1=0.2，t 1=727℃底面ε2=0.2，t 2=227℃。
其余四侧表面的温度及黑度均相同，为简化计算， 可将它看成整体看待，统称F3，F3是地面绝热表面，试计算1，2面之间的辐射换热量及表面 3的温度t 3,已知1，2面之间的角系数X 1，2=0.24、凝结液膜的流动和换热符合边界层的薄层性质，若把坐标X 取为重力方向（见附图），则竖壁膜状凝结换热时的边界层微分方程组可表示为：22)(yu g d dp y u u u lll∂∂++-=∂∂+∂∂μρχνχρ第 4 题 附 图y χu ∞第 6 题 附 图（1）22yt a y t t u ∂∂=∂∂+∂∂νχ式中，下角标l 表示液相，努谢尔特在作了若干简化假定后，将将上述方程组简化为：022=+g dyu d llρμ （2）022=dytd相应的边界条件为：y=0, u=0, t=t w (3)y=δ，0=dydut=t s(1) 试扼要说明，努谢尔特提出的简化假定有哪些？(2) 从（1）方程组简化成（2）方程组时，略去各项的依据（简化假定）以及边界条件（3）的依据（简化假定）分别是什么？（12分）5．一台逆流式换热器刚投入工作时的参数为t 1＇=360℃，t 1〃=300℃，t 2＇=30℃，t 2〃=200℃，G 1C 1=2500w/℃,k=800w/( m 2·℃),运行一年后发现，在G 1C 1，G 2C 2及t 1＇，t 2＇保持不便时，冷流体只能被加热到162℃。
水-乙二醇溶液在横向内肋圆管内强化传热和流动的数值分析
（7）
θ (0, r ) = θ ( p, r )
对于周期性的流动，压力分为两个部分如下式，
（8）
P( x, r ) = − ∆P( x / p) + Pl ( x, r )
（9）
等式右边第一项所涉及的是总流量，第二项所涉及的是局部流动特性，同时这两项必 须满足周期性的条件，用公式表示为，
水-乙二醇溶液在横向内肋圆管内强化传热和流动的数值分析
马军 1，赵红霞 1，韩吉田 1
（山东大学能源与动力工程学院，山东 济南 250061） 摘要：利用二维轴对称模型对在横向内肋圆管内流动的具有较大 Pr 的水-乙二醇溶液的强化传热和流动特 性进行了数值分析。 在表面等温的条件下， 利用 Fluent6.2 研究水-乙二醇溶液在圆管某段周期性流动的特性， 并且计算得到平均摩擦系数和传热系数。研究了不同的肋片高度 e 和节距 p 对强化传热的影响，结果发现 与具有小 Pr 的空气和水不同，水-乙二醇溶液没有像空气那样存在确定的 e/d 临界值，管内摩擦系数和肋高 直径比 e/d、肋高节距比 p/d 以及 Re 的值可以利用一个等式进行关联。而且不能利用计算空气的努塞尔数 时所采用的公式来确定其与 Re, e/d 和 p/d 的关系。 关键词：强化传热；水乙醇溶液；横向肋片；摩擦系数；努塞尔数
Numerical Study of Enhanced Heat Transfer and Flow of Water-glycol Mixture in Transversely Ribbed Circular Tubes
MA Jun1 ,ZHAO Hongxia1 ,HAN Jitian1
（1 School of Power and Energy, Shandong University, Jinan250061，Shandong，china）
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
“传热流动的数值分析”2015年大作业1. 2维条件下的无粘、不可压缩流体通过出口和入口流过箱体，具体情况如图所示，求该箱体内的流线情况，所有单位为厘米。
（1）流线方程为：22220x yψψ∂∂+=∂∂使用Gauss-Seidel 线迭代，0.1x y ∆=∆=，误差0.005ξ=，结果输出中，包括在y=0,1,2,3,4,5 处的所有X 处对应的流函数值。
（2）设出口处纵向速度V =0，试采用PSOR 方法，0.1x y ∆=∆=，计算在x=0,1,2,3,4,5 处的所有Y 处对应的流函数值，以及不同的松弛系数和迭代次数的关系曲线（至少三个系数）。
答：（1）该问题为稳态问题，流线方程为椭圆型方程，在求解方程时，首先对方程在计算域内进行离散。
计算域为：{}(,)05,05x y x y Ω=≤≤≤≤，在离散时，0.1x y h ∆=∆=∆=，因此可以得到流线方程的差分方程为：1,,1,,1,,122220i j i j i j i j i j i j h h ψψψψψψ+-+--+-++=∆∆ （1）整理后可得：1,1,,1,1,4i j i j i j i j i j ψψψψψ+-+-+++=（2）在本题中，采用Gauss-Seidel 线迭代方法进行求解，扫描方向选为自左向右，此时有111,1,11,1,1,4n n n n i j i j i j i jn i jψψψψψ++++--+++++=（3）由于是线推进，因此在当前线方向求解时，之前的扫描线上的参数已经得到更新，所以方程可改写为：11,1,111,1,,111444n ni j i j n n n i j i j i j ψψψψψ+-++++-++-+-=（4） 其中11,n i j ψ+-看着当前迭代层中的已知变量。
从方程（4）中可以明显看出，在扫描线方向上，离散方程组的系数矩阵为三对角矩阵，因此该方程组可以采用追赶法进行求解。
从左向右如此推进，在做完了全区内各列的求解后就完成了一轮迭代，具体的程序代码如附录1所示。
执行程序请参见文件包。
对于该方法，进出口边界的流函数假设为ax by ψ=+，对于进口，由于0,00,0.75120, 0ψψ==，因此()1600.75y ψ=-。
而对于出口边界，由于5,55,4.250, 120ψψ==，且进口流量相等，因此可以获得()1605y ψ=-。
计算输出在y=0,1,2,3,4,5 处的所有X 处对应的流函数值，具体数值请参见文件夹PSOR 中执行程序的输出文件RESULTS.TXT ，为了方便观察计算结果，本文采用Tecplot 作出云图，通过云图可以清晰的观察到计算域内的流线，如图1所示。
图中每条网格线的交点即为流函数值。
此方法的迭代次数为645次。
图1 Gauss-Seidel 线迭代方法求解结果（2）本题采用PSOR 方法进行计算。
在该方法中，需要对,k i j ψ进行预计算，本文采用Gauss-Seidel 点迭代对该值进行预计算，即：111111/L M k nnn kkl l kl l k kk l l k a a b a ψψψ⨯--==+⎛⎫=++ ⎪⎝⎭∑∑ （5）在本文中，Gauss-Seidel 点迭代实施步骤方式为先对同一条X 线逐点计算，然后再从左向右进行更新。
在程序计算时，方案实施非常简便，在同一轮迭代中，将计算域内所有的流函数值采用同一个矩阵进行存储，当计算更新某个点时，采用矩阵赋值的方法及时进行更新。
PSOR 方法的计算方程为：()11n n n k k k ψαψαψ-=-+ （6）和题（1）相似，边界条件也采用相同的，此外，出口纵向速度V=0进一步确定()ψ=-。
y1605在x=0,1,2,3,4,5 处的所有Y处对应的流函数值通过云图进行了展示，如图2所示。
此外详细的数值结果也可以在该题的程序文件夹PSOR中的RESULTS.TXT中获得，在此不做详细列出。
具体程序代码如附录2所示。
图2 PSOR方法求解结果此外，本文采用了7个不同的松弛系数进行了计算，分布是1.2,1.4,1.5,1.6,1.8,1.9,2.0,其中当松弛系数为2.0时，计算结果不再收敛，其余计算结果如图3所示。
随着松弛系数的增大，迭代次数减小。
图3 不同松弛系数下的迭代次数2. 试使用一阶波动方程，计算波在一个封闭管道内从t=0s 到t=0.15s 的传播过程。
假设声速V=200m/s ，在t=0s 的初始波形为如下图所示的三角波，请分别用以下方法和步长求解。
102030405060700510152025303540u (x )x（1） 一阶上风法 （2） Lax-Wendroff （3） BTCS 隐式法 使用以下的不同的步长： （a ） 1.0, 0.005x t ∆=∆= （b ） 1.0, 0.0025x t ∆=∆= （c ） 1.0, 0.00125x t ∆=∆=请以0.025s 为间隔图示0s 到0.15s 的波传递情况。
答：一阶波动方程为：u uV t x∂∂=-∂∂ （7）按照题目要求，对上述波动方程进行离散可以获得相应的离散方程并进行求解。
（1） 一阶上风法111n n ni i i V t V t u u u x x +-∆∆⎛⎫=-+⎪∆∆⎝⎭ （8） （2） Lax-Wendroff ()()22111112222n n n nn n n iii i i i i V t V t uu u u u u u x x++-+-∆∆=--+-+∆∆ （9） （3） BTCS 隐式法 1111122n n n n i i i i V t V t u u u u x x+++-+∆∆-++=∆∆ （10）从上面的方程中可以看出，方程（8）和（9）两种格式为显示格式，因此可以通过直接逐点求解，而方程（10）为隐式格式，需要对所有点进行联立求解，明显可以看出，该方程的系数矩阵为三对角矩阵，因此可以采用追赶法进行求解。
初值条件：()00 5.0210 5.015.0,025015.025.0025.070.0x x x u x x x x ≤≤⎧⎪-≤≤⎪=⎨-+≤≤⎪⎪≤≤⎩ ()0,00(70,0)0u u ==声速V=200.0m/s在本文中，将三种方法整理在同一个程序中，因此可以在相同程序中实现所有的功能，具体的程序代码如附录3所示。
通过计算，可以获得三种方法，三种不同步长的相关计算结果，一共9张结果，如图4-12所示。
从图中可以看出，对于不同的时间步长选择，三种离散格式都存在假扩散现象，即流动扩散。
对于一阶上风法，固定空间离散大小，减小时间步长会增强流动扩散现象，对于Lax-Wendroff 格式也是如此。
而对于BTCS 隐式法，减小时间步长会削弱流动扩散现象。
u (x )x图4 一阶上风法0.005T s ∆=的计算结果u (x )x图5 一阶上风法0.0025T s ∆=的计算结果u (x )x图6 一阶上风法0.00125T s ∆=的计算结果u (x )x图7 Lax-Wendroff 格式0.005T s∆=的计算结果u (x )x图8 Lax-Wendroff 格式0.0025T s ∆=的计算结果u (x )x图9 Lax-Wendroff 格式0.00125T s ∆=的计算结果u (x )x图10 BTCS 隐式法0.005T s ∆=的计算结果u (x )x图11 BTCS 隐式法0.0025T s ∆=的计算结果u (x )x图12 BTCS 隐式法0.00125T s ∆=的计算结果附录1 Gauss-Seidel线迭代程序C--------------------------------------------------------------------------C 主程序C--------------------------------------------------------------------------Program GSUSE MATHimplicit nonereal*8 A(1:51,1:51), B(1:51), X(1:51,1:51), Y(1:51)integer I,J,K,ITERREAL*8 M(49,4),N(49)REAL*8 E1,E2,EPSopen(6,file="RESULTS.TXT")open(7,file="tecplot.DAT")A=0.0d0B=0.0d0X=0.0d0Y=0.0D0M=0.0d0N=0.0D0E1=0.0D0E2=1.0D0EPS=0.005D0ITER=0C---------------------------------------------------------C 边界条件C---------------------------------------------------------DO I=1,51X(I,1)=120.0D0X(I,51)=0.0D0ENDDODO J=1,51X(1,J)=0.0D0X(51,J)=120.0D0ENDDOX(1,1)=120.0D0X(51,51)=0.0D0DO J=2,8X(1,J)=0.1D0*(9-J)*160.0D0 !进口边界，此处Y=a*x+(120-0)/0.75*y m/s ENDDODO J=50,44,-1X(51,J)=0.1D0*(51-J)*160.0D0 !出口边界，此处假设a=0.0d0m/sENDDOA=XDO WHILE(E2>EPS)E2=0.0D0DO I=2,50 !沿X方向线推进计算DO J=2,50IF(J.EQ.2)THENM(J-1,1)=0.0D0M(J-1,2)=1.0D0M(J-1,3)=-0.25D0M(J-1,4)=0.25*(X(I-1,J)+X(I,J-1)+X(I+1,J))ELSEIF(J.EQ.50)THENM(J-1,1)=-0.25D0M(J-1,2)=1.0D0M(J-1,3)=0.0D0M(J-1,4)=0.25*(X(I-1,J)+X(I+1,J)+X(I,J+1))ELSEM(J-1,1)=-0.25D0M(J-1,2)=1.0D0M(J-1,3)=-0.25D0M(J-1,4)=0.25*(X(I-1,J)+X(I+1,J))ENDIFENDDOCALL ZG(49,M,N)C DO K=1,49C WRITE(6,"(5F)")M(K,1),M(K,2),M(K,3),M(K,4),N(K) C ENDDOC PAUSEC------------------------------------------------------C 更新X方向的值，这样就可以线推进C------------------------------------------------------DO J=2,50X(I,J)=N(J-1)ENDDOENDDODO I=1,51DO J=1,51E1=ABS(X(I,J)-A(I,J))IF(E1>E2)THENE2=E1ELSEE2=E2ENDIFENDDOENDDOA=XITER=ITER+1WRITE(*,*)E2,ITERENDDOWRITE(6,"(51F)")XWRITE(7,*)'V ARIABLES="X","Y","F"'WRITE(7,*)"ZONE I= 51 ,J= 6 ,F=POINT"DO J=1,51IF(J.EQ.1 .OR. J.EQ.11 .OR. J.EQ.21 .OR.1 J.EQ.31 .OR.J.EQ.41 .OR. J.EQ.51) THENDO I=1,51WRITE(7,"(3F9.4)")(I-1)*0.1D0,(J-1)*0.1D0,X(I,J)ENDDOENDIFENDDOCLOSE(7)CLOSE(6)END PROGRAMC------------------------------------------------------------------------------C 追赶法C------------------------------------------------------------------------------ MODULE MATHCONTAINSSUBROUTINE ZG(N1,M,X)implicit real*8(a-h,o-z)INTEGER::N1REAL*8 MXLX(N1,4),X(N1),M(N1,4)MXLX=MDO i=1,N1IF(i==1)thenMXLX(i,2)=MXLX(i,2)MXLX(i,4)=MXLX(i,4)/MXLX(i,2)ELSEMXLX(i,2)=MXLX(i,2)-MXLX(i,1)*MXLX(i-1,3)MXLX(i,4)=(MXLX(i,4)-MXLX(i,1)*MXLX(i-1,4))/MXLX(i,2) ENDIFMXLX(i,3)=MXLX(i,3)/MXLX(i,2)ENDDODO i=N1,1,-1IF(i==N1)thenX(i)=MXLX(i,4)ELSEX(i)=MXLX(i,4)-MXLX(i,3)*X(i+1) ENDIFENDDOEND SUBROUTINEEND MODULE附录2 PSOR方法程序C---------------------------------------------------------------------C 主程序C---------------------------------------------------------------------Program PSORimplicit nonereal*8 A(1:51,1:51), B(1:51), X(1:51,1:51), Y(1:51),X1(1:51,1:51)integer I,J,K,ITERREAL*8 M(49,4),N(49)REAL*8 E1,E2,EPSREAL*8 AF,BFopen(6,file="RESULTS.TXT")open(7,file="tecplot.DAT")A=0.0d0B=0.0d0X=0.0d0Y=0.0D0M=0.0d0N=0.0D0E1=0.0D0E2=1.0D0EPS=0.005D0ITER=0AF=1.9C---------------------------------------------------------C 边界条件C---------------------------------------------------------DO I=1,51X(I,1)=120.0D0X(I,51)=0.0D0ENDDODO J=1,51X(1,J)=0.0D0X(51,J)=120.0D0ENDDOX(1,1)=120.0D0X(51,51)=0.0D0DO J=2,8X(1,J)=0.1D0*(9-J)*160.0D0 !进口边界，此处Y=a*x+(120-0)/0.75*y m/s ENDDODO J=50,44,-1X(51,J)=0.1D0*(51-J)*160.0D0 !出口边界，由于V=0，所以a=0.0ENDDOA=XX1=XDO WHILE(E2>EPS)E2=0.0D0DO I=2,50 !采用高斯赛德尔迭代高斯赛德尔迭代作为预计算值DO J=2,50BF=X(I,J)X(I,J)=0.25*(A(I,J+1)+A(I,J-1)+A(I-1,J)+A(I+1,J)) !此处用A代替X，这样可以通过在后续的设定A=X体现已考虑及时更新X的计算值X(I,J)=(1-AF)*BF+AF*X(I,J) !AF为松弛因子A=X !及时更新所计算的值，体现为高斯赛德尔迭代ENDDOENDDODO I=1,51DO J=1,51E1=ABS(X(I,J)-X1(I,J))IF(E1>E2)THENE2=E1ELSEE2=E2ENDIFENDDOENDDOX1=XA=XITER=ITER+1WRITE(*,*)E2,ITERENDDOWRITE(6,"(51F)")XWRITE(7,*)'V ARIABLES="X","Y","F"'WRITE(7,*)"ZONE I= 6 ,J= 51 ,F=POINT"DO J=1,51DO I=1,51IF(I.EQ.1 .OR. I.EQ.11 .OR. I.EQ.21 .OR.1 I.EQ.31 .OR.I.EQ.41 .OR. I.EQ.51) THENWRITE(7,"(3F9.4)")(I-1)*0.1D0,(J-1)*0.1D0,X(I,J)ENDIFENDDOENDDOCLOSE(7)CLOSE(6)END PROGRAM附录3 波动方程求解程序C----------------------------------------------------------------C 主程序C----------------------------------------------------------------PROGRAM MAINUSE MATHIMPLICIT NONEREAL*8 DT,DX,V !时间步长，空间步长，声速REAL*8 NP !方法选择，1为一阶上风法，2为Lax-Wendroff，3为BTCS隐式法REAL*8 UPRE(0:70),UNEX(0:70) !当前时层和下一时层的速度值INTEGER I,JREAL*8 TIMEREAL*8 XREAL*8 TIREAL*8 A(69,4),Y(69)TIME=0.0D0C---------------------------------------------------C 初值确定C---------------------------------------------------DO I=0,70X=DBLE(I)IF(X.LE.5.0D0)THENUPRE(I)=0.0D0ELSEIF(X.GE.5.0D0 .AND. X.LE.15.0D0)THENUPRE(I)=2.0D0*X-10.0D0ELSEIF(X.GE.15.0D0 .AND. X.LE.25.0D0)THENUPRE(I)=-2.0D0*X+50.0D0ELSEUPRE(I)=0.0D0ENDIFENDDOUNEX=UPREV=200.0D0c DT=0.005c DT=0.0025DT=0.00125DX=1.0NP=3OPEN(7,FILE="INITIAL.DAT")OPEN(8,FILE="RESULT.DAT")DO I=0,70WRITE(7,"(1F8.5,4X,1F8.5)")DBLE(I),UPRE(I)ENDDOCLOSE(7)IF(NP.EQ.1)THENWRITE(*,*)"数值格式为一阶上风法"WRITE(*,*)"DT=",DT,"DX=",DXELSEIF(NP.EQ.2)THENWRITE(*,*)"数值格式为Lax-Wendroff"WRITE(*,*)"DT=",DT,"DX=",DXELSEIF(NP.EQ.3)THENWRITE(*,*)"数值格式为BTCS隐式法"WRITE(*,*)"DT=",DT,"DX=",DXENDIFC-----------------------------------------------------TIME=TIME+DT !此处设置可以避免多算一步DO WHILE(TIME.LE.0.15D0)IF(NP.EQ.1)THEN !一阶上风法DO I=1,69UNEX(I)=(1.0D0-V*DT/DX)*UPRE(I)+V*DT/DX*UPRE(I-1) ENDDOELSEIF(NP.EQ.2)THENDO I=1,69UNEX(I)=UPRE(I)-V/2.0D0*DT/DX*(UPRE(I+1)-UPRE(I-1))+2 (V*DT/DX)**2/2*(UPRE(I+1)-2.0D0*UPRE(I)+UPRE(I-1))ENDDOELSEIF(NP.EQ.3)THENDO I=1,69IF(I.EQ.1)THENA(I,1)=0.0D0A(I,2)=1.0D0A(I,3)=V/2.0D0*DT/DXA(I,4)=V/2.0D0*DT/DX*UPRE(I-1)+UPRE(I)ELSEIF(I.EQ.69)THENA(I,1)=-V/2.0D0*DT/DXA(I,2)=1.0D0A(I,3)=0.0D0A(I,4)=-V/2.0D0*DT/DX*UPRE(I+1)+UPRE(I)ELSEA(I,1)=-V/2.0D0*DT/DXA(I,2)=1.0D0A(I,3)=V/2.0D0*DT/DXA(I,4)=UPRE(I)ENDIFENDDOCALL ZG(69,A,Y)DO I=1,69UNEX(I)=Y(I)ENDDOENDIFUPRE=UNEXTI=ABS(CEILING(TIME/0.025D0)-TIME/0.025D0)IF(TI.LE.1.0D-5)THENDO I=0,70WRITE(8,"(1F8.5,4X,1F8.5,4X,1F8.5)")TIME, DBLE(I),UNEX(I)ENDDOENDIFTIME=TIME+DTENDDOCLOSE(8)ENDPROGRAM其中，子程序ZG()见附录1。