传热学导热问题的数值解法
导热问题的数值解法
1 、重点内容：① 掌握导热问题数值解法的基本思路；
② 利用热平衡法和泰勒级数展开法建立节点的离散方程。
2 、掌握内容：数值解法的实质。
3 、了解内容：了解非稳态导热问题的两种差分格式及其稳定性。
由前述3 可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支——计算传热学（数值传热学），这些数值解法主要有以下几种：（１）有限差分法（ 2 ）有限元方法（ 3 ）边界元方法
数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。如：几何形状、边界条件复杂、物性不均、多维导热问题。
分析解法与数值解法的异同点：
相同点：根本目的是相同的，即确定① t=f(x ，y ，z) ；②。不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。§4-1 导热问题数值求解的基本思想及内节点离散方程的建立
实质
对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。
这些离散点上被求物理量值的集合称为该物理量的数值解。
2 、基本思路：数值解法的求解过程可用框图4-1 表示。
由此可见：
1 ）物理模型简化成数学模型是基础；
2 ）建立节点离散方程是关键；
3 ）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。
一数值求解的步骤
如图4-2 （a ），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下：
1 建立控制方程及定解条件
控制方程：是指描写物理问题的微分方程
针对图示的导热问题，它的控制方程（即导热微分方程）为：（a ）边界条件：x=0 时，
x=H 时，
当y=0 时，
当y=W 时，
区域离散化（确立节点）
用一系列与坐标轴平行的网格线把求解区域划分成若干个子区域，用网格线的交点作为需要确定温度值的空间位置，称为节点( 结点) ，节点的位置用该节点在两个方向上的标号m ，n 表示。
相邻两节点间的距离称步长。△x, △y 每个节点都可以看成是以它为中心的一个小区域的代表把节点代表的小区域称为元体（又叫控制容积），如图4-2(b) 。
2 建立节点物理量的代数方程（离散方程）
节点上物理量的代数方程称离散方程。其过程如下：
3 首先划分各节点的类型；
4 其次，建立节点离散方程；
5 最后，代数方程组的形成。
对节点(m,n) 的代数方程，当△x=△y 时，
有：（ b ）
6 设立迭代初场
代数方程组的求解方法有直接解法与迭代解法，传热问题的有限差分法中主要采用迭代法。采用迭代法求解时，需对被求的温度场预先设定一个解，这个解称为初场，并在求解过程中不断改进。
7 求解代数方程组
求解时遇到的问题：① 线性；② 非线性；③ 收敛性等。
如图4-2 （b ），除m=1 的左边界上各节点的温度已知外，其余(M-1)N 个节点均需建立离散方程，共有(M-1)N 个方程，则构成一个封闭的代数方程组。
1 ）线性代数方程组：代数方程一经建立，其中各项系数在整个求解过程中不再变化；
2 ）非线性代数方程组：代数方程一经建立，其中各项系数在整个求解过程中不断更新。
3 ）是否收敛判断：是指用迭代法求解代数方程是否收敛，即本次迭代计算所得之解与上一次迭代计算所得之解的偏差是否小于允许值。
关于变物性( 物性为温度的函数) 导热问题，建立的离散方程，四个邻点温度的系数不是常数，而是温度的函数。在迭代计算时，这些系数应不断更新，这是非线性问题。
解的分析
通过求解代数方程，获得物体中的温度分布，根据温度场应进一步计算通过的热流量，热应力及热变形等。因此，对于数值分析计算所得的温度场及其它物理量应作详细分析，以获得定性或定量上的结论。
二、稳态导热中位于计算区域内部的节点离散方程的建立方法
1 、基本概念
?内节点：位于计算区域内部的节点，称内节点；
?差分格式：差商中的差分可以用向前、向后、中心差分表示的格式称差分格式。
2 、基本方法
方法：① 泰勒级数展开法；② 热平衡法，以下分述之。
1 ）泰勒级数展开法
如图4-3 所示，以节点(m,n) 处的二阶偏导数为例，对节点(m+1,n) 及(m-1,n) 分别写出函数t 对(m,n) 点的泰勒级数展开式：
对(m+1,n) ：
（ a ）
对（m-1,n ）：
（ b ）
（a ）+ （ b ）得：
变形为的表示式得：
上式是用三个离散点上的值计算二阶导数的严格表达式，其中：―― 称截断误差，误差量级为，即表示未明确写出的级数余项中
的最低阶数为2 。
在数值计算时，用三个相邻节点上的值近似表示二阶导数的表达式即可，则相应的略去。于是得：（ 4 -1a ）
同理：（4-1b ）
根据导热问题的控制方程( 导热微分方程) 得：
(4-2)
若△x=△y 则有：
2) 平衡法：
其本质是傅里叶导热定律和能量守恒定律的体现。对每个元体，可用傅里叶导热定律写出其能量守恒的表达式。如图4-3 所示，元体在垂直纸面方向取单位长度，通过元体界面(w,e,n,s) 所传导的热流量可以对有关的两个节点根据傅里叶定律写出：
从节点(m-1,n) 通过界面w 传导到节点(m,n) 的热流量：
（ a ）
同理：通过界面e,n,s 传导给节点（m,n ）的热流量：
（ b ）
（ c ）
（ d ）
对元体(m,n). 根据能量守恒定律可知：（4-3 ）
其中，规定：导入元体（m,n ）的热流量为正；
导出元体（m,n ）的热流量为负。
说明：① 上述分析与推导是在笛卡儿坐标系中进行的；
② 热平衡法概念清晰，过程简捷；
③ 热平衡法与§2 — 2 建立微分方程的思路与过程一致，但不同的是前者是有限大小的元体，后者是微元体。
哈工程传热学数值计算大作业
传热学 二维稳态导热问题的数值解法 杨达文2011151419 赵树明2011151427 杨文晓2011151421 吴鸿毅2011151416
第一题： a=linspace(0,0.6,121); t1=[60+20*sin(pi*a/0.6)]; t2=repmat(60,[80 121]); s=[t1;t2]; %构造矩阵 for k=1:10000000 %理论最大迭代次数，想多大就设置多大S=s; for j=2:120 for i=2:80 S(i,j)=0.25*(S(i-1,j)+S(i+1,j)+S(i,j-1)+S(i,j+1)); end end if norm(S-s)<0.0001 break; %如果符合精度要求，提前结束迭代else s=S; end end S %输出数值解 数值解数据量太大，这里就不打印出来，只画出温度分布。 画出温度分布： figure(1) xx=linspace(0,0.6,121); yy=linspace(0.4,0,81); [x,y]=meshgrid(xx,yy); surf(x,y,S) axis([0 0.6 0 0.4 60 80]) grid on xlabel('L1') ylabel('L2') zlabel('t（温度）')
.60.66666777778L 1 L 2t （温度）
A0=[S(:,61)]; for k=1:81 B1(k)=A0(81-k+1); end B1 %x=L1/2时y方向的温度 A1=[S(41,:)] %y=L2/2时x方向的温度 x=0:0.005:0.6; y=0:0.005:0.4; A2=60+20*sin(pi*x/0.6)*((exp(pi*0.2/0.6)-exp(-pi*0.2/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6) )/2) %计算y=L2/2时x方向的解析温度 B2=60+20*sin(pi*0.3/0.6)*((exp(pi*y/0.6)-exp(-pi*y/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6))/ 2) %计算x=L1/2时y方向的解析温度 figure(2) subplot(2,2,1); plot(x,A1,'g-.',x,A2,'k:x'); %画出x=L1/2时y方向的温度场、画出x=L1/2时y方向的解析温度场曲线 xlabel('L1');ylabel('t温度'); title('y=L2/2'); legend('数值解','解析解'); subplot(2,2,2); plot(x,A1-A2); %画出具体温度场与解析温度场的差值曲线 xlabel('L1');ylabel('差值'); title('y=L2/2时，比较=数值解-解析解'); subplot(2,2,3); plot(y,B1,'g-.',y,B2,'k:x'); %画出y=L2/2时x方向的温度场、画出y=L2/2时x方向的解析温度场曲线 xlabel('L2');ylabel('t温度'); title('x=L1/2'); legend('数值解','解析解'); subplot(2,2,4); plot(y,B1-B2); %画出具体温度场与解析温度场的差值曲线 xlabel('L2');ylabel('差值'); title('x=L1/2时，比较=数值解-解析解'); y=L2/2时x方向的温度： 60 60.1635347276130 60.3269574318083 60.4901561107239 60.6530189159961 60.8154342294146 60.9772907394204 61.1384775173935 61.2988840936779 61.4584005332920 61.6169175112734 61.7743263876045 61.9305192816696 62.0853891461909 62.2388298405943 62.3907362037523 62.5410041260577 62.6895306207746 62.8362138946214 62.9809534175351 63.1236499915702 63.2642058188844 63.4025245687647 63.5385114436490 63.6720732440951 63.8031184326565 63.9315571966177 64.0573015095482 64.1802651916318 64.3003639687311 64.4175155301449 64.5316395850212 64.6426579173846 64.7504944397430 64.8550752452343 64.9563286582797 65.0541852837075
传热学数值计算大作业2014011673
数值计算大作业 一、用数值方法求解尺度为100mm×100mm 的二维矩形物体的稳态导热问题。物体的导热系数λ为1.0w/m·K。边界条件分别为： 1、上壁恒热流q=1000w/m2; 2、下壁温度t1=100℃； 3、右侧壁温度t2=0℃; 4、左侧壁与流体对流换热，流体温度tf=0℃，表面传热系数 h 分别为1w/m2·K、10 w/m2·K、100w/m2·K 和1000 w/m2·K; 要求： 1、写出问题的数学描述； 2、写出内部节点和边界节点的差分方程； 3、给出求解方法； 4、编写计算程序(自选程序语言)； 5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图； 6、就一个工况下（自选）对不同网格数下的计算结果进行讨论； 7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu 时间，迭代次数）进行讨论； 8、对4个不同表面传热系数的计算结果进行分析和讨论。 9、自选一种商业软件（fluent 、ansys 等）对问题进行分析，并与自己编程计算结果进行比较验证（一个工况）。（自选项） 1、写出问题的数学描述 设H=0.1m 微分方程 22220t t x y ??+=?? x=0，0 y=H ，0 取步长δx=0.02。已知x=0,Φ=0;x=1,Φ=1.令k=ρｕ/Γ计算结果图表： 程序及数据结果： 追赶法： #include a[i]=2+0.02*k; b[i]=4; c[i]=2-0.02*k; f[i]=0; } tdma(a,b,c,f,x); for(i=0;i 第四章 复习题 1、 试简要说明对导热问题进行有限差分数值计算的基本思想与步骤。 2、 试说明用热平衡法建立节点温度离散方程的基本思想。 3、 推导导热微分方程的步骤和过程与用热平衡法建立节点温度离散方程的过程十分相似， 为什么前者得到的是精确描述，而后者解出的确实近似解。 4、 第三类边界条件边界节点的离散那方程，也可用将第三类边界条件表达式中的一阶导数 用差分公式表示来建立。试比较这样建立起来的离散方程与用热平衡建立起来的离散方程的异同与优劣。 5．对绝热边界条件的数值处理本章采用了哪些方法？试分析比较之． 6．什么是非稳态导热问题的显示格式？什么是显示格式计算中的稳定性问题？ 7．用高斯－塞德尔迭代法求解代数方程时是否一定可以得到收敛德解？不能得出收敛的解时是否因为初场的假设不合适而造成？ 8．有人对一阶导数()()()2 21,253x t t t x t i n i n i n i n ?-+-≈ ??++ 你能否判断这一表达式是否正确，为什么？ 一般性数值计算 4-1、采用计算机进行数值计算不仅是求解偏微分方程的有力工具，而且对一些复杂的经验公式及用无穷级数表示的分析解，也常用计算机来获得数值结果。试用数值方法对Bi=0.1,1,10的三种情况计算下列特征方程的根:)6,2,1( =n n μ 3,2,1,tan == n Bi n n μμ 并用计算机查明，当2 .02≥=δτ a Fo 时用式（3-19）表示的级数的第一项代替整个级数（计 算中用前六项之和来替代）可能引起的误差。 解：Bi n n =μμtan ，不同Bi 下前六个根如下表所示： Bi μ 1 μ2 μ3 μ 4 μ 5 μ 6 0.1 0.3111 3.1731 6.2991 9.4354 12.5743 15.7143 1.0 0.8603 3.4256 6.4373 9.5293 12.6453 15.7713 10 1.4289 4.3058 7.2281 10.2003 13.2142 16.2594 Fo=0.2及0.24时计算结果的对比列于下表： Fo=0.2 δ=x Bi=0.1 Bi=1 Bi=10 第一项的值 0.94879 0.62945 0.11866 前六和的值 0.95142 0.64339 0.12248 比值 0.99724 0.97833 0.96881 Fo=0.2 0=x Bi=0.1 Bi=1 Bi=10 第一项的值 0.99662 0.96514 0.83889 前六项和的值 0.994 0.95064 0.82925 比值 1.002 1.01525 1.01163 Fo=0.24 δ=x 【2-1】一食品冷藏室由内层为19 mm 厚的松木，中层为软木层，外层为51 mm 厚的混凝土所组成。内壁面温度为-17.8 ℃,混凝土外壁面温度为29.4 ℃。松木、软木和混凝土的平均热导率分别为, 3, W/(m ·K)，要求该冷藏室的热损失为15W/m 2。求所需软木的厚度及松木和软木接触面处的温度。 解：三层平壁的导热。 1）所需软木的厚度2b 由 ∑=-=3141i i i b T T q λ 得 151 .0019.00433.0762.0051.08.174.29152+++=b 解得： m b 128.02= 2）松木和软木接触面处的温度3T 由 151 .0019 .08.17153+==T q 解得：9.153-=T ℃ 解题要点：多层平壁热传导的应用。 【2-2】为减少热损失，在外径为150 mm 的饱和蒸汽管道外加有保温层。已知保温材料的热导率λ=+ 198 T(式中T 为℃),蒸汽管外壁温度为180 ℃,要求保温层外壁温度不超过50 ℃,每米管道由于热损失而造成蒸汽冷凝的量控制在1×10-4 kg/(m ·s)以下,问保温层厚度应为多少(计算时可假定蒸汽在180 ℃下冷凝)。 解：保温层平均热导率为： )./(126.02 501801098.1103.04K m W =+??+=-λ 由于本题已知的是蒸汽管道外壁面温度，即保温层内壁面温度，故为一层导热。 由 )()(21 221r r Ln T T L Q -=λπ 得： )()(21 221r r Ln T T L Q -=πλ （1） 式中：m W L Wr L Q /9.2011 103.20191013 4=???==- 将其及其它已知数据代入式（1）得： )075 .0()50180(126.029.2012r Ln -??=π 解得：m r 125.02= mm m 5005.0075.0125.0==-=∴δ壁厚 解题要点：单层圆筒壁热传导的应用。 【2-8】烤炉内在烤一块面包。已知炉壁温度为175 ℃，面包表面的黑度为,表面温度为100 ℃,表面积为 5 m 2，炉壁表面积远远大于面包表面积。求烤炉向这块面包辐射 传递的热量。 解：两物体构成封闭空间，且21S S <<，由下式计算辐射传热量： W T T S Q 0.65)448373(0645.085.01067.5) (448424111012-=-????=-=-εσ 负号表示炉壁向面包传递热量。 解题要点：辐射传热的应用，两个灰体构成的封闭空间。 【2-10】在逆流换热器中，用初温为20 ℃的水将1.25 kg/s 的液体［比热容为 kJ/(kg ·K)、密度为850 kg/m 3 ］由80 ℃冷却到30 ℃。换热器的列管直径为Φ25 mm ×2.5 mm,水走管内。水侧和液体侧的对流传热系数分别为850 W/（m 2·K ）和1 700W/（m 2·K ）,污垢热阻可忽略。若水的出口温度不能高于50 ℃，求水的流量和换热器的传热面积。 传热学大作业报告二维稳态计算 院系：能源与环境学院 专业：核工程与核技术 姓名：杨予琪 学号：03311507 一、原始题目及要求 计算要求： 1. 写出各未知温度节点的代数方程 2. 分别给出G-S 迭代和Jacobi 迭代程序 3. 程序中给出两种自动判定收敛的方法 4. 考察三种不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m ℃）） 6. 绘出最终结果的等值线 报告要求： 1. 原始题目及要求 2. 各节点的离散化的代数方程 3. 源程序 4. 不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m ℃）） 6. 计算结果的等温线图 7. 计算小结 二、各节点的离散化的代数方程 左上角节点 )(21 1,22,11,1t t t += 右上角节点 )(2 15,24,15,1t t t += 左下角节点 C t ?=1001,5 右下角节点 )2(211,24,55,5λ λ x h t t x h t ?++?+= 左边界节点 C t i ?=1001,，42≤≤i 上边界节点 C t j ?=200,1，42≤≤j 右边界节点 )2(415,15,14,5,+-++= i i i i t t t t ，42≤≤i 下边界节点 )42()2(211,51,5,4,5∞+-?+++?+=t x h t t t x h t j j j j λλ ，42≤≤j 内部节点 )(2 1,1,11,1,,j i j i j i j i j i t t t t t +-+-+++= ，4,2≤≤j i 三、源程序 1、G-S 迭代法 t=zeros(5,5); t0=zeros(5,5); dteps=0.0001; for i=2:5 %左边界节点 t(i,1)=100; end for j=2:4 %上边界节点 t(1,j)=200; end t(1,1)=(t(1,2)+t(2,1))/2; t for k=1:100 for i=2:4 %内部节点 for j=2:4 t(i,j)=(t(i-1,j)+t(i+1,j)+t(i,j-1)+t(i,j+1))/4; end end t(1,5)=(t(1,4)+t(2,5))/2;%右上角节点 for i=2:4;%右边界节点 t(i,5)=(2*t(i,4)+t(i-1,5)+t(i+1,5))/4; end for j=2:4; %下边界节点 二维导热物体温度场的数值模拟 一、物理问题 有一个用砖砌成的长方形截面的冷空气通道， 于纸面方向上用冷空气及砖墙的温度变化很小， 可以近似地予以忽略。 在下列两种情况下试计算： 砖墙横截面上的温度分布；垂直于纸面方向的每 米长度上通过砖墙的导热量。 第一种情况：内外壁分别均匀维持在 0℃及 30℃； 第二种情况：内外壁均为第三类边界条 件， 且已知： t 1 30 C,h 1 10.35W / m 2 K 2 t 2 10 C, h 2 3.93W / m 2 K 砖墙导热系数 0.35/ m K 二、数学描写 由对称的界面必是绝热面， 态、无内热源的导热问题。 控制方程： 22 tt 22 xy 边界条件： 第一种情况： 由对称性知边界 1 绝热： 边界 2 为等温边界，满足第一类边界条件： t w 0 C ； 边界 3 为等温边界，满足第一类边界条件： t w 30 C 。 第一种情况： 由对称性知边界 1 绝热： q w 0； 边界 2 为对流边界，满足第三类边界条件： q w ( t )w h 2(t w 可取左上方的四分之一墙角为研究对象， 该问题为二维、 稳 图1- t f )； n t 边界3 为对流边界，满足第三类边界条件：q w ( ) w h 2 (t w t f )。 w n w 2 w f 0,m 6,n 1~ 7;m 7 ~ 16,n 7 30,m 1,n 1~12;m 2 ~ 16,n 12 三、方程离散 用一系列与坐标轴平行的间隔 0.1m 的二维网格线 将温度区域划分为若干子区域，如图 1-3 所示。 采用热平衡法， 利用傅里叶导热定律和能量守恒定 律，按照以导入元体（ m,n ）方向的热流量为正，列写 每个节点代表的元体的代数方程， 第一种情况： 边界点： 1 边界 绝热边界） ： 边界 图1-3 t m ,1 t 16,n 等温内边界） ： 14 (2t m,2 1 4 (2t 15,n t m 1,1 t m 1,1)，m 2 ~ 5 t 16,n 1 t 16,n 1)， n 8 ~ 11 边界 等温外边界） ： 内节 点： 1 (t t t t ) 4 m 1,n m 1,n m ,n 1 m,n 1 m 2 ~ 5,n 2 ~11;m 6 ~ 15,n 8 ~ 11 t m,n 第二种情况 边界点： 边界 1（绝热边界） ： t m ,1 1 4 (2t m,2 t m 1,1 t m 1,1)，m 2 ~ 5 t 16,n 1 4 (2t 15,n t 16,n 1 t 16,n 1)， n 8 ~11 4 边界 2（内对流边界） ： t6,n 2t 5,n t 6,n 1 t 6,n 1 2Bi 1t 1 ，n 1~ 6 6,n 2(Bi 2) t m,n t m,n 1、Jacobi 迭代 在Jacobi 迭代法中任一点上未知值的更新是用上一轮迭代中所获得的各邻 点之值来计算的，即 kk k k l l n l k n k a b T a T /)(1)1()(+=∑≠=- k=1,2,...,L 1×M 1 这里带括号的上角标表示迭代轮数。所谓一轮是指把求解区域中每一节点之值都更新一次的运算环节。显然，采用Jacobi 迭代式，迭代前进的方向（又称扫描方向）并不影响迭代收敛速度。这种迭代法收敛速度很慢，一般较少采用。但对强烈的非线性问题，如果两个层次的迭代之间未知量的变化过大，容易引起非线性问题迭代的发散。在规定每一层次计算的迭代轮次数的情况下，有利于Jacobi 迭代有利于非线性问题迭代的收敛。 2、Gauss-Seidel 迭代 在这种迭代法中，每一种计算总是取邻点的最新值来进行。如果每一轮迭代按T 的下角标由小到大的方式进行，则可表示为： kk k M L k l n l kl k l l n l kl n k a b T a T a T /)(1 11 ) 1(1 1) ()(++ =∑∑?+=--≠= 此时迭代计算进行的方向（即扫描方向）会影响到收敛速度，这是与边界条件的影响传入到区域内部的快慢有关的。 3、例题： 一矩形薄板几何尺寸如图所示，薄板左侧的边界温度T L =100K ，右侧温度T R =300K ，上侧温度T T =200K ，下侧温度T B =200K ，其余各面绝热，求板上个节点的温度。要求节点数目可以变化，写出程序。 解析： ⑴列出描述问题的微分方程和定解条件。 22 220t t x y ??+=??；对于离散化的问题，其微分方程根据热平衡原理得到： 四、非线笥问题迭代式解法的收敛性 每一层次上满足迭代法求解的收敛条件+相邻次间代数方程的系数变化不太大（亦即未知量的变化不太大←多数情形下非线性问题迭代式解法是可以收敛的）。 使相邻两层次间未知量变化不太大的措施： 1、欠松弛迭代 常用逐次欠弛线迭法（SLUR ）：一组临时系数下逐线迭代求解+对所得的解施以欠松弛，再用欠松弛后的解去计算新的系数，常数，以进入下一层次的迭代。 实施：常把欠松弛处理纳入迭代过程，而不是在一个层次迭代完成后再行欠松弛。 )( ) ()()1(n p p n n n p n p t a b bt a t t -∑+=+ω )()1() 1()( n p p n n n p p t a b b t b a t a ω ωω -+++∑=+ ∑+=+')1('b b bt a t a n n n p p )('))(1(',n p p p p t a b b a a ωωω-+==，用交替方向线迭代法求解这一方程，就实现了SLUR 的迭代求解。为一般化起见，上式中b t n 上没有标以迭代层次的符号（J ，GS 时不相同）。 2、采用拟非稳态法 前面已指出，稳态问题的迭代解法与非稳态问题的步进法十分相似。对于非线性稳态问题，从代数方程的一组临时系数进入到另一组临时系数亦好象非稳态问题前进了一个时间层，非稳态问题的物理特性：系数热惯性越大（↑??=τρ/v c a o p ），温度变化越慢，仿此，对稳态非线性问题，可在离散方程中加入拟非稳态项，以减小未知量托两个层次间的变化，即 由 )()1()1()()(n p o p n n n p o p p n n n n p p n t a b b bt a t a V S b a b b bt a t V S b a ++∑=+?-∑?+∑=?-∑++ o p p n n p o p n n n p a V S b a t a b b bt a t +?-∑++∑= +) ()1( 一直进行到b t t n p ,收敛，虚拟时间步τ?的大小通过计算实践确定。 3、采用Jacobi 点迭代法 中止迭代的判据（该层次迭代）除前述变化率判据外，还可以规定迭代的轮数，例如规定进行4-6次ADI 线迭代就结束该层次上的计算。此时，用收敛速度低的丁迭代也就起到了欠松弛的作用。 五、迭代法的收敛速度 1、收敛速度 对给定的代数方程组（包括是临时系数的情形），采用不同的迭代方法求解时，使一定的初始误差缩小成α倍所需要的迭代轮数K 是不相的。1<α 《传热学》上机大作业 二维导热物体温度场的数值模拟 学校：西安交通大学 姓名：张晓璐 学号:10031133 班级：能动A06 一．问题（4-23） 有一个用砖砌成的长方形截面的冷空气通道，形状和截面尺寸如下图所示，假设在垂直纸面方向冷空气和砖墙的温度变化很小，差别可以近似的予以忽略。在下列两种情况下计算：砖墙横截面上的温度分布；垂直于纸面方向上的每米长度上通过墙砖上的导热量。 第一种情况：内外壁分别维持在10C ?和30C ? 第二种情况：内外壁与流体发生对流传热，且有C t f ?=101， )/(2021k m W h ?=，C t f ?=302，)/(422k m W h ?=，K m W ?=/53.0λ 二．问题分析 1.控制方程 02222=??+??y t x t 2.边界条件 所研究物体关于横轴和纵轴对称，所以只研究四分之一即可，如下图： 对上图所示各边界： 边界1：由对称性可知：此边界绝热，0=w q 。 边界2：情况一：第一类边界条件 C t w ?=10 情况二：第三类边界条件 )()( 11f w w w t t h n t q -=??-=λ 边界3：情况一：第一类边界条件 C t w ?=30 情况二：第三类边界条件 )()( 22f w w w t t h n t q -=??-=λ 三：区域离散化及公式推导 如下图所示，用一系列和坐标抽平行的相互间隔cm 10的网格线将所示区域离散化，每个交点可以看做节点，该节点的温度近似看做节点所在区域的平均温度。利用热平衡法列出各个节点温度的代数方程。 第一种情况： 内部角点： 数值计算大作业 题目一、非线性方程求根 1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。 （1）如果令()N t 表示在t 时刻的人口数目，β 表示固定的人口出生率，则人口数目满足微分方程() ()dN t N t dt β=，此方程的解为0()=t N t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成() ()dN t N t v dt β=+， 此方程的解为 0()=+ (1) t t v N t N e e βββ -； 假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到 410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。 435000 1564000=1000000(1) e e βββ + - 2.数学原理 采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果) (x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。 设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x 进行泰勒展开，有 . ))(()()(???+-'+≈k k k x x x f x f x f 于是方程0)(=x f 可近似地表示为 ))(()(=-'+k k x x x f x f 这是个线性方程，记其根为1k x +，则1k x +的计算公式为 传热学作业数值计算 数值计算matlab程序内容：>> tw1=10; % 赋初值tw2=20; c=1.5; p2=20; p1=c*p2; L2=40; L1=c*L2; deltaX=L2/p2; a=p2+1; b=p1+1; ti=ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn=ti; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n=ti(i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw1+tw2*sin(pi*(j-1)/(b-1)); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end ti=tn; end k ti max1 t2=ones(a,b); %求解析温度场 for i=a:-1:1 for j=1:1:b y=deltaX*(a-i); x=deltaX*(j-1); t2(i,j)=tw1+tw2*sin(pi*x/L1)*(sinh(pi*y/L1))/(sinh(pi*L2/L1)); end end t2 迭代次数k =706 数值解温度场ti 第一章: 1-1 对于附图所示的两种水平夹层，试分析冷、热表面 间热量交换的方式有何不同？如果要通过实验来测定夹层中流体的导热系数，应采用哪一种布置？ 解：（a ）中热量交换的方式主要有热传导和热辐射。 （b ）热量交换的方式主要有热传导，自然对流和热辐射。 所以如果要通过实验来测定夹层中流体的导热系数，应采用（ a ）布置。 1-2 一炉子的炉墙厚13cm ，总面积为20m 2 ，平均导热系数为 1.04w/m 〃k ，内外壁温分别是520 ℃及50 ℃。试计算通过炉墙的热损失。如果所燃用的煤的发热量是 2.09 ×10 4 kJ/kg ，问每天因热损失要用掉多少千克煤？ 解：根据傅利叶公式 每天用煤 1-3 在一次测定空气横向流过单根圆管的对流换热实验中，得到下列数据：管壁平均温度t w = 69 ℃，空气温度t f = 20 ℃，管子外径d= 14mm ，加热段长80mm ，输入加热段的功率8.5w ，如果全部热量通过对流换热传给空气，试问此时的对流换热表面传热系数多大？ 解：根据牛顿冷却公式 1-4宇宙空间可近似的看作0K 的真空空间。一航天器在太空中飞行，其外表面平均温度为250K ，表面发射率为0.7 ，试计算航天器单位表面上的换热量？ 解：航天器单位表面上的换热量 1-5附图所示的空腔由两个平行黑体表面组成，孔腔内抽成真空，且空腔的厚度远小于其高度与宽度。其余已知条件如图。表面 2 是厚δ= 0.1m 的平板的一侧面，其另一侧表面 3 被高温流体加热，平板的平均导热系数λ=17.5w/m ? K ，试问在稳态工况下表面3 的t w3 温度为多少？ 解： 表面1 到表面2 的辐射换热量= 表面2 到表面3 的导热量 第二章: 二维稳态导热的数值计算 2.1物理问题 一矩形区域，其边长L=W=1，假设区域内无内热源，导热系数为常数，三个边温度为T1=0，一个边温度为T2=1，求该矩形区域内的温度分布。 2.2 数学描述 对上述问题的微分方程及其边界条件为：2222T T 0x y ??+=?? x=0，T=T 1=0 x=1，T=T 1=0 y=0，T=T 1=0 y=1，T=T 2=1 该问题的解析解：112121(1)sin n n n sh y T T n L x n T T n L sh W L ππππ∞=??? ?---??? ?=? ?-????? ??? ∑ 2.3数值离散 2.3.1区域离散 区域离散x 方向总节点数为N ，y 方向总节点数为M ，区域内任一节点用I,j 表示。 2.3.2方程的离散 对于图中所有的内部节点方程可写为：2222,,0i j i j t t x y ??????+= ? ??????? 用I,j 节点的二阶中心差分代替上式中的二阶导数，得： +1,,-1,,+1,,-1222+2+0i j i j i j i j i j i j T T T T T T x y --+= 上式整理成迭代形式：()()22 ,1,-1,,1,-12222+2() 2()i j i j i j i j i j y x T T T T T x y x y ++=++++ (i=2,3……,N-1)，(j=2,3……,M-1) 补充四个边界上的第一类边界条件得：1,1j T T = (j=1,2,3……,M) ,1N j T T = (j=1,2,3……,M) ,1i j T T = (i=1,2,3……,N) 主讲陶文铨 西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER 2009年9月7日,西安 数值传热学 第一章绪论 课程简介 1. 教材－《数值传热学》第二版，2001 2. 学时－45学时理论教学；10学时程序教学 3. 考核－平时作业/计算机大作业： 考试－40/60；考查－60/40 4. 方法－开放，参与，应用 5. 助手－郭东之，周文静，李兆辉 有关的主要国外期刊 1.Numerical Heat Transfer, Part A-Applications; Part B- Fundamentals 2.International Journal of Numerical Methods in Fluids. https://www.wendangku.net/doc/101207505.html,puter & Fluids 4.Journal of Computational Physics 5.International Journal of Numerical Methods in Engineering 6.International Journal of Numerical Methods in Heat and Fluid Flow https://www.wendangku.net/doc/101207505.html,puter Methods of Applied Mechanics and Engineering 8.Engineering Computations 9.Progress in Computational Fluid Dynamics 10. Computer Modeling in Engineering & Sciences (CMES) 11.ASME Journal of Heat Transfer 12.International Journal of Heat and Mass Transfer 13.ASME Journal of Fluids Engineering 14.International Journal of Heat and Fluid Flow 15.AIAA Journal 东南大学能源与环境学院 课程作业报告 作业名称：传热学大作业——利用matlab程序解决热传导问题 院系：能源与环境学院 专业：建筑环境与设备工程 学号： 姓名： 2014年11月9日 一、题目及要求 1.原始题目及要求 2.各节点的离散化的代数方程 3.源程序 4.不同初值时的收敛快慢 5.上下边界的热流量（λ=1W/(m℃）） 6.计算结果的等温线图 7.计算小结 题目：已知条件如下图所示： 二、各节点的离散化的代数方程 各温度节点的代数方程 ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4 ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4 tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 三、源程序 【G-S迭代程序】 【方法一】 函数文件为： function [y,n]=gauseidel(A,b,x0,eps) D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; y=G*x0+f; n=1; while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end 命令文件为： A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0; -1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 水—水蒸汽两相相变界面的数值模拟 ——两相流动与热物理大作业 姓名张蛟龙_______ 学号201328013524021__ 班级物理308_____ 指导教师刘捷__ 完成时间_2014.5.8_ 水—水蒸汽两相相变界面的数值模拟报告 一．文献综述 作为化石资源的替代产品，核能的高效，清洁一直备受青睐，然而光环之下，核废料的处理不禁让人黯然神伤。强致命性辐射，动辄千年的半衰期，惯用的办法只能是深埋，等待下一代的聪明才智。与此同时，核废料的利用和加速衰减一直是核能大国们的研究重点。欧洲的ADS系统第六代散裂靶模型计划的目标就是要验证高水平的核废料转换的可行性。散裂靶作为连接加速器和核废料的装置需要工作在高辐射和高热流密度的条件下，因此散裂靶的设计是ADS系统研制最有挑战的部分。由加速器产生的高能质子流轰击靶核产生中子作为外源中子驱动和维持次临界堆的运行。散裂靶在极小的空间内需承受极大的热负荷，质子束通道与靶核的自由面相邻更加剧了设计难度。受材料限制，流体的温度不能超过550度，因此必须保证流体维持在一定的流量。但同时又要考虑高流速带来的飞溅和回流造成的局部温度过高。这一装置在水作为散裂靶的实验中获得了成功。二．问题描述 2.1.模型及尺寸 图1、欧洲液态金属散裂靶V0.10示意图[1] 如图1所示的欧洲加速器驱动次临界堆（ADS ）之无窗散裂靶示意图，液态铅铋合金从上方管间流下并汇合，形成两相界面，质子束由中间的真空管进入打在自由面上。此次模拟用的是水，详细物理背景见文献[1]。 2.2. 控制方程 连续性方程 动量方程 能量方程 三． Openfoam 求解 有关Openfoam 的下载和安装在老师给的安装指导的推荐网站上有详细的操作，在此就不赘述。网址为：https://www.wendangku.net/doc/101207505.html,/download/ubuntu.php 。 3.1. OpenFoam 求解简述 Openfoam 是一款基于linex 的开源可编程软件，其求解过程的关键是三个文件夹的设置，即0，constant 和system 。0文件夹里存放的是初始条件和边界条件设置文件；constant 文件夹里存放的是网格文件，物性参数和求解器模型；system 文件夹里存放的是求解过程控制，差分格式和代数方程求解器设置文件。以下就三个文件的设置展开简述初始条件、边界条件、物性参数，网格个数、疏密设置差分格式、界面捕获算法、气蚀模型等的选择和设置。 3.2. 0文件夹 包含有5个文件，分别为alph-water ，p_rgh ，U ，epsilon ，k ，详细设置见附录1，这里只着重强调在大作业完成过程中几个曾经连续考虑的点。 首先是参数的量纲设置。在Openfoam 文件中常会见到这样一行代码：dimensions [0 0 0 0 0]，这便是量纲，单位顺序依次是 [质量，长度，时间，温度，物质的量，电流，光强]。 其次是边界条件和初始条件的设置。在alph-water 中，alpha 代表水所占比例，参照userguide ，1时表示全部为液相，0时表示全部为气相。初始内部场的设置均为1，即起始时刻，散裂靶内部充满水。水入口是边界类型为“定值”，即 0)(=?? +??i i u x t ρρi b j ij i j i j i F x x p x u u t u +??+??-=??+??τρρετδρρ+=+-++??++??j j b j c ij k i i ij i i j j i i u F Q u p u u u e u x u u e t ] )2 1 ([)]21([ 热物理过程的数值模拟Numerical Simulation of Thermophysics Process 讲稿 主讲：李隆键 第一章概论 1．1流动与传热过程的予测方法及特点 流动、传热、燃烧问题是热工类各专业和机械类动力机械专业所研究和解决的主要问题之一，燃烧问题实际上是有化学反应的流动与传热问题，推而广之，在所有热物理过程中，几乎都涉及到流动、传热问题。 预测的重要性： ①在规定设计参数的相应的结构下，热物理过程是否满足要求，达到预定的指 标？要预测； ②优化设计，不同方案的比较，要预测； ③减少设计、生产、再设计和再生产的费用； ④减少设计更改； ⑤减少试验和测量次数。 问题的核心：速度场、温度场（传热量）、浓度场等。 一、热物理问题的予测方法：理论分析法、实验测定、数值模拟 1、理论分析 以数学分析为基础，求解描述热物理过程的定解问题，获得函数形式的解，表示求解区域内物理量连续分布的场（速度场、温度场、浓度场……）。 控制方程+单值条件（数学模型）→理论解（分析解，解析解） 根据解的准确程度，又可再分为： （1）精确分析解（严格解） 特点：函数形式的解；它在求解区域精确地满足定解问题。 具体解法：直接积分法、分离变量法、积分变换法、热源法、映射法。 （2）近似分析解法 特点：函数形式的解，在求解区域上近似地满足定解问题（但在总量上满足相应的守恒原理，动量守恒、动量守恒、能量守恒、质量守恒）。 具体解法：积分法（从积分方程出发） 变分近似解法 摄动法（从微分方程出发） 2、实验测定 （1）纯实验法 （2）相似理论实验法：同类相似，减少变量数目→减少工作量，得到规律性结计算传热学中国石油大学(华东)第四章大作业
传热学第四版课后题答案第四章
第二章 传热习题答案
传热学大作业报告 二维稳态导热
西安交通大学传热学大作业二维温度场热电比拟实验1
计算传热学数值模拟
热物理过程的数值模拟-计算传热学3.(DOC)
西安交通大学传热学大作业
数值计算大作业
哈尔滨工程大学传热学大作业数值计算matlab程序内容
传热学_杨茉_部分习题与解答
传热学上机C程序源答案之二维非稳态导热的数值计算
数值传热学chapter_1
传热学MATLAB温度分布大课后复习完整编辑版
两相流大作业
热物理过程的数值模拟-计算传热学1