数值传热学陶文铨第四章作业(完整资料).doc

1-9 一砖墙的表面积为122m ，厚为260mm ，平均导热系数为1.5W/（m.K ）。

设面向室内的表面温度为25℃，而外表面温度为-5℃，试确定次砖墙向外界散失的热量。

解：根据傅立叶定律有：WtA9.207626.05)(25125.1=--⨯⨯=∆=Φδλ1-12 在一次测定空气横向流过单根圆管的对流换热实验中，得到下列数据：管壁平均温度t w =69℃，空气温度t f =20℃，管子外径 d=14mm ，加热段长 80mm ，输入加热段的功率8.5w ，如果全部热量通过对流换热传给空气，试问此时的对流换热表面传热系数多大？ 解：根据牛顿冷却公式()f w t t rlh q -=π2所以()f w t t d qh -=π＝49.33W/(m 2.k)1-20 半径为0.5 m 的球状航天器在太空中飞行，其表面发射率为0.8。

航天器内电子元件的散热总共为175W 。

假设航天器没有从宇宙空间接受任何辐射能量，试估算其表面的平均温度。

解：电子原件的发热量＝航天器的辐射散热量即：4T Q εσ=4A QT εσ=∴ =187K 热阻分析 ;;2-4 一烘箱的炉门由两种保温材料A 及B 组成，且B A δδ2=(见附图)。

已知)./(1.0K m W A =λ,)./(06.0K m W B =λ,烘箱内空气温度4001=f t ℃，内壁面的总表面传热系数)./(501K m W h =。

为安全起见，希望烘箱炉门的 外表面温度不得高于50℃。

设可把炉门导热作为一维问题处理，试决定所需保温材料的厚度。

环境温度=2f t 25℃，外表面总传热系数)./(5.922K m W h =。

解：热损失为()()22111f f BBA A fwf t t h t t h t t q -+-=+-=λδλδ又50=fw t ℃；B A δδ=联立得m m B A 039.0;078.0==δδ2-9 双层玻璃窗系由两层厚为6mm 的玻璃及其间的空气隙所组成，空气隙厚度为8mm 。

习题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 分。

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 x T T S δδλλ-=++V 所以 434111. 1.36311T T T =++ 即 43122293T T -=采用区域离散方法B 22d T T=0dx - 由控制容积法 0w e dT dT T x dT dT ⎛⎫⎛⎫--∆= ⎪ ⎪⎝⎭⎝⎭ 所以代入2点4点有322121011336T T T T T ----= 即 239028T T -= 544431011363T T T T T ----= 即 34599 02828T T T -+=对3点采用中心差分有 432322+T 013T T T --=⎛⎫ ⎪⎝⎭ 即 2349901919T T T -+= 对于点5 由x=1 1dT dx =，得 5416T T -= （1）精确解求左端点的热流密度由 ()21x x e T e e e -=-+ 所以有 ()220020.64806911x x x x dTe e q e e dx e e λ-====-+=-=++ （2）由A 的一阶截差公式（3）由B 的一阶截差公式（4）由区域离散方法B 中的一阶截差公式：通过对上述计算结果进行比较可得：区域离散B 有控制容积平衡法建立的离散方程与区域离散方程A 中具有二阶精度的格式精确度相当！ 4-3解：将平板沿厚度方向3等分，如图由题可知该导热过程可看作无限大平板的一维稳态有源导热问题，则控制方程为x=0， T 0=75℃x=0.1 dT =h(T-T )dx f λ-1点 ，2点采用中心差分有21022+T 0T T S xλ-+=∆ （1） 32122+T 0T T S x λ-+=∆ （2） 右端点采用一阶截差的离散231f hx T T T x h λλ⎡⎤+⎢⎥⎣⎦=⎛⎫+ ⎪⎝⎭V （3） 右端点采用二阶截差的离散代入（1）（2）（3）得1223132280.62 5.67625T T T T T T T -=--=-= 解得123278.87769.9T T T ===代入（4）得解得 12380.6380.6675.1T T T ===精确解 22d T +S=0dxλ （4） x=0， T 0=75℃ （5） x=0.1 dT =h(T-T )dxf λ- （6）代入数据积分的将 x 1=10.13⨯，x 2=20.13⨯, x 3=0.1T 1=80.56 T 2=80.56 T 3=75.1通过比较可得右端点采用二阶截差的离散更接近真实值。

传热学思考题参考答案第一章：1、用铝制水壶烧开水时，尽管炉火很旺，但水壶仍安然无恙。

而一旦壶内的水烧干后水壶很快就被烧坏。

试从传热学的观点分析这一现象。

答：当壶内有水时，可以对壶底进行很好的冷却（水对壶底的对流换热系数大），壶底的热量被很快传走而不至于温度升得很高；当没有水时，和壶底发生对流换热的是气体，因为气体发生对流换热的表面换热系数小，壶底的热量不能很快被传走，故此壶底升温很快，容易被烧坏。

2、什么是串联热阻叠加原则，它在什么前提下成立以固体中的导热为例，试讨论有哪些情况可能使热量传递方向上不同截面的热流量不相等。

答：在一个串联的热量传递过程中，如果通过每个环节的热流量都相同，则各串联环节的总热阻等于各串联环节热阻的和。

例如：三块无限大平板叠加构成的平壁。

例如通过圆筒壁，对于各个传热环节的传热面积不相等，可能造成热量传递方向上不同截面的热流量不相等。

第二章：1、扩展表面中的导热问题可以按一维问题处理的条件是什么有人认为,只要扩展表面细长,就可按一维问题处理,你同意这种观点吗|答：条件：（1）材料的导热系数，表面传热系数以及沿肋高方向的横截面积均各自为常数（2）肋片温度在垂直纸面方向（即长度方向）不发生变化，因此可取一个截面（即单位长度）来分析（3）表面上的换热热阻远远大于肋片中的导热热阻，因而在任一截面上肋片温度可认为是均匀的（4）肋片顶端可视为绝热。

并不是扩展表面细长就可以按一维问题处理，必须满足上述四个假设才可视为一维问题。

2、肋片高度增加引起两种效果:肋效率下降及散热表面积增加。

因而有人认为随着肋片高度的增加会出现一个临界高度，超过这个高度后，肋片导热热流量会下降，试分析该观点的正确性。

答：的确肋片高度增加会导致肋效率下降及散热表面积增加，但是总的导热量是增加的，只是增加的部分的效率有所减低，所以我们要选择经济的肋片高度。

第三章：1、由导热微分方程可知，非稳态导热只与热扩散率有关,而与导热系数无关。

数值传热学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 对整个控制容积作能量平衡，有：2150)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中的相关项作局部线性化处理，然后迭代计算；求解结果：，（迭代精度为10-4）818.822=T 635.353=T 迭代计算的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)n gin th a r e 结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：习题4-14充分发展区的温度控制方程如下：)(1rTr r r x T uc p ∂∂∂∂=∂∂λρ对于三种无量纲定义、、进行分析如下w b w T T T T --=Θ∞∞--=ΘT T T T w ww T T T T --=Θ∞1）由得：wb wT T T T --=Θww 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 --=Θ∞ww T T T T +Θ-=∞)(由可得：T xT x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(])[(r T T r T T T r T 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)(θφλθφλφλφρθφρφρ、和分别是圆柱坐标的3个坐标轴，、和分别是其对应的速度分量，其中x r θ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=u ，；对称线上，0=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 θληλθηηλληηη上式化简得：011(1(1=+∂∂∂∂+∂∂∂∂θηθηηηηηU U 边界条件：，1=η0=U ，；对称线上，0=η0=∂∂ηU 0=∂∂θU定义无量纲温度：λ/0R q T T b-=Θ其中，是折算到管壁表面上的平均热流密度，即：；0q Rq q wπ=0由无量纲温度定义可得：bT Rq T +Θ=λ0将表达式和无量纲半径代入能量方程得：T η(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T uc b p 化简得：（1）)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p 由热平衡条件关系可以得：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 即得单值性条件：=Θ⎰⎰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∂∂Γ=∂∂φφρ边界条件如下：LL x x φφφφ====,;,00上述方程的精确解如下： 11)/(00--=--⋅Pe L x Pe L e e φφφφΓ=/uL Pe ρ2．将分成20等份，所以有：L ∆=P Pe 20 1 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 =i 2）一阶迎风中间节点： ∆-∆++++=P P i i i 2)1(11φφφ20,2 =i 3）混合格式当时，中间节点： 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 number tt=[1 5 10];%dimensionless length m=20;%mdim is the number of inner node mdim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculation y0=1;yL=2;%cal exact solution for n=1:size(tt,2) t=m*tt(1,n); if t==0 yval1(n,:)=eval(y1); else yval1(n,:)=eval(y); end end%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))==0 yval1(indexf(n),1)=yL; else yval1(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)-yval5com(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));a n d A l l t h i n g s i n t h ei r b e i n g a r e g 13精确解与数值解的对比图，其中边界条件给定，。

第四章复习题1、 试简要说明对导热问题进行有限差分数值计算的基本思想与步骤。

2、 试说明用热平衡法建立节点温度离散方程的基本思想。

3、 推导导热微分方程的步骤和过程与用热平衡法建立节点温度离散方程的过程十分相似，为什么前者得到的是精确描述，而后者解出的确实近似解。

4、 第三类边界条件边界节点的离散那方程，也可用将第三类边界条件表达式中的一阶导数用差分公式表示来建立。

试比较这样建立起来的离散方程与用热平衡建立起来的离散方程的异同与优劣。

5．对绝热边界条件的数值处理本章采用了哪些方法？试分析比较之．6．什么是非稳态导热问题的显示格式？什么是显示格式计算中的稳定性问题？7．用高斯－塞德尔迭代法求解代数方程时是否一定可以得到收敛德解？不能得出收敛的解时是否因为初场的假设不合适而造成？8．有人对一阶导数()()()221,253x t t t x t i n i n i n in ∆-+-≈∂∂++你能否判断这一表达式是否正确，为什么？ 一般性数值计算4-1、采用计算机进行数值计算不仅是求解偏微分方程的有力工具，而且对一些复杂的经验公式及用无穷级数表示的分析解，也常用计算机来获得数值结果。

试用数值方法对Bi=0.1,1,10的三种情况计算下列特征方程的根:)6,2,1( =n n μ3,2,1,tan ==n Binn μμ并用计算机查明，当2.02≥=δτa Fo 时用式（3-19）表示的级数的第一项代替整个级数（计算中用前六项之和来替代）可能引起的误差。

Bi n n =μμtanFo=0.2及0.24时计算结果的对比列于下表：第一项的值 0.94513 0.61108 0.10935 前六项的值 0.94688 0.6198 0.11117 比值 0.998140.986940.98364 Bi=0.1 Bi=1 Bi=10 第一项的值 0.99277 0.93698 0.77311 前六项和的值0.99101 0.92791 0.76851 比值1.001771.009781.005984-2、试用数值计算证实，对方程组⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=++=++=-+5223122321321321x x x x x x x x x用高斯-赛德尔迭代法求解，其结果是发散的，并分析其原因。

主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER 2010年9月27日, 西安数值传热学第四章扩散方程的数值解及其应用(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调和平均已经广泛为国内外学术界所接受。

≤1数学上的稳定未必导致物理上有意义的解无内热源一维非稳态导热，初场均匀，两表面0]T +代入下式：P（全隐格式）才能满足。

结论：数学上的稳定未必导致物理上有意义的解；推=xΔa TP P极坐标均可以表示成为：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 =−Δ∑代数方程迭代求解收敛的充分条件是，因为可以确保代数方程迭代求解收敛。

第四章 传热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）即得。

主讲陶文铨西安交通大学能源与动力工程学院热流中心CFD-NHT-EHT CENTER 2010年9月27日, 西安数值传热学第四章扩散方程的数值解及其应用(2)4.4 求解代数方程的TDMA 及ADI 方法4.4.1求解一维导热问题代数方程的三对角阵算法4.4.2求解多维非稳态导热全隐格式的ADI 方法1.求解方法概述1.一维导热问题代数方程通用形式2.Peaceman-Rachford 的ADI 迭代2.Thomas 算法3.第一类边界条件的处理4.4 方法4.4.1 求解一维导热问题代数方程的三对角阵算法Thomas算法的一般形式将上式改写为：所谓消元就是要找出系数间的关系。

(b)乘以CiP=11;i−−1,i=B D3.第一类边界条件下Thomas 算法的实施第一类边界条件下，求解区域为i=2,….M1-1=M2。

将消元公式用于i=1, 注意T 1是给定的：1121T PT Q =+10;P =11Q T =因T M1已知，消元从M 2开始：2212M M M T P T Q =+注意：采用附加源项法来处理第二类，第三类边界条件时，均将第二类，第三类边界条件问题视为第一类边界条件问题，数学上的处理与此相同。

4.4.2 求解多维导热问题代数方程的方法求解二维非稳态导热全隐格式代数方程的方法(1) 五对角阵算法(Penta-diagonal ,PDMA)(2) 交替方向隐式方法(Alternative-directionImplicit, ADI)2. 3-D Peaceman-Rachford方法将tΔ三等分：tΔX方向为隐式，第一个/3y,z方向为显式式；第二，三个/3tΔ分别在y,z方向实施隐式;2－D交替方向隐式u i,j,k , v 表示方向二阶导数的中心差分；nT用von Neumann分析方法可以证明稳定性条件为：4.5 管道内充分发展对流换热概说4.5.1管道内充分发展对流换热的定义1. 简单的充分发展对流换热2. 复杂的充分发展对流换热4.5.2能实现充分发展对流换热的边界条件4.5.3部分算例汇总4.5 管道内充分发展对流换热概说4.5.1管道内充分发展对流换热的定义平直通道中的充分发展对流换热属于这一类。