数值传热学部分习题答案

习题4－2一维稳态导热问题的控制方程：022=+∂∂SxTλ依据本题给定条件，对节点2采用二阶精度的中心差分格式，h,Tf321节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程：节点1：1001=T节点2：1505105321−=+−TTT节点3：75432=+−TT求解结果：，852=T403=T对整个控制容积作能量平衡，有：02150)4020(15)(3=×−−×=Δ+−=Δ+xSTThxSqffB即：计算区域总体守恒要求满足习题4－5在4－2习题中，如果，则各节点离散方程如下：25.03)(10fTTh−×=节点1：1001=T节点2：1505105321−=+−TTT节点3：25.03325.032)20(4015])20(21[−×+=−×++−TTTT对于节点3中的相关项作局部线性化处理，然后迭代计算；求解结果：，818.822=T635.353=T（迭代精度为10-4）迭代计算的Matlab程序如下：x=30;x1=20;whileabs(x1-x)0.0001a=[100;5-105;0-11+2*(x-20)^(0.25)];b=[100;-150;15+40*(x-20)^(0.25)];t=a^(-1)*b;x1=x;x=t(3,1);1endtcal=t习题4-12的Matlab程序%代数方程形式AiTi=CiTi+1+BiTi-1+Dimdim=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);forn=1:mdimcoematrix(n,n)=A(1,n);ifn=2coematrix(n,n-1)=-1*B(1,n);endifnmdimcoematrix(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);forn=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);forn=(mdim-1):-1:1Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=[T;Tcal];%绘图比较给定T值和计算T值plot(Tcal,'r*')holdonplot(T)结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：2节点1节点2节点3字段4字段5字段6字段7字段8字段9字段10T的初始值1.46869391.1594949.53424416-.50680737-2.0679442-4.2476615-7.1232765-10.72954-15.03053-19.884531T的计算值1.46869391.1594949.53424416-.50680737-2.0679442-4.2476615-7.1232765-10.72954-15.03053-19.884531习题4-14充分发展区的温度控制方程如下：)(1rTrrrxTucp∂∂∂∂=∂∂λρ对于三种无量纲定义wbwTTTT−−=Θ、∞∞−−=ΘTTTTw、wwTTTT−−=Θ∞进行分析如下1）由wbwTTTT−−=Θ得：wwbTTTT+Θ−=)(由T可得：xTxTxTTTxTwbwwb∂∂Θ−+∂∂Θ=∂+Θ−∂=∂∂)1(])[(rTrTTrTTTrTwwbwwb∂∂Θ−+∂Θ∂−=∂+Θ−∂=∂∂)1()(])[(由与bTr无关、与Θx无关以及xT∂∂、rT∂∂的表达式可知，除了均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的；wT2）由∞∞−−=ΘTTTTw得：∞∞+Θ−=TTTTw)(由T可得：xTxTTTxTww∂∂Θ=∂+Θ−∂=∂∂∞∞])[(rTrTTrTTTrT∂∂Θ+∂Θ∂−=∂+Θ−∂=∂∂∞∞∞)(])[(由与bTr无关、Θ与x无关以及xT∂∂、rT∂∂的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流的情况外，有constqw=0=∂∂rTw，则该无量纲温度定义是可以用分离变量法的；3）由wwTTTT−−=Θ∞得：wwTTTT+Θ−=∞)(由T可得：xTxTTTxT∂∂Θ−=∂+Θ−∂=∂∂∞)1(])[(3rTrTTrTTTrT∂∂Θ−+∂Θ∂−=∂+Θ−∂=∂∂∞∞)1()(])[(同2）分析可知，除了轴向及周向均匀热流constqw=的情况外，有0=∂∂rTw，该无量纲温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：Srrrrrrxxwrvrrrux+∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂+∂∂+∂∂)(1)(1)()(1)(1)(θφλθφλφλφρθφρφρx、r和θ分别是圆柱坐标的3个坐标轴，u、和分别是其对应的速度分量，其中vwx是管内的流动方向；对于管内的层流充分发展有：0=v、，0=w0=∂∂xu；并且x方向的源项：xpS∂∂−=Rr方向的源项：rpS∂∂−=θ方向的源项：θ∂∂−=prS1由以上分析可得到圆柱坐标下的动量方程：x方向：0)(1)(1=∂∂−∂∂∂∂+∂∂∂∂xpurrrurrrθλθλr方向：0=∂∂rpθ方向：0=∂∂θp边界条件：Rr=，0=u，0=r0=∂∂ru；对称线上，0=∂∂θu不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：)(1)(1θλθλρ∂∂∂∂+∂∂∂∂=∂∂TrrrTrrrxTucp边界条件：Rr=，wqrT=∂∂λ；0=r，0=∂∂rTπθ/0=，0=∂∂−θλT2）定义无量纲流速：Lq=0图4－244dxdpRuU2−=λ并定义无量纲半径：Rr/=η；将无量纲流速和无量纲半径代入x方向的动量方程得：0))1((1))1((122=∂∂−∂−∂∂∂+∂−∂∂∂xpUdxdpRRRRUdxdpRRRRθληλθηηλληηη上式化简得：01)1(1)(1=+∂∂∂∂+∂∂∂∂θηθηηηηηUU边界条件：1=η，0=U0=η，0=∂∂ηU；对称线上，0=∂∂θU定义无量纲温度：λ/0RqTTb−=Θ其中，是折算到管壁表面上的平均热流密度，即：0qRqqwπ=0；由无量纲温度定义可得：bTRqT+Θ=λ0将T表达式和无量纲半径η代入能量方程得：)(1)(100θληλθηηλληηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂RqRRRRqRRRxTucbp化简得：)1(1)(10θηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂xTucqRbp（1）由热平衡条件关系可以得：mmmbmpbppRUUqRuuRqAuudxdTAucxTucxTuc020221221)(===∂∂=∂∂ππρρρ将上式代入式（1）可得：5)1(1)(12θηθηηηηη∂Θ∂∂∂+∂Θ∂∂∂=mUU边界条件：0=η，0=∂Θ∂η；1=η，Rqqwπη10==∂Θ∂0=θ，0=∂Θ∂θ；πθ=，0=∂Θ∂θ单值条件：由定义可知：0/0=−=ΘλRqTTbbb且：∫∫Θ=ΘAAbUdAUdA即得单值性条件：0=Θ∫∫AAUdAUdA3）由阻力系数及定义有：fRe228)(21/Re⎟⎠⎞⎜⎝⎛=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡−=DDUDuudxdpDfememmeνρ且：mWbmWbmWRqTTDTTqNu,0,,0~2)/(2Θ=−=−=λλ65－21．一维稳态无源项的对流－扩散方程如下所示：xxu22∂∂Γ=∂∂φφρ（取常物性）边界条件如下：LLxxφφφφ====,;,00上述方程的精确解如下：11)/(00−−=−−⋅PeLxPeLeeφφφφΓ=/uLPeρ2．将L分成20等份，所以有：Δ=PPe20图示如下：123456………………………1718192021对于中心差分、一阶迎风、混合格式和QUICK格式分别分析如下：1）中心差分中间节点：2)5.01()5.01(11−Δ+Δ++−=iiiPPφφφ20,2=i2）一阶迎风中间节点：Δ−Δ++++=PPiii2)1(11φφφ20,2=i3）混合格式当时，中间节点：1=ΔP2)5.01()5.01(11−Δ+Δ++−=iiiPPφφφ20,2=i当时，中间节点：10,5=ΔP1−=iiφφ20,2=i4）QUICK格式*12111)35(8122121⎥⎦⎤⎢⎣⎡−−−++++++=+−−ΔΔ−ΔΔ+ΔiiiiiiiPPPPPφφφφφφφ2≠i*1111)336(8122121⎥⎦⎤⎢⎣⎡−−++++++=+−ΔΔ−ΔΔ+ΔiiiiiiPPPPPφφφφφφ2=i7数值计算结果与精确解的计算程序如下：%exceptforHS,anyotherschemedoesnttakePe0intoconsideration%expressionofexactsolutiony=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'))%inthecaseofPe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%gridPenumbertt=[1510];%dimensionlesslengthm=20;%mdimisthenumberofinnernodemdim=m-1;X=linspace(0,1,m+1);%initialvalueofvariableduringcalculationy0=1;yL=2;%calexactsolutionforn=1:size(tt,2)t=m*tt(1,n);ift==0yval1(n,:)=eval(y1);elseyval1(n,:)=eval(y);endend%extratreatmentbecausemaxnumberinMATLABis10^308ifmax(isnan(yval1(:)))yval1=yval1';yval1=yval1(:);indexf=find(isnan(yval1));forn=1:size(indexf,1)ifrem(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';8end%CDsolutiond=zeros(size(tt,2),mdim);a=repmat([1],size(tt,2),mdim);forn=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*(