上海交通大学硕士研究生课程《传热流动的数值分析大作业》
上海交通大学硕士研究生课程《传热流动的数值分析大作业》
“传热流动的数值分析”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）
《传热传质学》课程教学大纲 - 上海交通大学-机械与动力
《传热传质学》课程教学大纲课程名称：传热传质学课程代码：PO306学分/学时：3学分/51学时开课学期：秋季学期适用专业：机械工程及自动化、热能与动力工程、建筑环境与设备、核工程与科学及相关专业先修课程：流体力学、工程热力学、高等数学、大学物理后续课程：无开课单位：机械与动力工程学院一、课程性质和教学目标（需明确各教学环节对人才培养目标的贡献，专业人才培养目标中的知识、能力和素质见附表）课程性质：传热传质学是机械类专业的一门重要专业基础课，是机械、能源动力和相关专业的必修主干课。
教学目标：传热传质学是研究由温差引起的热量传递规律的科学。
本课程不仅为学生学习有关专业课程提供必要的基础理论知识，也为从事相关专业技术工作、科学研究工作及管理工作提供重要的理论基础。
（A5.2，A5.3，A5.4，B2，B3，B4，C1，C2，C4）本课程由基本概念、热传导、热对流、热辐射及应综合用五部分组成。
通过本课程教学，不仅使学生在热量传递过程的特点和规律、实际传热过程的综合分析等方面树立正确的概念，同时培养学生科学抽象、逻辑思维能力，进一步强化实践是检验理论的唯一标准的认识观。
具体来说：（1）掌握热传导、热对流和热辐射三种传热方式的基本规律、基本概念和相关能量守恒方程，并能用于实际传热问题的分析。
（B2、C1、C2）（2）初步掌握数值计算的基本分析过程、特点和实际应用能力以及商业数值分析软件。
（A5.2、B2）（3）初步掌握采用实验手段解决实际传热问题的技能，直观地认识传热过程的特点、测量传热参数的基本仪器。
（A5.2、A5.3、A5.4、B4）（4）能运用常用工质物性表、诺谟图、以及其他一些相关图表（如角系数图等）。
（A5.2）（5）初步具有综合分析实际传热问题的能力、从实际问题抽象为理论，并运用理论分析解决实际问题能力。
（B2、B3、C4）（6）强化理论来源于实践，实践是检验理论的唯一标准的认识观。
（A5.2，B4, C2）二、课程教学内容及学时分配（含实践、自学、作业、讨论等的内容及要求）1、绪论（2学时）：（B4, C2）本课程概论，并介绍热传导、热对流和热辐射的基本定义、基本计算公式、传热过程简单介绍以及热阻分析法。
[上海交通大学]上海交大传热学
![[上海交通大学]上海交大传热学](https://img.taocdn.com/s3/m/fee1f4f6227916888586d72d.png)
o
Φy
x dx
t t t
E o uE tin x( x ) y( y ) z( z) dd y x d z
Φxdx
§2-2 导热微分方程式及定解条件
c 内热源的生成热 Q gΦ dVΦ dxdydz
d 热力学能的增量 Qst Φdz ？
把Qin、Qout、Qg、Qst 带入前面的能量守恒方程
〔1〕画出等效热阻图
〔2〕推导总热阻〔等效热 阻〕、总传热系数以及总换 热量公式
〔3〕给出左侧外壁面温度 的计算公式
tf1, h1
1, d1
L1
2,
tf2, h2
d2
4,
d4
3, d3
L2
L3
13
Quick Review
1 传热学的定义和意义
2 三种传热方式及各自的特点和公式：
(1) 导热
ΦA dt
大纲要求
1. 热量传递的基本方式及传热机理。 2. 一维傅里叶定律的基本表达式及其中各物理量的定义、单位。 3. 牛顿冷却公式的基本表达式及其中各物理量的定义、单位。 4. 黑体辐射换热的四次方定律基本表达式及其中各物理量的定义、单位。 5. 传热过程及传热系数的定义及物理意义。 6. 热阻的概念，对流热阻、导热热阻的定义及基本表达式。 7. 接触热阻及污垢热阻的概念。 8. 使用串联热阻叠加的原那么和在换热计算中的应用。 9. 对流换热和传热过程的区别。表面传热系数〔对流换热系数〕和传热系数的
3 导热系数(热导率) q
- grad t
(1)物理意义：热导率的数值就是物体中单位温度梯度、单位时
间、通过单位面积的导热量 W(m C。) 热导率的数值表征物质
导热能力大小，由实验测定。 (2) 影响因素：物质的种类、材料成分、温度、湿度、压力、密度等
相变传热与流体流动数值分析作业1
相变传热与流体流动数值分析作业1学院（系）：能源与动力学院专业：能源与环境工程学生姓名：学号：指导教师：完成日期：大连理工大学Dalian University of TechnologyThe Finite Volume Method for Diffusion Problems Subjects：I.Consider the problem of source-free heat conduction in an insulated rod.The equation governing one-dimensional steady state conductive heat transfer is, where k equals 1000W/m/K. The ends are maintained at constant temperatures of 100℃ and 500℃,the cross-section area A is 0.01㎡.II.A large plate of thickness L=2cm,the thermal conductivity k=0.5 W/m/K, uniform heat generation q=1000kW/m3.The face A and B are at temperatures of 100℃and200℃respectively. The governing equation is .III.There is a cylindrical fin with uniform cross-section area A.The base is at a temperature of 100℃and the end is insulated. The fin is exposed to an ambient temperature of 20℃. The governing equation is.And m=hp/(kA)=25m-2,L=1m.Solution:// 王佳琪-作业1.cpp : 定义控制台应用程序的入口点。
西安交通大学西安交通大学《《《《数值传热学数值传热学
西安交通大学西安交通大学《《《《数值传热学数值传热学西安交通大学西安交通大学《《数值传热学数值传热学》》课程大作业20140114一. 题目（1）百叶窗翅片的二维模型如图1 所示。
在流动与换热已经进入周期性充分发展的阶段，可以取出一个翅片单元进行传热与流动阻力的分析计算。
在稳态，层流，常物性，翅片温度恒定的条件下，对于表1给定的几何尺寸，进行Re ＝10－500 范围内的数值模拟，揭示每个计算单元的平均Nu 数与阻力系数f 与Re 的关系；Nu ，f 以及Re 定为：112()Re ;;0.5pm m m dp dx L u L h L f Nu u νρλ==?=其中m u 为来流平均速度；m h 为每块条片的平均换热系数。
表1 几何参数L1/mmTp/mm Lp/mm Delta/mm /θ 30 18.6 30 1.5 25图1 百叶窗翅片二维模型图2 阶梯型逼近二. 建议建议与要求与要求1. 为便于处理流固耦合问题，计算可对图1中打阴影线的区域进行；2. 可采用图2 所示的阶梯型网格处理倾斜的翅片；3. 按照《西安交通大学学报》的论文格式撰写本报告；4. 2014年4月30号前交课程论文到东三楼204房间。
三. 参考文献[1] 陶文铨编著，数值传热学（第二版），2001，西安交通大学出版社，节11.2[2] Wang L B, Tao, W Q. Numerical analysis on heat transfer and fluid flow for arrays of non-uniform plate length aligned at angles to the flow direction. Int J Numerical Methods for Heat and Fluid Flow , 1997, 7(5,6):496[3] Gong L. Li Z Y, He Y L, Tao W Q. Discussion on numerical treatment of periodic boundary condition for temperature. Numerical Heat Transfer, Part B , 2007, 52(5):429-448。
热管传热性能检测系统及其检测评估方法
热管传热性能检测系统及其检测评估方法
摘要
计算机芯片的集成度倍增，使得芯片功率和散热量增加。利用液体工质的相变 传热的热管散热器被广泛运用于 CPU 冷却。为了满足热管大批量生产环节的传热 性能检验质量控制，结合热管的基本性质和实际工作特点，重新设计了热管传热性 能检测系统，简称为 Qgo-nogo 检测系统，它主要包括热管工件、加热系统装置、冷却 系统装置和测量系统四部分。
上海交通大学工程硕士学位论文
第一章 绪 论
芯吸收。并用径向雷诺数来判定蒸汽的径向流动的情况。当|Rer|<<1 时，蒸汽流动 中粘滞力起支配作用，速度分布曲线接近于通常的 Poiseuille 抛物线。在|Rer|较大 时，蒸发段和冷凝段的流动情况有所不同。对于热管内汽-液交界面压差与质量流之 间的关系，在热管的蒸发段，为了维持工作液体连续的蒸发，必须使 Pse 大于 Pve； 而在冷凝段为了维持工作液体连续的冷凝，必须使 Pvc 大于 Psc；Pse 是与蒸发段液体 液体表面温度 Tse 对应的蒸汽压力；Pve 是蒸发段的蒸汽压力；Pvc 是冷凝段的蒸汽压 力；Psc 是与冷凝段液体表面温度 Tsc 对应的蒸汽压力。在热管的具体条件下，不计 辐射，热流密度由对流和传导两部分组成。热管工作时具有以下特征：
图 1 CPU 散热功率逐年变化情况 Fig1 Time-line of power dissipation for CPU
1
上海交通大学工程硕士学位论文
第一章 绪 论
为了适应 CPU 性能的不断提高而对 CPU 的冷却散热手段提出新的要求，CPU 的冷却方式发生了许多变化。从最初的自然对流散热到在 CPU 芯片上安装自然对 流的散热器，到现在常规的使用风扇强迫冷却的散热器。常规的散热器的散热方式 都是采用铝制、铜制散热片外加风扇，依靠的是单相流体的强迫对流换热方法，这 些方法只能用于热流密度不大于 10W/cm2 的 CPU 散热，目前已经不能够满足 CPU 芯片稳定工作的需要，特别是随着内部散热空间的减小，已无法采用常规的散热方 式，必须采用新的技术来运用于 CPU 冷却中。由于热管具有极高的导热性、优良 的等温性、热流密度可变性、热流方向的可逆性、恒温环境的适应性等优良特点， 可以满足 CPU 对散热装置紧凑、可靠、控制灵活、高散热效率、不需要维修等要 求[2]。典型的热管散热器如图 2 所示。目前，热管技术已经广泛应用在电气设备散 热、CPU 和电子器件冷却、半导体组件以及大规模集成电路板的散热方面。
相变传热与流体流动数值分析作业5
相变传热与流体流动数值分析作业5学院（系）：能源与动力学院专业：能源与环境工程学生姓名：王佳琪学号：21110060指导教师：李维仲教授完成日期：2012.4.25大连理工大学Dalian University of TechnologyTransient convection-diffusion using QUICK differencingSubjects：Consider convection and diffusion in the one-dimensional domain sketched in figure 8.7. Calculate the transient temperature field if the initial temperature is zero everywhere and the boundary conditions are φ=0 atx=0 and ∂φ/∂x=0 at x=L. The data are L=1.5m,u=2m/s, ρ=1.0kg/m3 andΓ=0.03kg/m.s. The source distribution defined by Figure 8.8 applies at times t>0 with a=-200, b=100, x1=0.6m,x2=0.2m. Write a computer program to calculate the transient temperature distribution until it reaches a steady state using the implicit method for time integration and the Hayase et al of the QUICK scheme for the convective and diffusive terms, and compare this result with the analytical steady state solution.Solution:The following is the general programs:#include"stdafx.h"#include<iostream>#include<math.h>#include<stdlib.h>using namespace std;#define N 45#define M 200double aw[N][M],ap[N][M],ae[N][M],f[N][M],φ[N][M];Фφvoid tdma(int p){int i;double u[N],l[N],y[N];u[0]=ap[0][p];y[0]=f[0][p];for(i=1;i<N;i++){l[i-1]=aw[i][p]/u[i-1];u[i]=ap[i][p]-l[i-1]*ae[i-1][p];y[i]=f[i][p]-l[i-1]*y[i-1];}φ[N-1][p]=y[N-1]/u[N-1];for(i=N-2;i>=0;i--){φ[i][p]=(y[i]-ae[i][p]*φ[i+1][p])/u[i];}}int main(){double L=1.5,ρ=1,u=2,Γ=0.03,F,D,x,Da,Fa,dt=0.01,φa=0; double ap0[N][M],Sp[N][M];int i,j;x=L/N;D=Γ/x;F=ρ*u;Fa=F;Da=D;for(i=0;i<N;i++){φ[i][1]=0;}for(j=1;j<=M;j++){for(i=0;i<N;i++){if(i==0){ae[i][j]=-(D+Da/3-3.0/8*F);aw[i][j]=0;Sp[i][j]=-(8.0*D/3+Fa);f[i][j]=(8.0*Da/3+Fa)*φa+ap0[i][j]*φ[i][j-1]+(-x/2*200+100)*x;ap[i][j]=-aw[i][j]-ae[i][j]+ap0[i][j]-Sp[i][j];}else if(i==1){ae[i][j]=-D+3.0/8*F;aw[i][j]=-(D+F);Sp[i][j]=0;f[i][j]=ap0[i][j]*φ[i][j-1]+(-3*x/2*200+100)*x;ap[i][j]=-aw[i][j]-ae[i][j]+ap0[i][j]-Sp[i][j]-5.0/8*F;}else if(i==N-1){ae[i][j]=0;aw[i][j]=-(D+F-2.0/8*F);Sp[i][j]=0;f[i][j]=1.0/8*F*(-φ[i-2][j])+ap0[i][j]*φ[i][j-1];ap[i][j]=-aw[i][j]-ae[i][j]+ap0[i][j]-Sp[i][j]-3.0/8*F;}else if(i>1&&i<18){ae[i][j]=-D+3.0/8*F;aw[i][j]=-(D+F-1.0/8*F);Sp[i][j]=0;f[i][j]=1.0/8*F*(-φ[i-2][j])+ap0[i][j]*φ[i][j-1]+(-(x/2+i*x)*200+100)*x;}else if(i>=18&&i<24){ae[i][j]=-D+3.0/8*F;aw[i][j]=-(D+F-1.0/8*F);Sp[i][j]=0;f[i][j]=1.0/8*F*(-φ[i-2][j])+ap0[i][j]*φ[i][j-1]+((x/2+i*x)*100-80)*x;}else{ae[i][j]=-D+3.0/8*F;aw[i][j]=-(D+F-1.0/8*F);Sp[i][j]=0;f[i][j]=1.0/8*F*(-φ[i-2][j])+ap0[i][j]*φ[i][j-1];}ap0[i][j]=ρ*x/dt;tdma (j);}for(i=0;i<N;i++){cout <<"φ["<<i<<"]"<<"["<<M-1<<"]="<<φ[i][M-1]<<" ";}}}The result:Subject 1Comparison of the numerical results with the analytical solutionDiscussion:The analytical and numerical steady state solutions are compared in the picture above. As can be seen, the use of the QUICK scheme and a fine grid for spatial discretisation ensures near-perfect agreement.。
相变传热与流体流动数值分析作业3
相变传热与流体流动数值分析作业3学院（系）：能源与动力学院专业：能源与动力工程***名：***学号：********指导教师：宋永臣教授完成日期：2012.3.6大连理工大学Dalian University of TechnologyThe Finite Volume Method for Convection-DiffusionProblemsSubject:A propertyΦis transported by means of convection and diffusion through the one-dimensional domain sketched in Figure1.The governingequation is ddx (ρuΦ)=ddx(ΓdΦdx); boundary conditions are Φ0=1 atx=0 and ΦL=0 at x=L. Using five equally spaced cells for convection and diffusion calculate the distribution of Φa function of x for case:(i)Case1:u=0.1m/s;(ii)Case2:u=2.5m/s;(iii)Case3:u=2.5m/s with 20 grid nodes;The following data apply: Length L=1.0m, ρ=1.0kg/m3, Γ=0.1kg/m/s.Φ=1X=0Φ=0 x=1Solution:(I)The central differencing scheme:// 王佳琪-作业-中心差分.cpp : 定义控制台应用程序的入口点。
//#include<stdafx.h>#include<iostream>#include<math.h>#include<stdlib.h>#include<cstdlib>#include<iomanip>#include<fstream>#include<sstream>#include<string>#define N 5using namespace std;int i;double aw[N],b[N],ae[N],f[N],x[N];/*-------追赶法求解数组-------*/void tdma( ){double l[N],u[N],y[N];for(i=1;i<N;i++){u[0]=b[0];l[0]=0;l[i]=aw[i]/u[i-1];u[i]=b[i]-l[i]*ae[i-1];}y[0]=f[0];for(i=1;i<N;i++)y[i]=f[i]-l[i]*y[i-1];x[N-1]=y[N-1]/u[N-1];for(i=N-2;i>=0;i--)x[i]=(y[i]-ae[i]*x[i+1])/u[i];}void main(){void Output( );/*---------定义变量及边界条件---------*/double F,D,u,ρ,Γ,x,L,φA,φB,Sp[N];u=0.1;ρ=1;Γ=0.1;L=1;φA=1;φB=0;x=L/N; F=ρ*u;D=Γ/x;/*---------网格离散---------*/for(i=0;i<N;i++){if(i==0){aw[i]=0;ae[i]=-(D-F/2);Sp[i]=-(2*D+F);f[i]=(2*D+F)*φA;}else if(i==N-1){ae[i]=0;aw[i]=-(D+F/2);Sp[i]=-(2*D-F);f[i]=(2*D-F)*φB;}else{aw[i]=-(D+F/2);ae[i]=-(D-F/2);Sp[i]=0;f[i]=0;}b[i]=-aw[i]-ae[i]-Sp[i];}tdma();Output( );}void Output( ){/*---------后处理文件生成---------*/ostringstream name;name.str("");name << "central_.plt";ofstream out(name.str().c_str());out<< "Title = central"<<endl<< "VARIABLES = X , Y1 , Y2 \n"<<"ZONE I="<<N<<",J="<<1<<",F=POINT"<<endl;for(i = 0;i <N ;i++){out<< double(0.2*i+0.1) << " " << x[i] << " "<<(2.7183-exp(0.2*i+0.1))/1.7183 << endl;}}Results:(i)case1: u=0.1m/s; N=5;（ii）Case2: u=2.5m/s; N=5;（iii）case3:u=2.5m/s; N=20;Discussion:From the figure, we can see the agreements with the analytical solution are very good in case 1 and case 3. While in case 2, because of the finite grid with a high speed, the numerical solution appears to oscillate, called ‘wiggles’. In one word, the central differencing scheme doesn’t adapt to the high peclet number (F/D>2).(II)The upwind differencing scheme:// 王佳琪-作业3-迎风.cpp : 定义控制台应用程序的入口点。
【参考文档】上海交通大学考研大纲-范文模板 (14页)
本文部分内容来自网络整理，本司不为其真实性负责，如有异议或侵权请及时联系，本司将立即删除！== 本文为word格式，下载后可方便编辑和修改！ ==上海交通大学考研大纲篇一：上海交通大学810传热学考研大纲上海交通大学动力工程与工程热物理考试大纲与解析一、专业科目与代码：810传热学二、指定参考书《传热学》（第4版）杨世铭陶文铨高等教育出版社 201X.8《传热学重点难点及典型题精解》王秋旺西安交通大学大学出版社 201X.10三、 810动力工程与工程热物理考试大纲(传热学)与解析一、序论1.热量传递的基本方式及传热机理[conduction, convection, radiation，总结三种方式传热原理以及区别]2.一维傅立叶定律的基本表达式及其中各物理量的定义.单位。
[??=λ????公式]3.牛顿冷却公式的基本表达式及其中各物理量的定义.单位。
[??=?Δ??，h为过程量，区别状态量λ]4.黑体辐射換热的四次方定律基本表达式及其中各物理量的定义.单位。
[??=ε????4,黑度ε，玻尔兹曼常量??，热力学温度??]5.传热过程及传热系数的定义及物理意义。
[传热过程概念、区别传热过程系数和表面传热系数]6.热阻的概念，对流热阻.导热热阻的定义及基本表达式。
[???????7.接触热阻及污垢热阻的概念。
8.使用串联热阻叠加的原则和在換热计算中的应用。
[原理与电路相似]9.对流热换和传热过程的区别。
表面传热系数（对流換热系数）和传热系数的区别。
10.导热系数，表面传热系数和传热系数之间的区别。
[过程量与状态量，物性参数相同（温度压力）导热系数一定，表面传热系数和流动过程量（流动速度、状态等）有关，过程不同大小不同，不是恒定的]二、导热基本定律及稳态导热1.矢量傅立叶定律的基本表达式及其中各物理量的定义.单位。
[负号表示热流方向与温度升高方向相反]2.温度场.等温面.等温线的概念。
[等温线（面）上温度相同、区域内温度分布叫做温度场] 1δ????3.利用能量守恒定律和傅立叶定律推导导热微分方程的基本方法。
- 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。