数值传热学作业-第四章
《传热学》课后习题答案-第四章
t k i，j 1 t k i，j t k i，j 1 t k i , j r r rj rj r 2 r 2 rj r
并简化，可以得出与上式完全一样相同的结果。
4-7、 一金属短圆柱在炉内受热厚被竖直地移植到空气中冷却， 底面可以认为是绝热的。为用数值法确定冷却过程中柱体温 度的变化， 取中心角为 1rad 的区域来研究 （如本题附图所示） 。 已知柱体表面发射率，自然对流表面传热系数，环境温度， 金属的热扩散率，试列出图中节点（1，1） ， （M,1）(M,n)及 （M,N） 的离散方程式。 在 r 及 z 方向上网格是各自均分的。 解：应用热平衡法来建立四个节点点离散方程。 节点（1，1） ：
， 离散方程的建立 4-5、试将直角坐标中的常物性无内热源的二维稳态导热微分方程化为显式差分格式，并指 出其稳定性条件（ x y) 。 解：常物性无内热源二维非稳态方程微分方程为
4.3636t 2 2.53t1 1.8336t f
t2
2.53t f 1.8336t f
2t 2t t a x 2 y 2
Bi=0.1,1,10 的三种情况计算下列特征方程的根
n (n 1,2,6) :
n a Fo 2 0.2 并用计算机查明，当 时用式（3-19）表示的级数的第一项代替整个级数（计
算中用前六项之和来替代）可能引起 的误差。 解： n Bi 0.1 1.0 10
tan n
第四章
复习题 1、 试简要说明对导热问题进行有限差分数值计算的基本思想与步骤。 2、 试说明用热平衡法建立节点温度离散方程的基本思想。 3、 推导导热微分方程的步骤和过程与用热平衡法建立节点温度离散方程的过程十分相似， 为什么前者得到的是精确描述，而后者解出的确实近似解。 4、 第三类边界条件边界节点的离散那方程，也可用将第三类边界条件表达式中的一阶导数 用差分公式表示来建立。试比较这样建立起来的离散方程与用热平衡建立起来的离散方 程的异同与优劣。 5．对绝热边界条件的数值处理本章采用了哪些方法？试分析比较之． 6．什么是非稳态导热问题的显示格式？什么是显示格式计算中的稳定性问题？ 7．用高斯－塞德尔迭代法求解代数方程时是否一定可以得到收敛德解？不能得出收敛的解 时是否因为初场的假设不合适而造成？
数值传热学chapter_4a
主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER 2009年9月28日, 西安数值传热学第四章扩散方程的数值解及其应用(1)4.1 一维导热问题4.1.1一维稳态导热的通用控制方程4.1.3界面导热系数的确定方法4.1.4 一维非稳态导热控制方程的离散化4.1.2通用控制方程控制容积积分法的离散4.1.5 数学上的稳定未必导致物理上有意义的解一维稳态导热问题不同坐标系通用控制方程0 P P()0P x x Δ=i数学上的稳定未必导致物理上有意义的解无内热源一维非稳态导热，初场均匀，两表面0]T +代入下式：P（全隐格式）才能满足。
结论：数学上的稳定未必导致物理上有意义的解；推=xΔP Pa T+极坐标均可以表示成为：2.解决通用化的一种方案为写出适合于三种坐标系中系数的通用表达式，特引进两个辅助变量：（1）x –方向标尺因子，scaling factor ，x-方向的距离表示成为sx x δi 。
对直角、圆柱坐标规定1;sx ≡（2）y-方向引入一个名义半径，R 。
对直角坐标R ＝1，据此，东西导热距离为：sx xδi 东西导热面积为：R /y sxΔ对极坐标取;sx r =对圆柱与极坐标R ＝r三种二维正交坐标系中离散方程的统一表达式按这种方式编制程序时，只要设置一个变量MODE，4.3 源项与边界条件的处理4.3.1非常数源项的线性化处理1. 线性化方法4.3.2第二、三类边界条件的处理2. 线性化方法讨论3. 线性化方法应用实例1. 补充以边界节点代数方程的方法2. 附加源项法S= P2. 线性化方法讨论（1）对与被求解变量有关的非常数源项，线性化比假定为常数更合理：用*()PS f T =来表示P 的源项比落后一个迭代步；P C P T S S S =+（2）任何复杂的函数总可以用线性函数来近似逼近；线性又是建立线性代数方程所必须的；（3）是为保证代数方程迭代求解收敛所必须；0P S ≤P P nb nb a a b φφ=+∑P nb a a ≥∑P nb P a a S V =−Δ∑代数方程迭代求解收敛的充分条件是，因为可以确保代数方程迭代求解收敛。
数值传热学第四章编程题
4-5迭代法求解节点温度。
说明：此处给出的是C++程序代码，使用牛顿迭代法，迭代收敛精度1.0e-6；程序运行结果附后。
/*NHT 4-5 newton*created on 2012-10-19 by Sanye*/#include <cmath>#include <stdio.h>#include <iostream>using namespace std;int main(){ double funT=1.0,dfunT=1.0,temp1=1.0,temp2=1.0;double T=20.0;//primary valueint i=0; //for TEST!cout<<"primary T= "<<T<<" ;"<<endl;while(fabs(funT/dfunT)>=1.0e-6){ i++;if(i==1&&(T<=20.0))T=100.0;//in case unreasonable T;temp1=pow(T-20.0,0.25);temp2=pow(T-20.0,-0.75);funT=0.5*T-80+2*T*temp1-40*temp1;dfunT=0.5+2*temp1+0.5*T*temp2-10*temp2;T=T-funT/dfunT;cout<<" step "<<i<<";\tT3'-T3= "<<funT/dfunT<<" ;"<<endl;}cout<<"total steps: "<<i<<endl;cout<<"\tT[1]= "<<100<<"\tT[2]= "<<0.5*(T+130)<<"\tT[3]= "<<T<<endl;return 0;}运行结果：4-12 编写TDMA算法程序验证其正确性。
传热学第五第1-4章习题解答
《传热学》（第五版）第0章-第3 章习题解答第0章 绪论0-4、解答题略。
0-6 答:对流换热和对流不是同一现象.热对流是指:若流体有宏观运动，且内部存在温差，则由于流体各部分之间发生相对位移，冷热流体相互掺混而产生的热量传递现象，简称对流.如热空气往上升时,把热量传给上部空间的冷空气的流动属于对流.对流换热是指流体在与它温度不同的壁面上流动时，二者之间（流体与壁面之间）产生的热量交换现象。
它是导热与热对流同时存在的复杂热传递过程。
如暖气片周围的空气受热后,沿着教室墙壁的流动;热水在热力管道内的流动等属于对流换热.0-6 答:首先,冬季和夏季的最大区别在于室外温度的不同.夏季室外温度比室内温度高,因此通过墙壁的热量传递是由室外传向室内.而冬季室外气温比较比室内低, 通过墙壁的热量传递是由室内传向室外.因此冬季和夏季墙壁内表面的温度不同,夏季高而冬季低.因此人体在冬季通过辐射而与墙壁的散热比夏季高很多.人对冷暖感觉的衡量指标是散热量的大小而不是温度的高低,即当人体散热量少时感到热, 人体散热量多时感到冷.拉上窗帘后顿觉暖和,是因为窗帘起到了保温层的作用,减少了通过窗户向外散失的热量,故顿觉暖和!0-9 答:真空玻璃夹层:阻止热传导和对流换热;夹层内镀银:反射辐射热;热量如何通过瓶胆传到外界: 略瓶胆的玻璃尖嘴打破变得很差,因为空气进入夹层后,会由于空气与瓶胆壁面之间的对流换热而引起热量散失. 0-13：解： 61.0124161.036.08711121=++=++=h h R k λδ(m 2·K)/W 64.1610.011===k R k W/(m 2·K) 92.45)1018(64.1)(21=+=-=f f t t k q W/m 2 ∵)(111w f t t h q -= ∴47.178792.4518111=-=-=h q t t f w ℃ 又∵)(222f w t t h q -= ∴63.912492.4510222-=+-=+=h q t t f w ℃38.292.45⨯⨯==ΦqA =385.73 W0-14：解：4104.723452.0-⨯=⨯⨯==A R A λδ K/W (面积为A 2的平板表面上的热阻) 3104.4452.0-⨯===λδR (m 2·K)/W （单位面积热阻）431007.3104.4150285⨯=⨯-=∆=-R t q W/m 2 541084.161007.3⨯≈⨯⨯==ΦqA W0-15：解： ∵)(f w t t h q -= ∴15573511085=+=+=h q t t f w ℃ W7.20065.214.31050511023=⨯⨯⨯⨯=⋅⋅==Φ-lR q qA π0-17：解： （1）012.0851500011121=+=+=h h R (m 2·K)/W 3.83012.011===R k W/(m 2·K) 90963624)45500(3.83=⨯-⨯=∆=ΦtA k W（2）92820024)45500(85'=⨯-⨯=∆=ΦtA k W误差%2909636909636928200%100'≈-=⨯ΦΦ-Φ=ε （3）可以忽略，因为厚度很小，金属的导热系数较大，则导热热阻λδ很小。
第四章-传热习题参考答案
第四章 传热3 直径为φ60×3mm 的钢管用30mm 厚的软木包扎，其外又用100mm 厚的保温灰包扎，以作为绝热层。
现测得钢管外壁面温度为-110℃，绝热层外表面温度为10℃。
软木和保温灰的导热系数分别为0.043和0.07W/(m. ℃)，试求每米管长的冷损失量。
解：m W r r r r t L Q /25207.0)60/160ln(2043.0)30/60ln(101102)/ln(2)/ln(223112-=⨯⨯+⨯⨯--=+∆=πππλπλ4 蒸汽管道外包扎有两层导热系数不同而厚度相同的绝热层，设外层的平均直径为内层的两倍。
其导热系数也为内层的两倍。
若将两层材料互换位置，而假定其它条件不变，试问每米管长的热损失将改变多少？说明在本题情况下，哪一种材料包扎在内层较为合适？解：内层管内径r 1,外径r 2,外层管外径r 32(221r r +)=(223r r +) , 2312r r r r -=- 23125,3r r r r ==⇒ 122λλ=πλπλπλπλ4)3/5ln(2)3ln(2)/ln(2)/ln(11223112+∆=+∆=t r r r r t L Qπλπλπλπλ2)3/5ln(4)3ln(2)/ln(2)/ln('11123212+∆=+∆=tr r r r t L Q 25.1)3/5ln(3ln 2)3/5ln(23ln '=++=⇒Q Q 所以导热系数小的应该包扎在内层。
7 在并流换热器中，水的进出口温度分别为15℃和40℃，油的进、出口温度分别为150℃和100℃。
现因生产任务要求油的出口温度降至80℃，假设油和水的流量、进口温度和物性均不变，若原换热器的管长为1m 。
试求此换热器的管长增至若干米才能满足要求。
设换热器的热损失可忽略。
解：''''''''m m m m m m t tQ Q S S t S t S t KS t KS Q Q ∆∆⋅=⇒∆∆=∆∆= （1） 其中：4.110015080150)()(''21'2'1=--=--=T T C W T T C W Q Q ph h ph h （2） 又由：C t t t t t T T T T t t C W T T C W t t C W T T C W pc c ph h pc c ph h ︒=⇒--=--⇒⎭⎬⎫-=--=-50''')'()'()()(21212212112211221 C t m ︒=-=∆∴5.9260/135ln 60135C t m ︒=-=∆8.6930/135ln 30135' （3）将（2）（3）代入（1）即得。
数值传热学--作业
数值传热学大作业—淬火过程的瞬态热分析专业：材料工程班级：研1303班学号：S2*******指导教师：孙斌煜姓名：李康一、问题描述某零件材料为45钢，按照国标GB/T6912-1999规定的45钢推荐热处理制度为840C 。
淬火.600C 。
回火，淬火介质为水，试计算零件温度随时间的我变化曲线和最后时刻的温度场云图 （1）45钢弹性模量:200GPa 泊松比：0.3质量密度：78503/m kg膨胀系数：15.5e-6m/C 。
比热：448C J kg / 导热系数：70()C m W .*/ （2）水 密度：9963/m kg 比热：4185C J kg / 导热系数：2()C m W .*/水沸腾对流换热系数：1200()C m W .*2/初始45钢温度840，水的初始温度为20C 。
，水槽宽1m,中间位零件最大截面60mm ，下图为淬火过程的零件截面。
图-1二、创建模型1.建立分析项目（1）在Windows系统下执行“开始”—“所有程序”—“ANSYS14.0”—“Mechanical APDL（ANSYS）14.0”命令，启动Mechanical APDL（ANSYS）14.0，进入主界面。
（2）选择热分析过滤菜单GUI：选择菜单Main Menu —Preprocessor，弹出分析项目对话框，选择Thermal 热分析，如图2 所示，完成后单击OK按钮结束。
2.更改分析名称和标题（1）改变工作项目标题GUI:File→Change Title,弹出对话框，输入“Thermal01”如下图，单击OK结束。
（2）更改项目名称GUI:File→Change Title,弹出对话框，输入“Thermal01”下方的复选框，如下图所示，完成单击OK完成。
3.创建材料模型要点：创建模型顺序依次为工件，水（1）添加导热系数GUI:Main Menu →Preprocessor→Material Prop→Material Models→Thermal→Conductivity→Isotropic,弹出对话框,输入导热系数70，如下图完成后单击OK 结束输入（2）添加比热容GUI:Main Menu →Preprocessor→Material Prop→Material Models→Thermal→Specific Heat.弹出比热容输入对话框，在文本框中输入工件比热容448，如下图，成后单击OK结束输入（3）添加密度GUI:Main Menu →Preprocessor→Material Prop→Material Models→Thermal→Density，弹出密度输入对话框，在文本框中输入工件比热容7785，如下图，成后单击OK结束输入(4)创建材料2依照上述步骤添加水的比热容4185，密度996，导热系数25.选择单元GUI:Main Menu →Preprocessor→Element type→add/Edit/Delete→add,弹出下图所示单元对话框，选择Thermal Solide的Quad 8node77单元，按OK键结束设置单元选项GUI:Main Menu →Preprocessor→Element type→add/Edit/Delete→add，弹出Element type对话框，单击对话框中的Option，弹出设置单元对话框，在单元形状K3文本框选择Plane Thickness,如下图，单击OK结束关闭对话框。
数值传热学第4章作业
习题4－2图4-22 习题4-2插图[解]一维稳态导热问题的控制方程为：0=+⎪⎭⎫ ⎝⎛S dx dT dx d λ 4-2-1 该问题的边界条件为：()⎪⎩⎪⎨⎧=-=-==2,0,1001x T T h dx dT x T f λ 4-2-2分别对节点2，3进行离散，将已知数据代入离散格式中，得到方程组：130232=-T T 4-2-375432=+-T T 4-2-4 联立式（4-2-3）、式（4-2-4），可以解出2T ，3T ： 852=T ，403=T 。
下面验证总体守恒性：4-2-5右端3放出的热量为：()()30020401533=-⨯=-=f T T h Q 4-2-6在总体容积内部产生的热量为：2.0150 2.0300S Q S x =⨯∆=⨯=还需要证明左端是绝热条件： 节点2的热平衡为：21851000.550.5150757501T T xS x λ--+∆=+⨯=-+=∆ 左端绝热，所以计算结果符合总体能量守恒。
习题 4-5[解] 根据习题4-2的分析，可以得到节点2的离散方程：130232+=T T 4-5-1对于节点3，应用边界条件：()()1324330.510f f T T S T T T T xλδ--+⨯=-- 4-5-2式（4-5-2）可以整理成：()5432355751020T T T =+-- 4-5-3采用局部线性化方法，可以得到：()()()()515***444333331020102012.520T T T T T -=-+-- 4-5-4节点3的离散方程表示成：()()()51**44323335575 2.52012.52020T T T T T =++---- 4-5-5迭代求解得出：2382.82;35.64T T == 检验热平衡：内热源生成热1300φ=； 右端散热5/4210(35.6420)311.0h T φ=∆=-=左端散热382.821000.5150510.91φ-=⨯+⨯=-所以123()30031110.90φφφ-+=-+≅不作热平衡扣0.5 分。
数值传热学习题答案(汇总版)
2-4-9
= rP rS
式（2-4-9）也可以写成 a PTP = a E TE + aW TW + b 的形式。而且两种结果是一致的。
2—6：
n n TE −TW dT P , n = 解：将 ， dx 2x n n TE −2TPn + TW d 2T P , n = , dx2 x 2
dk = f (x ) 代入原方程，得： dx
令
2-4-4
rk rk a E = ， aW = ， a P = a E + aW ， b x w x e
= SrP r ，
式（2-4-4）可以写成 a PTP = a E TE + aW TW + b 的形式。 2. 再用 Taylor 展开法导出 k
2 2 uE + uP u = ， 2 2 e
2 2 uW + uP u = 2 2 w
t u ut N − uP y = (y ) ， n n
t
t ut u p − uS y = (y ) 。 s s
t
(y ) n = (y ) s = y
n n n n TE −TW TE −2TPn + TW k + f (x ) +S=0 整理得： 2x x 2
4kT P= 2k + xf ( x)T E+2k − xf ( x)T W +2x 2 S
− 2k 时， a E 会成为负值， x 2k 当 f(x）> 时， aW 会成为负值。 x
rk dr = rk r r dr dr dr
w
e
1 d
传热学课后习题第四章答案
选择步长 x y ，又边界为绝热 h 0 ，整理有
t i 1, j t i , j
即
t i , j 1 t i , j 2
t i , j 1 t i , j 2
0
2ti 1, j ti, j 1 ti, j 1 4ti, j 0
4-2 解：根据热平衡方程有
t 2 t10 t5 t 7 4t 6 0
对节点 7
t 3 t11 t 6 t8 4t 7 0
对节点 10
t 6 t14 t9 t11 4t10 0
对节点 11
t 7 t15 t10 t12 4t11 0
对节点 14
t10 t18 t13 t15 4t14 0
hxt f t i , j
t i , j 1 t i , j 2
hx
0
(t i 1, j t i , j 1 ) 2(1
h
x)t i , j 2
tf 0
4-3 解：根据已知条件，划分网格如图所示 第一类边界条件
t w 200C
第二类边界条件
3 200 （2t 4 t 8 100 ） 14 3 3 200 t8 （2t 5 t 7 t 9 ） 14 3 t7
节点 9 为第三类边界条件下的外拐角边界节点
ti1, j ti, j 1 (2Bi 2)ti, j 2Bi t f 0
对节点 20
t16 t 20 t t 24 20 (t19 t 20 ) hy t f t 20 0 2 2
对节点 16
数值传热学陶文铨主编第二版习题答案
数值传热学4-9章习题答案习题4－2一维稳态导热问题的控制方程：022=+∂∂S xTλ 依据本题给定条件，对节点2节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程： 节点1： 1001=T节点2： 1505105321-=+-T T T 节点3：75432=+-T T求解结果：852=T ，403=T对整个控制容积作能量平衡，有：02150)4020(15)(3=⨯+-⨯=∆+-=∆+x S T T h x S q f f B即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果25.03)(10f T T h -⨯=，则各节点离散方程如下：节点1： 1001=T节点2： 1505105321-=+-T T T节点3：25.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T对于节点3中的相关项作局部线性化处理，然后迭代计算； 求解结果：818.822=T ，635.353=T （迭代精度为10-4）迭代计算的Matlab 程序如下： x=30; x1=20;while abs(x1-x)>0.0001a=[1 0 0;5 -10 5;0 -1 1+2*(x-20)^(0.25)]; b=[100;-150; 15+40*(x-20)^(0.25)]; t=a^(-1)*b;x1=x;x=t(3,1);endtcal=t习题4-12的Matlab程序%代数方程形式A i T i=C i T i+1+B i T i-1+D imdim=10;%计算的节点数x=linspace(1,3,mdim);%生成A、C、B、T数据的基数；A=cos(x);%TDMA的主对角元素B=sin(x);%TDMA的下对角线元素C=cos(x)+exp(x); %TDMA的上对角线元素T=exp(x).*cos(x); %温度数据%由A、B、C构成TDMAcoematrix=eye(mdim,mdim);for n=1:mdimcoematrix(n,n)=A(1,n);if n>=2coematrix(n,n-1)=-1*B(1,n);endif n<mdimcoematrix(n,n+1)=-1*C(1,n);endend%计算D矢量D=(coematrix*T')';%由已知的A、B、C、D用TDMA方法求解T%消元P(1,1)=C(1,1)/A(1,1);Q(1,1)=D(1,1)/A(1,1);for n=2:mdimP(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1));Q(1,n)=(D(1,n)+B(1,n)*Q(1,n-1))/(A(1,n)-B(1,n)*P(1,n-1)); end%回迭Tcal(1,mdim)=Q(1,mdim);for n=(mdim-1):-1:1Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=[T;Tcal];%绘图比较给定T值和计算T值plot(Tcal,'r*')hold onplot(T)结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ 对于三种无量纲定义w b w T T T T --=Θ、∞∞--=ΘT T T T w 、w w T T T T --=Θ∞进行分析如下1）由wb wT T T T --=Θ得：w w b T T T T +Θ-=)(由T 可得：x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(])[(rT r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，除了w T 均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的； 2）由∞∞--=ΘT T T T w 得： ∞∞+Θ-=T T T T w )(由T 可得：xT x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(rT r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[( 由b T 与r 无关、Θ与x 无关以及x T ∂∂、rT∂∂的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，则该无量纲温度定义是可以用分离变量法的； 3）由wwT T T T --=Θ∞得： w w T T T T +Θ-=∞)(由T 可得：xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(T r T T r T T T r T w w w w ∂-+∂Θ∂-=∂+Θ-∂=∂∂∞∞1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：S r r r r r r x x w r v r r r u x +∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂+∂∂+∂∂)(1)(1)()(1)(1)(θφλθφλφλφρθφρφρ x 、r 和θ分别是圆柱坐标的3个坐标轴，u 、v 和w 分别是其对应的速度分量，其中x 是管内的流动方向；对于管内的层流充分发展有：0=v 、0=w ，0=∂∂xu； 并且x 方向的源项：x pS ∂∂-=r 方向的源项：r pS ∂∂-= θ方向的源项：θ∂∂-=pr S 1 由以上分析可得到圆柱坐标下的动量方程： x 方向： 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x pu r r r u r r r θλθλ r 方向：0=∂∂r pθ方向：0=∂∂θp边界条件： R r =，0=u0=r ，0=∂∂r u ；对称线上，0=∂∂θu 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂Tr r r T r r r x T uc p 边界条件： R r =，w q r T =∂∂λ；0=r ，0=∂∂rTπθ/0=，0=∂∂-θλT2）定义无量纲流速：dxdp R uU 2-=λ并定义无量纲半径：R r /=η；将无量纲流速和无量纲半径代入x 方向的动量方程得：0))1((1))1((122=∂∂-∂-∂∂∂+∂-∂∂∂xp U dx dp R R R R U dx dp R RR R θληλθηηλληηη 上式化简得：01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：1=η，0=U0=η，0=∂∂ηU ；对称线上，0=∂∂θU定义无量纲温度：λ/0R q T T b-=Θ其中，0q 是折算到管壁表面上的平均热流密度，即：Rq q wπ=0； 由无量纲温度定义可得：b T Rq T +Θ=λ将T 表达式和无量纲半径η代入能量方程得：)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p （1）由热平衡条件关系可以得：mm m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T uc 020221221)(===∂∂=∂∂ππρρρ 将上式代入式（1）可得：)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=m U U 边界条件：0=η，0=∂Θ∂η；1=η，R q q w πη10==∂Θ∂0=θ，0=∂Θ∂θ；πθ=，0=∂Θ∂θ单值条件： 由定义可知：0/0=-=ΘλR q T T b b b 且： ⎰⎰Θ=ΘAAb UdAUdA即得单值性条件：0=Θ⎰⎰AA UdAUdA 3）由阻力系数f 及Re 定义有：228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m me νρ 且：m W b m W b m W R q T T D T T q Nu ,0,,0~2)/(2Θ=-=-=λλ5－21．一维稳态无源项的对流－扩散方程如下所示：xx u 22∂∂Γ=∂∂φφρ （取常物性）边界条件如下：L L x x φφφφ====,;,00上述方程的精确解如下：11)/(00--=--⋅PeL x Pe L e e φφφφ Γ=/uL Pe ρ 2．将L 分成20等份，所以有：∆=P Pe 201 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下： 1） 中心差分中间节点： 2)5.01()5.01(11-∆+∆++-=i i i P P φφφ 20,2 =i2） 一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ 20,2 =i3） 混合格式当1=∆P 时，中间节点：2)5.01()5.01(11-∆+∆++-=i i i P P φφφ20,2 =i当10,5=∆P 时，中间节点： 1-=i i φφ 20,2 =i 4） QUICK 格式*12111)35(8122121⎥⎦⎤⎢⎣⎡---++++++=+--∆∆-∆∆+∆i i i i i i i P P P P P φφφφφφφ 2≠i *1111)336(8122121⎥⎦⎤⎢⎣⎡--++++++=+-∆∆-∆∆+∆i i i i i i P P P P P φφφφφφ 2=i数值计算结果与精确解的计算程序如下：%except for HS, any other scheme doesnt take Pe<0 into consideration %expression of exact solutiony=dsolve('a*b*Dy=c*D2y','y(0)=y0,y(L)=yL','x')y=subs(y,'L*a*b/c','t')y=simple(subs(y,'a*b/c*x','t*X'));ysim=simple(sym(strcat('(',char(y),'-y0)','/(yL-y0)')))y=sym(strcat('(',char(ysim),')*(yL-y0)','+y0'))% in the case of Pe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%grid Pe numbertt=[1 5 10];%dimensionless lengthm=20;%mdim is the number of inner nodemdim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculationy0=1;yL=2;%cal exact solutionfor n=1:size(tt,2)t=m*tt(1,n);if t==0yval1(n,:)=eval(y1);elseyval1(n,:)=eval(y);endend%extra treatment because max number in MATLAB is 10^308if max(isnan(yval1(:)))yval1=yval1';yval1=yval1(:);indexf=find(isnan(yval1));for n=1:size(indexf,1)if rem(indexf(n,1),size(X,2))==0yval1(indexf(n),1)=yL;elseyval1(indexf(n),1)=y0;endendyval1=reshape(yval1,size(X,2),size(yval1,1)/size(X,2));yval1=yval1';end%CD solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*tt(1,n))*y0;d(n,mdim)=0.5*(1-0.5*tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval2=TDMA(a,b,c,d,mdim);yval2=[repmat([1],size(tt,2),1),yval2,repmat([2],size(tt,2),1)]; Fig(1,X,yval1,yval2,tt);title('CD Vs. Exact Solution')% FUS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval3=TDMA(a,b,c,d,mdim);yval3=[repmat([1],size(tt,2),1),yval3,repmat([2],size(tt,2),1)]; Fig(2,X,yval1,yval3,tt);title('FUS Vs. Exact Solution')% HS solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);if t>2b(n,:)=repmat([0],1,mdim);c(n,:)=repmat([1],1,mdim);d(n,1)=y0;elseif t<-2b(n,:)=repmat([1],1,mdim);c(n,:)=repmat([0],1,mdim);d(n,mdim)=yL;elseb(n,:)=repmat([0.5*(1-0.5*t)],1,mdim);c(n,:)=repmat([0.5*(1+0.5*t)],1,mdim);d(n,1)=0.5*(1+0.5*t)*y0;d(n,mdim)=0.5*(1-0.5*t)*yL;endendc(:,1)=0;b(:,mdim)=0;% numerical cal by using TDMA subfuctionyval4=TDMA(a,b,c,d,mdim);yval4=[repmat([1],size(tt,2),1),yval4,repmat([2],size(tt,2),1)]; Fig(3,X,yval1,yval4,tt);title('HS Vs. Exact Solution')%QUICK Solutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);for n=1:size(tt,2)t=tt(1,n);b(n,:)=repmat([1/(2+t)],1,mdim);c(n,:)=repmat([(1+t)/(2+t)],1,mdim);d(n,1)=(1+tt(1,n))/(2+tt(1,n))*y0;d(n,mdim)=1/(2+tt(1,n))*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval5=zeros(size(tt,2),mdim);yval5com=yval5+1;counter=1;%iterativewhile max(max(abs(yval5-yval5com)))>10^-10if counter==1yval5com=TDMA(a,b,c,d,mdim);endfor nn=1:size(tt,2)for nnn=1:mdimif nnn==1d(nn,nnn)=((6*yval5com(nn,nnn)-3*y0-3*yval5com(nn,nnn+1))*tt(1,nn))/(8*(2+tt(1, nn)))+((1+tt(1,nn))/(2+tt(1,nn))*y0);elseif nnn==2d(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-y0)*tt (1,nn))/(8*(2+tt(1,nn)));elseif nnn==mdimd(nn,nnn)=((5*yval5com(nn,nnn)-3*yL-yval5com(nn,nnn-1)-yval5com(nn,nnn-2))*tt (1,nn))/(8*(2+tt(1,nn)))+(1/(2+tt(1,nn))*yL);elsed(nn,nnn)=((5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-yval5 com(nn,nnn-2))*tt(1,nn))/(8*(2+tt(1,nn)));endendendyval5=TDMA(a,b,c,d,mdim);temp=yval5;yval5=yval5com;yval5com=temp;counter=counter+1;endyval5=yval5com;yval5=[repmat([1],size(tt,2),1),yval5,repmat([2],size(tt,2),1)];Fig(4,X,yval1,yval5,tt);title('QUICK Vs. Exact Solution')%-------------TDMA SubFunction------------------function y=TDMA(a,b,c,d,mdim)%form a b c d resolve yval2 by using TDMA%eliminationp(:,1)=b(:,1)./a(:,1);q(:,1)=d(:,1)./a(:,1);for n=2:mdimp(:,n)=b(:,n)./(a(:,n)-c(:,n).*p(:,n-1));q(:,n)=(d(:,n)+c(:,n).*q(:,n-1))./(a(:,n)-c(:,n).*p(:,n-1));end%iterativey(:,mdim)=q(:,mdim);for n=(mdim-1):-1:1y(:,n)=p(:,n).*y(:,n+1)+q(:,n);end%-------------ResultCom SubFunction------------------ function y=ResultCom (a,b,c)for n=1:max(size(c,2))y(2*n-1,:)=a(n,:);y(2*n,:)=b(n,:);end%-------------Fig SubFunction------------------ function y=Fig(n,a,b,c,d)figure(n);plot(a,b);hold onplot(a,c,'*');str='''legend(';for n=1:size(d,2)if n==size(d,2)str=strcat(str,'''''Pe=',num2str(d(1,n)),''''')''');elsestr=strcat(str,'''''Pe=',num2str(d(1,n)),''''',');endendeval(eval(str));精确解与数值解的对比图，其中边界条件给定10=φ，2=L φ。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
4-1解：采用区域离散方法A 时；内点采用中心差分123278.87769.9T T T ===22d T T=0dx - 有 i+1i 122+T 0i i T T T x---=∆ 将2点,3点带入 321222+T 0T T T x --=∆ 即321209T T -+= 432322+T 0T T T x --=∆4321322+T 0T T T x --=∆ 即4321209T T T -+-= 边界点4（1）一阶截差 由x=1 1dT dx =，得 4313T T -=（2）二阶截差 11B M M q x x xT T S δδλλ-=++V所以 434111. 1.36311T T T =++即 43122293T T -=采用区域离散方法B22d TT=0dx - 由控制容积法 0w edT dT T x dT dT ⎛⎫⎛⎫--∆= ⎪ ⎪⎝⎭⎝⎭ 所以代入2点4点有322121011336T T T T T ----= 即 239028T T -=544431011363T T T T T ----= 即3459902828T T T -+= 对3点采用中心差分有432322+T 013T T T --=⎛⎫⎪⎝⎭即2349901919T T T -+= 对于点5 由x=11dT dx =，得 5416T T -= （1）精确解求左端点的热流密度由 ()21x x eT e e e -=-+所以有 ()2220.64806911x xx x dT e e q e e dxe e λ-====-+=-=++ （2）由A 的一阶截差公式210.247730.743113x T T dT q dxλ=-=-==⨯= （3）由B 的一阶截差公式0.216400.649213x dTq dxλ=-=-== （4）由区域离散方法B 中的一阶截差公式：210.108460.6504()B BT T dT dx x δ-⎛⎫==⨯= ⎪⎝⎭ 通过对上述计算结果进行比较可得：区域离散B 有控制容积平衡法建立的离散方程与区域离散方程A 中具有二阶精度的格式精确度相当！4-3解: 对平板最如下处理:1 2 3 4由左向右点分别表述为1、2、3、4点，x 的正方向为由左向右； 控制方程为λd 2tdx +S =0 （1）边界条件为X=0，T=75℃；X=0.1，λdTdx +ℎ(T −T f )=0；则2、3点采用二阶截差格式，有 则有以下两式：λT3−2T2+T1∆x+S=0（2）λT4−2T3+T2∆x2+S=0（3）一阶截差公式可由λdTdx+ℎ(T−T f)=0变形得到λ(T4−T3∆x)=h（T4−T f）再变形得到T4=[T3+h×∆xλT f]/（1+h×∆xλ）（4）二阶截差公式可以联立λT5−2T4+T3∆x2+S=0和λ(T5−T32∆x)=h（T4−T f），可得以下公式T4=[T3+∆x2S2λ+h×∆xλ]/（1+h×∆xλ）（5）分别联立2、3、4式与2、3、5式，把S=50×103W/m3，λ=10W/m∙℃，h=50 W/m∙℃，T f=25℃，T1=75℃，∆x= 1/30带入到式子中，则有联立2、3、4式的解为：T2=78.58℃，T3=76.59℃，T4=69.03℃联立3、4、5式的解为：T2=80.42℃，T3=80.28℃，T4=74.58℃对控制方程进行积分，并将边界条件带入，则有关于T的方程T=−2500x2+250x+75（6）把x2=130,x3=230,x3=0.1代入上述6式则有：T2=80.56℃，T3=80.56℃，T4=75.1℃相比之下，对右端点采用二阶截差的离散更接近真实值4-4解：对平板作如下分析：1 2 3 4 5 由左向右分别对点编号为1、2、3、4、5 控制方程与4-3相同，为λd 2tdx +S =0 （1）边界条件为X=0，T=75℃；X=0.1，λdTdx +ℎ(T −T f )=0；设1点和2点的距离为∆x ，另1点对2点进行泰勒展开，有d 2t dx =（T 1−T 2+dT dx ∆x ）2∆x其中dT dx=T 3−T 22∆x，则有λ2T 1−3T 2+T 3∆x 2+S =0 （2）对3点进行离散有λT 4−2T 3+T 2∆x 2+S =0 （3）对右端点有： [a p +A 1ℎ+(δx )5λ]T 4=a w T 3+[S/∆x +AT f 1ℎ+(δx )5λ]代入数据有T 3−3T 2+155.56=0 T 4−2T 3+T 2=−5.56342.85T4-300T3=1681解得：T2=78.1℃，T3=78.7℃，T4=73.8℃由导热定律有T4−T3∆x =2T5−T4∆x则有T5=71.35℃4—12编写程序：M=rand(10,3)A=M(:,1);B=M(:,2);C=M(:,3);B(10)=0;C(1)=0;T=12:21;D(1)=A(1)*T(1)-B(1)*T(2)for i=2:9;D(i)= A(i)*T(i)-B(i)*T(i+1)-C(i)*T(i-1)endD(10)= A(10)*T(10)-C(10)*T(9);P(1)=B(1)/A(1);Q(1)= D(1)/A(1);for i=2:10;P(i)=B(i)/(A(i)-C(i)*P(i-1));Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1)); endfor i=10:-1:2;t(10)=Q(10);t(i-1)=P(i-1)*t(i)+Q(i-1);enddisp(D(1:10))disp(T(1:10))disp(t(1:10))运行结果：由运行结果可知：无论系数怎样变化，T与t都是一致的。
编程题：具体编程如下：n=input('请输入网格数：n=');dt=input('请输入时间步长： dt =');% 设定区域的长宽L=input('请输入边界长度：L =');a=input('请输入材质热扩散率：a =');S=input('请输时间层数：S =');N=n+1 ;% 节点数dx=L/n;% 空间步长F=a*dt/(dx^2);% 傅里叶数% %%%%%%%%%%%%%%%%%%%%%%初始条件%%%%%%%%%%%%%%%%%%%%%%%%% for t=1:1:Sfor i=1:1:Nif t==1T(i,t)=100;else if i==1||i==NT(1,t)=0;T(N,t)=10;endendendend% %%%%%%%%%%%%%%%%%%%%%初始条件%%%%%%%%%%%%%%%%%%%%%%%%%A=2*F+1;B=F;C=F;for t=2:1:S% %%%%%%%%%%%%%%%%%%%计算各个时层的P、Q%%%%%%%%%%%%%%%%%%% for i=1:1:N-1D(i,2)=T(i,1);if i==1P(1,t)=0;Q(1,t)=T(1,t);elseP(i,t)=B/(A-C*P((i-1),t));Q(i,t)=(D(i,t)+C*Q((i-1),t))/(A-C*P((i-1),t));endend% %%%%%%%%%%%%%%%%%%%%%计算各个时层的P、Q%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%计算各个时层的温度%%%%%%%%%%%%%%%%%%%for i=N-1:-1:2T(i,t)=P(i,t)*T((i+1),t)+Q(i,t);D(i,(t+1))=T(i,t);endendx=0:0.01:1;% 生成其中几个时层图像plot(x,T(1:101,1),':bo',x,T(1:101,2),':ro',x,T(1:101,3),':bo',x,T(1:1 01,4),':ro',x,T(1:101,5),':bo',x,T(1:101,6),':ro')xlabel('x');ylabel('T');title('一维非稳态热传导方程的数值解')输入节点数为100,时间步长为0.1,边界长度为1,热扩散率为1,时间层数为50输入节点数为100,时间步长为1,边界长度为1,热扩散率为1,时间层数为50与显式比较输入节点数为100,时间步长为0.1例题：温度为f T =90 80 ℃的水进入一管径为0.1 0.12m 的长圆管作层流流动，平均速度为0.6m/s 。
管外受到温度为∞T =30 ℃的流体的冷却，e h =8w/(2m ℃)。
设换热已进入充分发展阶段：1、试确定管内对流换热系数e h ；2、给出x=1m 截面上水温沿半径方向的分布；3、平均水温沿x 方向的分布。
要求：确定网格独立解的节点数。
由题意可列出方程：)(1r T r r r x T u c p ∂∂∂∂=∂∂λρ边界条件：⎪⎩⎪⎨⎧-=∂∂-==∂∂=∞)(,r 0,0T T h r T R rT r e λ无量纲化：∞∞--=ΘT T T T b ，R r =η，PeR x X ⋅=,a Ru P m e 2=,)1(22η-=m u u ⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⋅-⋅=Λ-====-Λ+⎰=10212))1(4/(1)0,00)1()(ηηηφφηφηφηηφηηφηηηd Bi d d d d d d d d w对控制方程积分，写成追赶法的形式：r T r r T r r T r r i i i ∆-Λ+-∆+∆+=∆-+)1()21()21(22*11ηφη 边界条件：⎩⎨⎧=∆⋅+=-11121)1(L L i T T r B T T 查得90 80 ℃水的导热系数λ=0.68，导温系数71062.1-⨯=α,u m =0.6,R=0.05,设初始的Φ*=1，代入求解得到Φ，直到 3*10-≤-φφφ )(∞-Λ-=T T dX dT b bXb e C T Λ-⋅+=30 由边界条件X=0时，T b =90 C=60 51037.06030⨯Λ-⋅+=xb eT 取x=1时，解得T b 30)30(+-Λ=b T T φ编程为： %phi 为φ，eta 为η,phi1为φ*L1=input('输入节点个数：');R=0.05;%外节点法划分网格n=R/(L1-1);i=2:L1;X(i)=(i-1)*n;X(1)=0;X(L1)=R;eta=X/R;eta(1)=0;%迭代部分phi1=ones(1,L1);for i=2:L1-1A=2*X/n;B=0.5+X/n;C=X/n-0.5;Lambda=1/(4*sum(phi1.*eta.*(1-eta.^2)*(n/R))); S=Lambda.*eta.*(1-eta.^2).*phi1;D=S.*n;endA(1)=1;B(1)=1;C(1)=0;D(1)=0;P(1)=B(1)/A(1);Q(1)=D(1)/A(1);A(L1)=1+8*n/0.68;B(L1)=0;C(L1)=1;D(L1)=0;%追赶法解方程组for i=2:L1P(i)=B(i)/(A(i)-C(i)*P(i-1));Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1));endphi(L1)=Q(L1);for i=L1:-1:2phi(i-1)=P(i-1)*phi(i)+Q(i-1);endwhile abs((phi-phi1)./phi)>=10e-3phi1=phi;Lambda=1/(4*sum(phi1.*eta.*(1-eta.^2).*(n/R))); S=Lambda.*eta.*(1-eta.^2).*phi1;D=S.*n;for i=2:L1P(i)=B(i)/(A(i)-C(i)*P(i-1));Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1));endphi(L1)=Q(L1);for i=L1:-1:2phi(i-1)=P(i-1)*phi(i)+Q(i-1);endendTb=30+60*exp(-Lambda/3.7037e+5);T=phi*Lambda*(Tb-30)+30;h=-0.68*(T(L1-1)-T(L1))/n/(T(L1)-Tb) figure(1);plot(T,X)title('1m处水温延半径的分布');xlabel('温度');ylabel('半径');%平均水温x=linspace(0,10);Tb=30+60*exp(-Lambda.*x/3.7037e+5); figure(2);plot(x,Tb)title('平均水温延x的分布');xlabel('圆管长度');ylabel('平均水温');输入节点个数：100可得：h =29.8341x=1m截面上水温沿半径方向的分布图如下：平均水温沿x方向的分布。