数值传热_传热学作业_matlab

合集下载

数值传热学部分习题答案

习题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(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT w，该无量纲温度定义是可以用分离变量法的；习题4－181）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：r r x x w r v r r r u x ∂∂+∂∂∂∂=∂∂+∂∂+∂∂1)()(1)(1)(φλφρθφρφρx 、r 和θ分别是圆柱坐标的3个坐标轴，u 、v 和w 管内的流动方向；对于管内的层流充分发展有：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 φ。

(完整版)传热学MATLAB温度分布大作业完整版

东南大学能源与环境学院课程作业报告作业名称：传热学大作业——利用matlab 程序解决热传导问题院系：能源与环境学院专业：建筑环境与设备工程学号：姓名：2014年11月9日一、题目及要求1. 原始题目及要求2. 各节点的离散化的代数方程3. 源程序4. 不同初值时的收敛快慢5. 上下边界的热流量(入=1W/(m C))6. 计算结果的等温线图7. 计算小结题目：已知条件如下图所示：200 °Ch=10W/二、各节点的离散化的代数方程各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24;100'C绝热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)>=epsx0=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;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746 【方法2 】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4; t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4; t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4; t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4; t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4; t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4; t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endco ntour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000100.0000 136.8905 146.9674 149.8587150.744 4100.0000 102.3012 103.2880 103.8632104.349 6100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117Jacobi 迭代程序】函数文件为：function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*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;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6);xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746四、不同初值时的收敛快慢1、 [方法1]在Gauss 迭代和Jacobi 迭代中，本程序应用的收敛条件均为norm(y-xO)>=eps即使前后所求误差达到 e 的-6次方时，跳出循环得出结果。

数值传热_数值传热学大作业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%)

传热大作业-数值解法-清华-传热学

一维非稳态导热的数值解法一、导热问题数值解法的认识（一）背景所谓求解导热问题，就是对导热微分方程在规定的定解条件下的积分求解。

这样获得的解称为分析解。

近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 。

数值传热学习题答案

数值传热学习题答案数值传热学习题答案数值传热学是热力学的一个重要分支，主要研究热量在物质中传递的机理和规律。

在实际工程中，我们经常会遇到各种与传热有关的问题，通过数值计算可以得到准确的答案。

下面我将为大家提供一些数值传热学习题的答案，希望能够帮助大家更好地理解和应用这门学科。

1. 一个铝制热交换器的表面积为10平方米，其表面温度为100摄氏度，环境温度为20摄氏度。

已知铝的导热系数为200 W/(m·K)，求热交换器的传热速率。

答：根据传热定律，传热速率与传热面积、传热系数和温度差之间成正比。

传热速率 = 传热系数× 传热面积× 温度差。

将已知数据代入公式中，可得传热速率= 200 × 10 × (100 - 20) = 160,000 W。

2. 一个房间的尺寸为5米× 5米× 3米，墙壁和天花板的厚度为0.2米，墙壁和天花板的导热系数为0.5 W/(m·K)，室内温度为25摄氏度，室外温度为10摄氏度。

求房间的传热损失。

答：房间的传热损失可以通过计算墙壁和天花板的传热速率来得到。

墙壁和天花板的传热速率 = 传热系数× 传热面积× 温度差。

墙壁和天花板的传热面积 = 2 × (5 × 5) + 2 × (5 × 3) = 70平方米。

将已知数据代入公式中，可得墙壁和天花板的传热速率= 0.5 × 70 × (25 - 10) = 525 W。

因此，房间的传热损失为525瓦特。

3. 一个水箱的体积为1立方米，初始温度为20摄氏度，水的密度为1000千克/立方米，比热容为4186 J/(千克·摄氏度)，水箱的表面积为2平方米，表面温度为100摄氏度。

已知水的传热系数为0.6 W/(m^2·K)，求水箱内水的温度随时间的变化。

matlab传热计算程序

matlab传热计算程序
传热计算在工程学和科学领域中是一个重要的应用。

Matlab是一个功能强大的工程计算软件，可以用于传热计算。

在Matlab中，你可以使用各种方法来进行传热计算，比如有限元法、差分法、有限体积法等。

以下是一些常见的传热计算程序的示例：
1. 热传导方程求解，你可以编写一个Matlab程序来求解热传导方程，根据给定的边界条件和初始条件，使用差分法或有限元法来离散方程，并进行时间步进求解，得到温度场的分布。

2. 对流换热计算，对于流体内部的对流换热问题，你可以编写一个Matlab程序来求解Navier-Stokes方程和能量方程，结合有限体积法来进行流场和温度场的耦合求解。

3. 辐射换热计算，针对辐射换热问题，你可以编写一个Matlab程序来计算辐射传热，比如使用辐射传热方程和辐射传热模型，结合离散方法进行求解。

4. 传热系统优化，除了单一的传热计算，你还可以使用Matlab进行传热系统的优化设计，比如通过建立传热模型和耦合其
他工程模型，使用优化算法来寻找最优的传热系统设计参数。

总之，Matlab提供了丰富的工具和函数，可以用于传热计算的各个方面。

通过编写程序，你可以灵活地进行传热计算，并且可以根据具体的问题需求进行定制化的计算和分析。

希望这些信息对你有所帮助。

传热学数值计算大作业

数值计算大作业一、用数值方法求解尺度为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 tx y∂∂+=∂∂x=0，0<y<H ：()f th t t xλ∂-=-∂ 定解条件 x=H ，0<y<H ：t=t 2 y=0，0<x<H ：t=t1t 1t 2h ；t fq=1000 w/m 2y=H ，0<x<H ：tq yλ∂-=∂ 2、写出内部节点和边界节点的差分方程内部节点：()()1,,1,,1,,122220m n m n m nm n m n m n t t t t t t x y -+-+-+-++=∆∆左边界： (),1,,1,1,,,022m n m n m n m nm n m n f m n t t t t t t x x h y t t y y y xλλλ-++---∆∆∆-+++∆=∆∆∆右边界： t m,n =t 2上边界： 1,,1,,,1,022m n m n m n m nm n m n t t t t t t y y q x x x x yλλλ-+----∆∆∆+++∆=∆∆∆ 下边界： t m,n =t 13、求解过程利用matlab 编写程序进行求解，先在matlab 中列出各物理量，然后列出内部节点和边界节点的差分方程，用高斯-赛德尔迭代法计算之后用matlab 画图。

数值传热学部分习题问题详解2

习题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); end tcal=t习题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(])[(rT r T T r T T T r T w w w w ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂∞∞)1()(])[( 同2）分析可知，除了轴向及周向均匀热流const q w =的情况外，有0=∂∂rT 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 p S ∂∂-= 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等份，所以有：∆=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Λ=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 5－3乘方格式：⎪⎪⎩⎪⎪⎨⎧<-≤≤--+≤≤->=∆∆∆∆∆∆∆∆10,010,)1.01(100,)1.01(10,055P P P P P P P P D a e E当1.0=∆P 时有：951.0)1.01.01()1.01(55=⨯-=-=∆P D a eE因为：301.0/3)()()()()()(===Γ=Γ=∆eee e e e e e e P u x u u x D ρδρρδ 所以：5297.2830951.0951.0=⨯==e E D a由系数关系式∆=-P D a D a eEw W 可得： 53.3130)951.01.0()(=⨯+=⨯+=∆w eEW D D a P a 且： 205.01.010=⨯=∆∆=txa P p ρ 当采用隐式时1=f ，因此可得：0597.62253.315297.280=++=++=P W E P a fa fa a同理可得当10=∆P 时有：0=E a ，3=W a ，5=P a5－5二维稳态无源项的对流－扩散问题的控制方程：)()()()(yy x x y v x u ∂∂Γ∂∂+∂∂Γ∂∂=∂∂+∂∂φφφρφρφφ 对于一阶迎风、混合、乘方格式的通用离散方程：S S N N W W E E P P a a a a a φφφφφ+++=其中：[]0,)(e e e E F P A D a -+=∆ []0,)(w w w W F P A D a +=∆ []0,)(n n n N F P A D a -+=∆ []0,)(s s s S F P A D a +=∆5－71）QUICK 格式的界面值定义如下：⎪⎪⎩⎪⎪⎨⎧-+=-+=)36(81)36(81WW P W w W E P e φφφφφφφφ0>u 对（5－1）式dxdx d d dx u d )()(φφρΓ=积分可得： w e w e dxd dx d u u )()()()(φφφρφρΓ-Γ=-对流项采用QUICK 格式的界面插值，扩散项采用线性界面插值，对于0>u 及均分网格有：)]()([]))(36())(36[(81x x u u W P w P E e w WW P W e W E P ∆-Γ-∆-Γ=-+--+φφφφρφφφρφφφ 整理得：WW w W w e w E e e P w e w e u u u x u x x x u u φρφρρφρφρρ)(81])(43)(81[])(83[)]()(83)(43[-++∆Γ+-∆Γ=∆Γ+∆Γ+-上式即为QUICK 格式离散得到的离散方程；2）要分析QUICK 格式的稳定性，则应考虑非稳平流方程：xut ∂∂-=∂∂φφ 在t ∆时间间隔内对控制容积作积分：⎰⎰⎰⎰∆+∆+∂∂-=∂∂t t t e w e w tt tdxdt x u dtdx xφφ得：dt u dx tt tw e e wttt ⎰⎰∆+∆+--=-)()(φφφφφ随时间变化采用阶梯显式，随空间变化采用QUICK 格式得：t u x WW P W W E P tP t t P ∆+---+-=∆-∆+)]3636(81[)(φφφφφφφφ整理得：xu t ni n i n i n i ni n i ∆+-+-∆---++87332111φφφφφφ对于初始均匀零场，假设在),(n i 点有一个扰动n i ε；对1+i 点写出QUICK 格式的离散方程：xu tni n i n i n i n i n i ∆+-+-∆--+++++8733121111φφφφφφ可得：ni n i xt u εφ∆∆=++8711 对1-i 点分析可得：ni n i xt u εφ∆∆-=+-8311 由于扩散对扰动的传递恒为正，其值为ni x t ερ2∆Γ∆，所以根据符号不变原则有： 0)/)83(2≥∆Γ∆+∆∆-ni n i n i xt x t u εερε 整理得到QUICK 格式的稳定性条件为：38≤∆P 5－91）三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形式，所以定义如下：⎩⎨⎧<++=>++=00u c b a u c b a EEE P e W P E e φφφφφφφφ根据上述定义，在0>u 时对控制容积内的对流项作积分平均可得：])()([1)(11WW W P E e w w e c b c a b a xxdx x x φφφφφφφ--+-+∆=-∆=∂∂∆⎰由表2－1式可知三阶迎风格式的差分格式：xxni n i n i n i ni ∆+-+=∂∂--+1221264211,φφφφφ 由控制容积积分法得到的对流项离散格式应与Taylor 离散展开得到的离散格式具有相同的形式和精度，所以比较可得：61,65,31-===c b a所以三阶迎风格式的函数插值定义为：⎪⎪⎩⎪⎪⎨⎧<-+=>-+=06165310616531u u EE E P e W P E e φφφφφφφφ2）由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数表达式的定义在形式上显然是一致的；6－1二维直角坐标中不可压缩流体的连续方程及动量方程如下：⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂=∂∂+∂∂)3()()()()()()2()()()()()()1(0vu S y y v x x v y p y vv x vu t v S y y u x x u x p y uv x uu tu y v x u ηηρρρηηρρρ假设常粘性，则0==v u S S ；对公式（2）及（3）分别对y x ,求偏导得：⎪⎪⎩⎪⎪⎨⎧∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂⎪⎪⎭⎫⎝⎛∂∂∂∂+∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂33222233)()()()()()(y v x v y y p y y vv y x vu y t v y y u x x u x p x y uv x x uu x t u x ηηρρρηηρρρ 两式相加得并变换积分顺序有：⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂y v x u yy v x u xy p x p x v u x u v y v v y y v u y u v x u u x y v x u t 2222222222ηηρρ利用连续方程有：⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂2222y p x p x v u y v v y y u v x u u x ρ ⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂+∂∂+∂∂+∂∂∂∂22222222222y p xp y v x u y v x u y v x u x v y u ρ 最后即得：⎥⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂x v y u y v x u y p x p ρ222226－4假设5*=P p ，则有：5105*-=-=e u 5.3)05(7.0*=-⨯=n v由连续性条件有：s w n e v u v u +=+按SIMPLE 算法有：'''*5)(P E P e e e p p p d u u +-=-+= '''*7.05.3)(P n P n n n p p p d v v +=-+=将上两式代入连续性方程中有：20507.05.35''+=+++-P P p p计算得：06.42'=P p所以：06.4706.425'*=+=+=P P P p p p06.371006.47=-=-=E P e p p u 94.32)006.47(7.0)(7.0=-⨯=-=N P n p p v6－5假设250*3=p ，150*6=p ，所以各点的流量为：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=11)15040(1.020)150250(2.024)25010(1.04)270250(2.010)250275(4.0*****E DC B A Q Q Q Q Q 上述流量满足动量方程，但并不满足连续性方程，所以对流量修正：⎪⎪⎪⎩⎪⎪⎪⎨⎧-⨯+-=-⨯+=-⨯+-=-⨯+-=-⨯+=)(1.011)(2.020)(1.024)(2.04)(4.010'6'5'6'3'3'4'2'3'3'1p p Q p p Q p p Q p p Q p p Q ED C B A 对节点3作质量守恒有：B DC A Q Q Q Q +=+即得：)(2.04)(2.020)(1.024)(4.010'2'3'6'3'3'4'3'1p p p p p p p p -⨯+--⨯+=-⨯+--⨯+对节点3作质量守恒有：F E D Q Q Q =+即得：20)(1.011)(2.020'6'5'6'3=-⨯+--⨯+p p p p联立求解上两式有：70.48'3-=p ，13.69'6-=p修正后的压力为：3.20170.48250'3*33=-=+=p p p 87.8013.69150'6*66=-=+=p p p修正后的流量为：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=09.4)87.8040(1.009.24)87.803.201(2.013.19)3.20110(1.074.13)2703.201(2.048.29)3.201275(4.0ED C B A Q Q Q Q Q由)(76p p C Q F F -=。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

节点 i 温度 t(℃)
1 250.00
2 232.01
3 213.59
4 194.71
5 175.35
6 155.46
7 134.98
8 113.87
9 92.06
10 69.47
11 46.00
热流密度 q= -4.47E+04 w/m2
（6）温度分布图
10.一砖墙厚 240mm，内、外表面的表面传热系数分别为 6.0 w/(m2·k)和 15 w/(m2·k)，墙 3 体材料的导热系数λ=0.43 w/(m·k)，密度ρ=1668kg/m ，比热容 c=0.75kJ/(kg·K)，室内空气温度保持不变为 20℃，室外空气温度周期性变化，中午 12 点温度最高为 3℃，晚上 12 点温度最低为-9℃，试用数值计算方法确定内、外墙壁面温度在一天内的变化。解：（1）建立模型本题属于常物性，无热源的一维非稳态导热过程将平壁沿厚度方向（x 方向）划分为 N 个均匀相等的间距，时间从τ=0 开始，按Δτ分割为 k 段，令 i 表示内节点位置，k 表示 kΔτ时刻，节点布置如图所示
开始
输入 n, np, tm, Bia, Bib, Fo
n=10, tm=480, tfa=20, Bia=6*0.24/(n-1)/0.43, Bib=15*0.24/(n-1)/0.43; Fo=0.43/1668/750*86400/tm/(0.24/(n-1))^2
No 打印“不稳定”
Yes k=k+1 k=1
1
2
3
4
N
N+1
本题给出 N=10。（2）通过热流法建立离散方程 a) 内节点离散方程
L
P
R
对节点 P(i)所代表的微元体，在 x 方向上与节点 P 相邻的节点分别为 L(i-1)和 R(i+1)。由于节点之间的间距很小，可以认为相邻节点间的温度分布是线性的。节点 P 所代表的网格单元与它周围各网格单元之间的导热量可根据傅里叶定律直接写为:
k ha ( t k f − t1 ) − λ k t1k − t2 t k +1 − t1k ∆x = ρc 1 ∆x ∆τ 2
整理上式，可得
ha ∆x k k 1 ρ c∆x 2 k +1 k t −t + ( t f − t1 ) = 2 λ∆τ ( t1 − t1 ) λ
k 2 k 1
上式中
τ k x i i=n+1
0 i=1
本题给出 N=10。（2）通过热流法建立离散方程（采用显式差分格式） a) 内节点离散方程导热微分方程式为
∂t ∂ 2t =a 2 ∂τ ∂x
温度对时间的一阶导数，采用向前差分，则
(1)
tik +1 − tik ∂t = ∆τ ∂τ i ,k
外壁面温度(℃) -1.39 -1.86 -2.14 -2.30 -2.36 -2.32 -2.20 -2.01 -1.76 -1.47 -1.13 -0.76 -0.37 0.04 0.46 0.89 1.31 1.72 2.11 2.49 2.83 3.15 3.44 3.68
时间 12:30 13:00 13:30 14:00 14:30 15:00 15:30 16:00 16:30 17:00 17:30 18:00 18:30 19:00 19:30 20:00 20:30 21:00 21:30 22:00 22:30 23:00 23:30 0:00
Φ LP = λLP
其中
ti −1 − ti ×1 ∆x
Φ RP = λRP
ti +1 − ti ×1 ∆x
λRP = 43 + 0.08 × λLP
(ti + ti +1 ) 2 (t + t ) = 43 + 0.08 × i −1 i 2
对节点 P(i)所代表的微元写热平衡式，即可得节点 P(i)温度的离散方程
tfb=-3-6*cos((s-1)/tm*2*pi)
控制打印 NP No No k>tm
Yes
打印 tn, tw Yes
打印 Bia, Bib, Fo
停机
（4）编写程序本题使用 matlab 软件，所编写的程序如下： clear;clc; n=10 % x 节点数 tm=480 % 时间节点 tfa=20; % 室内温度 np=10; Bia=6*0.24/(n-1)/0.43; % 计算 Bi 及 Fo 准则数初始值 Bib=15*0.24/(n-1)/0.43; Fo=0.43/1668/750*86400/tm/(0.24/(n-1))^2; p=(2*Bib+2)*Fo % 显式差分格式的稳定性条件 P<1 t=ones(1,n); % 输入墙体内部温度的初始值，各节点均设为 1 tw=ones(1,tm+1); tn=ones(1,tm+1); if(p>1) % 迭代程序开始 disp('不稳定') else for s=1:tm+1 tfb=-3-6*cos((s-1)/tm*2*pi); % 室外温度 for k=2:n-1 t(1)=2*Fo*(t(2)+tfa*Bia)+(1-2*Fo-2*Bia*Fo)*t(1); t(n)=2*Fo*(t(n-1)+tfb*Bib)+(1-2*Fo-2*Bib*Fo)*t(n-1); t(k)=Fo*(t(k-1)+t(k+1))+(1-2*Fo)*t(k); end tn(s)=t(1); tw(s)=t(n); end % 迭代程序结束 Bia % 输出计算结果 Bib Fo for k=1:tm/np tnn(1,k)=tn(1,np*k); tww(1,k)=tw(1,np*k); end tnn tww w=0:tm; % 绘制温度变化图 plot(24*w/tm,tn(w+1),'-r',24*w/tm,tw(w+1)),grid axis([0,24,-10,22]) set(gca,'xtick',[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]) set(gca,'ytick',[-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16,18,20,22,24]); xlabel('时间'),ylabel('温度（℃）') legend('内壁面温度','外壁面温度',0) end
温度对 x 的二阶导数，采用中心差分表达式应为
(2)
∂ 2t tik−1 − 2tik + tik+1 = 2 ∆x 2 ∂x i ,k tik +1 − tik t k − 2tik + tik+1 = a i −1 ∆τ ∆x 2
将上式移项整理得
(3)
将式(3)和式(2)代入式(1)，就得到内节点(i, k)的节点离散方程
Φ LP + Φ RP = 0
λLP ti =
b)
ti −1 − ti t −t ×1 + λRP i +1 i × 1 = 0 ∆x ∆x λLP ti −1 + λLP ti +1 λLP + λRP
边界节点离散方程由于本题的壁面温度属于第一类边界条件，因此
t1 = 250 tn +1 = 46
（5）计算结果准则数初始值：Bia = 0.3721，Bib = 0.9302，Fo = 0.0870
时间内壁面温度(℃) 0:30 8.22 1:00 9.45 1:30 10.28 2:00 10.89 2:30 11.37 3:00 11.76 3:30 12.09 4:00 12.36 4:30 12.59 5:00 12.79 5:30 12.97 6:00 13.13 6:30 13.27 7:00 13.40 7:30 13.52 8:00 13.64 8:30 13.75 9:00 13.86 9:30 13.96 10:00 14.06 10:30 14.16 11:00 14.26 11:30 14.36 12:00 14.46 （6）温度变化图
5. 一厚度为 250mm 无限大平壁，其导热系数λ=43+0.08t w/(m· k)，平壁一侧温度为 250℃，另一侧温度为 46℃，试用数值方法确定平壁内的温度分布，并确定通过该平壁的热流密度。解：（1）建立模型本题属于非常物性，无热源的一维稳态导热过程将平壁沿厚度方向（x 方向）划分为 N 个均匀相等的间距。节点布置如图所示
输入 n, e, k, t1, tn+1, ti
N=10, t1=250, tn+1=46, ti=1, e=0.01, k=10000
迭代次数 k=0
tti=ti
No
k=k+1
Yes No k>m Yes
打印 m, ti, q
打印“不收敛”
停机
（4）编写程序本题使用 matlab 软件，所编写的程序如下： clear;clc; t=ones(1,11); % 设定各项初始值 q=ones(1,11); t(1)=250; t(11)=46; e=0.01; k=1; while 1 % 迭代程序 tt=t; for i=2:10 a=43+0.08*(t(i-1)+t(i))/2; b=43+0.08*(t(i)+t(i+1))/2; t(i)=(a*t(i-1)+b*t(i+1))/(a+b); end k=k+1; ttt=abs(t(5)-tt(5)); if(ttt<e) break; end end for i=2:11 % 计算热流密度 a=43+0.08*(t(i-1)+t(i))/2; q(i)=q(i-1)+a*(t(i)-t(i-1)); end k t q=q(i)/0.25 w=0:25:250;v=w/25+1;y=t(v); %绘制温度分布图 plot(w,y),xlabel('位置（mm）'),ylabel('温度（℃）') legend('平壁内的温度分布',0),grid （5）计算结果迭代次数 k=75 平壁温度分布见下表