数值传热_传热学作业_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 -=。
matlab在传热学例题中的应用
matlab在传热学例题中的应用
MATLAB 是一种广泛使用的数学软件，可以在传热学等领域中用于数值计算和可视化。
以下是 MATLAB 在传热学例题中的应用:
1. 求解热传导方程
热传导方程是传热学中的基本方程之一，可以用于描述热量在固体表面上的传递。
MATLAB 可以用于求解热传导方程，例如可以使用Navier-Stokes 方程求解器来求解热传导方程。
2. 模拟热传导过程
通过使用 MATLAB 中的数值积分方法，可以模拟热传导过程，例如在求解热传导问题时，可以使用有限差分法 (FDM) 来求解热传导问题。
3. 可视化传热过程
MATLAB 可以用于可视化传热过程，例如可以使用 MATLAB 中的图像处理工具箱来绘制热传导过程的可视化图像。
此外，MATLAB 还可以用于制作动画，以展示传热过程的变化。
4. 研究热传导特性
通过使用 MATLAB 进行传热模拟，可以研究热传导特性，例如可以研究热传导率、热阻等特性。
此外，MATLAB 还可以用于研究热传导的非线性特性，例如可以使用非线性优化工具箱来求解最优热传导率。
MATLAB 在传热学中的应用非常广泛，可以帮助传热学者更好地理解和研究传热过程。
数值传热学习题集
简答题集锦1.流动与传热数值模拟的基本任务是什么？（把原来在时间域及空间域上连续的物理量的场，如速度场和压力场，用一系列有限个离散点上的变量值的集合来代替，通过一定的原则和方式建立起关于这些离散点上场变量之间关系的代数方程组，然后求解代数方程组获得场变量的近似值CFD可以看做是在流动基本方程（质量守恒方程飞动量守恒方程、能量守恒方程）控制下对流动的数值模拟。
通过这种数值模拟，我们可以得到极其复杂问题的流场内各个位置上的基本物理量（如速度、压力、温度、浓度等）的分布，以及这些物理量随时间的变化情况，确定旋涡分布特性、空化特性及脱流区等。
）2.数值模拟过程如何实现，主要步骤是那些？（建模、网格划分、坐标系、数学方程、求解、后处理）a.建立反映工程问题或物理过程本质的数学模型;b.选择与计算区域的边界相适应的坐标系；c.建立网格；d.建立离散方程；e.求解代数方程组；f.后处理，显示计算结果3.建立离散方程有哪些主要方法？比较说明各种方法的优缺点？（有限差分、有限体积、有限元、有限分析等）4什么叫控制方程？常见的控制方程有哪几个？各用在什么场合？5试写出控制方程的通用形式，并说明通用形式中各项的意义？（写明通式，以及各个方程中通式的表达形式）6推导x 方向的动量控制方程中的源项u S 的表达式。
由此证明当密度和黏度为常数时，u S 变为0。
X 方向N-S 方程：Mx S xw z u z x v y u y divu x u x x p Dt Du +∂∂+∂∂∂∂+∂∂+∂∂∂∂++∂∂∂∂+∂∂-=)][()]([)2(μμλμρ)()())()())())()()()()()][()]([)2(gradu div divu xz w y v x u x gradu div S divu xz w y v x u x S S divu xz w y v x u x gradu div S xw z x v y x u x z u z y u y x u x S xw z u z x v y u y divu x u x Mx u Mx Mx Mx μλμμλμλμμμμμμμμμμλμ+∂∂+∂∂+∂∂+∂∂∂∂=++∂∂+∂∂+∂∂+∂∂∂∂=+∂∂+∂∂+∂∂+∂∂∂∂+=+∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂∂∂=+∂∂+∂∂∂∂+∂∂+∂∂∂∂++∂∂∂∂（（（）（）（ 因为密度为常数，所以连续性方程：0=∂∂+∂∂+∂∂z w y v x u ρρρ 推 得：0=∂∂+∂∂+∂∂z w y v x u 所以：Su= 0)()=∂∂+∂∂+∂∂+∂∂∂∂divu x z w y v x u x λμ（7区域离散为分几种，说明各自的特点。
数值传热_传热学作业_matlab
开始
输入 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
ha ∆x ρ c∆x 2 1 = ，于是 = Bia ， λ λ ∆τ Fo 1 ( t1k +1 − t1k ) 2 Fo 1 2 Bia + 2
k k t2 − t1k + Bia ( t k f − t1 ) =
移项整理得
k k + Bia t k t1k +1 = 2 Fo ( t2 f ) + (1 − 2 Bia Fo − 2 Fo ) t1 ， Fo ≤
tik +1 =
a∆τ k a ∆τ t + tik+1 ) + 1 − 2 2 2 ( i −1 ∆x ∆x
k ti
上式可写作
tik +1 = Fo ( tik−1 + tik+1 ) + (1 − 2 Fo ) tik
Fo ≤
b)
1 2
边界节点离散方程 由于本题的壁面温度属于第三类边界条件，因此
（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） 温度变化图
数值传热学答案
例题3-3：算例的matlab程序clear%************一维非稳态热传导方程的数值解法与解析解对比程序**************** %****************************前处理（定义网格与节点）********************x=linspace(0,1,11);t=linspace(0,0.2,101);%****************************前处理（定义边界）************************** % 初值条件for i=1:11if i*0.1<=0.5T(i+1,1)=2*i*0.1;elseT(i+1,1)=2*(1-i*0.1);endend% 边界条件T(1,:)=0;T(11,:)=0;%****************************求解器（差分算法）*************************F=input('please input 网格fourier数：F=')%采用时间一阶向前差分，空间的的二阶中心差分的显式离散格式for j=2:101for i=2:10T(i,j)=F*(T(i+1,j-1)+T(i-1,j-1))+(1-2*F)*T(i,j-1);endend%****************************后处理（结果输出）*************************%输出x和T的关系的二维图形axis([0 1 0 0.6]);plot(x,T([1:11],41),':*r',x,T([1:11],21),':*b',x,T([1:11],11),':*k')hold on%**************************与解析解进行比较**************************%方程的解析解：T(x,t)=(8/pi^2)*exp(-((pi^2)*t))*sin(pi*x)%*******************************************************************axis([0 1 0 0.6]);y1=(8/pi^2)*exp(-(pi^2)*0.192)*sin(pi*x);% t=0.192y2=(8/pi^2)*exp((-pi^2)*0.096)*sin(pi*x);% t=0.096y3=(8/pi^2)*exp((-pi^2)*0.048)*sin(pi*x);% t=0.048xlabel('x');ylabel('T');title('一维非稳态热传导方程的数值解法与解析解');plot(x,y1,'r',x,y2,'b',x,y3,'k')3-9.试证明扩散项的中心差分格式具有守恒性。
数值传热学作业-第一章
1、二维非稳态导热微分方程：S YT X T t T p +∂∂+∂∂=)2222c （λδδρ。
对于时间步进（x 方向，y 方向）及空间而言，该方程为何种类型的方程？解： 将二维非稳态导热微分方程化为：0c 2222=+-∂∂+∂∂S tT Y T X T p δδρλλ （1）x 方向：0,0,a ===c b λ。
则：04b 2=-=∆ac ，所以该二维非稳态导热方程为抛物型方程。
（2）y 方向：0,0,a ===c b λ。
则：04b 2=-=∆ac ，所以该二维非稳态导热方程也为抛物型方程。
（3）对于空间而言，二维非稳态导热方程可知：,0,a b c λλ===则：2240b ac λ∆=-=-<，所以该二维非稳态导热方程为椭圆型方程。
2、（补充不可压、常物性的条件。
写出守恒型和非守恒型控制方程，并推导二者关系。
） 解：由题可知，该流体为不可压缩、常物性流体，而且是有内热源的二维问题。
守恒型控制方程： 质量守恒方程：0=∂∂+∂∂yv x u ； 由于流体自身条件，使得0==v u S S ，得动量守恒方程：()()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂-=∂∂+∂∂22221y u x u v x p y vu x uu ρ ()()⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂-=∂∂+∂∂22221y v x v v y p y vv x uv ρ ； 能量守恒方程：()()T S y T x T a y vT x uT +⎪⎪⎭⎫ ⎝⎛∂∂+∂∂=∂∂+∂∂2222 . 非守恒型控制方程：质量守恒方程：无非守恒型 动量守恒方程：⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂-=∂∂+∂∂22221y u x u v x p y u v x u u ρ ⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂-=∂∂+∂∂22221y v x v v y p y v v x v u ρ ；能量守恒方程：T S y T xT a y T v x T u +⎪⎪⎭⎫ ⎝⎛∂∂+∂∂=∂∂+∂∂2222流速及温度的边界条件：进口截面：c T in =，()y u u =，0=v ；在两平板界面上：0=∂∂yT ，0=v ，()y u u =; 出口截面：0=∂∂x T ，0=∂∂x u ，0=∂∂y v .。
热传导问题的MATLAB数值计算
收稿日期:2002-05-09.作者简介:李 灿(1968-),女,副教授;株州,湖南冶金职业技术学院冶金系(412000).¹Partial different ial equation toolbox user c s guide.T he M ath Works,Inc.,2000.热传导问题的M AT LAB 数值计算李 灿湖南冶金职业技术学院冶金系高彦栋 黄素逸华中科技大学能源与动力工程学院摘要:分析了应用M AT LAB 中PDE 工具箱解热传导问题的方法和步骤,编制了三个难以用解析方法求解的算例.采用有限元法求解导热偏微分方程,应用PDE 工具箱得到数值解.对适合圆柱坐标描述的问题,通过公式变化将其转换为能用PDE 工具箱求解的形式.算例表明,用M AT LAB 对复杂形状和有内热源的非稳态导热问题进行数值计算和图形处理是方便高效的.关 键 词:热传导;非稳态导热;M AT L AB;数值计算中图分类号:T K 124 文献标识码:A 文章编号:1671-4512(2002)09-0091-03许多工程问题需要确定物体内部的温度场或确定其内部温度到达某一限定值所需要的时间,因此研究导热问题特别是非稳态导热问题十分重要.目前非稳态导热问题的描述方程为多维非线性的偏微分方程,这些方程只在几何形状与边界条件都较简单的情况下才能求得理论解,而对于几何形状和边界条件复杂的情况多用数值解法,需借助于计算机将时间和空间坐标划分成数量巨大的网格才能得到较精确的数值解.本文应用M ATLAB 中PDE 工具箱,求解复杂边界条件的热传导问题.1 求解方法求解方法是基于数值解法中的有限元法[1],其基本原理是把计算区域划分成一系列的三角形单元,每个单元上取一个节点,选定一个形状函数(抛物线形或双曲线形),并通过单元中节点上的被求解变量值表示该函数.通过对控制方程作积分来获得离散方程.有限元法的最大优点是对不规则区域适应性好,故用MATLAB 方法求解的结果在边界上也较精确.对于适合圆柱坐标和球坐标描述的问题,通过对其热传导方程的变换,也能在MATLAB 中求解.应用MATLAB 的PDE (Partial Differential Equation)工具箱可以解如下四类偏微分方程¹-$#(c $u)+au =f ;d(5u /5t)-$#(c $u)+au =f ;d(52u/5t 2)-$#(c $u)+au =f ;-$#(c $u)+au =E du,(1)式中,u 为域8上的求解变量;E 为特征值;d,c,a,f 为常数或变量;t 为时间变量.前3个方程分别称为椭圆方程、抛物线方程和双曲线方程,第4个方程称为特征值方程.导热问题的通用微分方程可写成[2]Q c p (5u /5t )=$(K $u )+q v ,(2)式中,u 为求解变量,此处表示被求解物体内的温度;K 为导热系数;q v 为热源的发热率密度;Q 为密度;c p 为定压比热容.可以看出,式(1)和式(2)中的抛物线方程有着类似的形式.其中,求解变量为区域的温度,d 与Q c p ,c 与K ,f -au 与q v 可以一一对应.M ATLAB 中的PDE 工具箱定义了两类边界条件hu =r ;n #(c $u)+qu =g ,(3)式中,n 为垂直于边界的单位矢量;h ,r ,q 和g 为常量或与u 有关的变量.方程(3)中的第1个方程称为狄利克雷(Dirichlet)边界条件,第2个方程称为纽曼(Neumann)边界条件.可以看出,导热问题中的第一类边界与狄利克雷边界条件对应,第二类和第三类边界条件与纽曼边界条件对应.这些对应关系可以使用MATLAB 中的PDE 工具箱第30卷第9期 华 中 科 技 大 学 学 报(自然科学版) V ol.30 No.92002年 9月 J.Huazhong U niv.of Sci.&T ech.(Nature Science Edition)Sep.2002来求解.对一个导热问题的计算可以按图1的步骤进行.图1 M AT L AB 计算流程图2 算例2.1 三维非稳态无内热源的导热问题边长为0.5m,0.7m 和1.0m 的长方形钢锭,置于炉温u f =1200e 的加热炉内,计算5h 后钢锭的温度.已知钢锭的K =40.5W/(m #e ),A =0.722@10-5m 2/s ,u 0=25e ,钢锭与外界的对流换热系数h =348W/(m 2#e ).由对应关系可得d =Q c p =K /A ;c =K ;f =0;a =0,边界条件为纽曼边界条件,且钢锭的6个边界条件均相同,由对应关系有:q =h ; g =hf .求得5h 后钢锭内部的温度分布如图2,温度梯度如图3.两图还显示了有限元求解的网格,图3底平面的箭头方向为热流密度方向.图2 5h时刻钢锭的温度分布云图图3 5h 时刻钢锭的温度梯度云图如果导热体物性系数K 为温度u 的函数,只要写出K (u)的函数关系式,就可以得到解.2.2 有内热源的圆柱体非稳态导热问题有一半径为0.2m,长为3m 的圆柱形核电站用燃烧棒置于u f =100e 的水中,由于链式反应,棒内有恒定的产热率密度q v =20000W/m 3,计算10h 后燃烧棒内的温度分布.已知,燃烧棒的密度Q =7800kg /m 3,c p =500W #s/(kg #e ),K =40W/(m #e ),u 0=0e ,燃烧棒右端恒温t r =100e ,左端有一恒热流q l =5000W/m 2,燃烧棒外表面与外界水的对流换热系数h =50W/(m 2#e ).此问题宜采用圆柱坐标,由于燃烧棒内温度沿半径对称分布,因此可以转换为(r ,z )坐标的二维问题.将圆柱坐标内的热传导方程改写为Q c p r 5u 5t -55r K r 5u 5r -55z K r 5u 5z =q v r ,(4)以使其形式与式(1)拟合.式(4)与式(1)中的抛物线方程对比可以得出:d =Q c p r ;c =K r ;a =0;f =q v r ,式中,z 对应第一个坐标方向(在直角坐标中为x 方向);r 对应第二个坐标方向(在直角坐标中为y ).燃烧棒左端的边界条件为:n #(K $u )=q r ,为纽曼边界条件,由对应关系得:q =0; g =q l r ,燃烧棒右端为狄利克雷边界条件u =100.燃烧棒上(外)边界条件n #(K $u)=h(u f -u)为纽曼边界条件,由对应关系得q =hr ;g =hu f r.解析域的下边界为棒的中心,其边界条件为n #(K $t)=0,也为纽曼边界条件.q =0,则把q 和g 都设为0即可.求得10h 时刻燃烧棒内部的温度分布如图4所示,热流密度分布如图5所示.图4 10h 时刻燃烧棒的温度分布云图92 华 中 科 技 大 学 学 报(自然科学版) 第30卷图5 10h 时刻燃烧棒的热流密度云图如果内热源是时间或空间的函数,写出函数关系式,也可以得到解.2.3 复杂边界的热传导问题考虑这样一个问题:一个正方形内嵌一菱形,其中正方形区域的密度为2W/m 2,导热系数为10W/(cm #e ),菱形的密度为1W/m 2,导热系数为5W/(cm #e ),并有内热源的发热率密度为10W/m 2,两个区域的定压比热容均为0.1J/(kg #e ).初始温度为0e .计算0.1s 之后的等温图如图6所示,箭头所指为热流方向,热流密度图如图7所示.由此可见,应用MATLAB 可以方便快捷地解出复杂几何形状和复杂边界条件的非稳态问题.并且其强大的图形可视化功能使计算结果形象、直观而便于理解.图6 0.1s时刻等温图图7 0.1s 时刻热流密度云图参考文献[1]陶文铨.数值传热学(第2版).西安:西安交通大学出版社,2001.[2]程俊国,张洪济.高等传热学.重庆:重庆大学出版社,1990.Numerical simulation of problems in heat conduction using MATLABL i Can Gao Yandong H uang S uyiAbstract:The method and steps for finding the solutions for problems in heat conduction w ith the PDE toolbox in M ATLAB are described.T hree ex amples difficult to resolve w ith the analy tical method are g iven.The partial differential equation (PDE)for heat conduction is solved w ith the finite element method and the PDE toolbox is adopted to obtain the num erical simulation.Problems suitable for description w ith cylindrical coordinates are transformed into forms that are capable of solution w ith the PDE toolbox through formula variation.Examples of calculation show that M ATLAB is convenient and highly efficient for numerical simulation and graphic processing of com plex g eometry and non -steady -state heat conduction problems w ith internal thermal source.Key words:heat conduction;non -steady state heat conduction;MATLAB;numerical simulationLi Can Assoc.Prof.;Dept.of M etallurgy,Hu c nan Metallurg y Professional and Technical College,Zhuzhou 412000,Hu c nan,China.93第9期 李 灿等:热传导问题的M AT LAB 数值计算。
传热学数值计算大作业
传热学数值计算大作业传热学是研究物体内部和之间热量传递的科学，其应用范围广泛，例如在工程领域中，传热学的数值计算被广泛用于优化热传递过程，提高能源利用效率。
本文将介绍传热学数值计算的大作业，主要内容包括问题陈述、计算方法和结果分析等。
问题陈述：本次大作业的问题是研究一个热管的热传递特性。
具体来说，热管由内外两个半圆形的金属管组成，内管壁与外管壁之间是一种导热的传热介质。
问题要求计算热管内外壁的温度分布，并分析传热过程的效率和优化热管的设计。
计算方法：计算热传递过程需要运用一些热传导定律和传热方程。
首先，根据Fourier 热传导定律，可得到内外壁的温度梯度。
然后，使用热传导方程来描述热传递过程，其中包括热扩散项和传热源项。
在计算热传导时需要注意材料的热导率、导热介质的热传导性质等参数。
在计算中，可以使用一些数值方法来离散化热传导方程，例如有限差分法、有限元法等。
其中，有限差分法是一种常见的数值方法。
通过将热传导方程中的导数用差分表达式替代，可以将偏微分方程转化为代数方程。
然后，可以使用迭代方法求解代数方程，得到温度分布的数值解。
结果分析：通过数值计算，可以得到热管内外壁的温度分布。
根据温度分布，可以分析热传递过程中的热流分布和传热效果。
例如，可以计算内外壁之间的热传导率，评估热管的热传递效率。
同时，可以对热管的设计进行优化。
例如，可以通过改变热导率高低、加大导热介质的厚度等方式，来提高热传递效果。
此外，对于热管的材料选择和导热介质的设计，还可以进行参数敏感性分析。
通过改变各个参数的数值，可以研究其对热传递过程的影响程度。
这有助于优化热管的设计，并提供一些实际应用方面的建议。
总结：传热学的数值计算是研究热传递现象的重要工具，可以帮助我们深入了解传热过程，优化传热装置的设计。
通过本次大作业，我们可以学习和练习传热学数值计算的方法和技巧，提升对传热现象的理解和分析能力。
希望通过这次大作业，能够更好地应用所学知识，解决实际问题。
哈尔滨工程大学传热学大作业数值计算matlab程序内容
传热学作业数值计算1 程序内容： 数值计算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; =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= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (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 =tn; end k 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 数值解温度场 数值解每次迭代的最大误差max1 =9.8531e-07 解析温度场 t2 解析温度场取第11行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差数值计算matlab 程序内容：程序内容： >> tw1=10; tw2=20; c=1.5; p2=20; p1=c*p2; L2=20; deltaX=L2/p2; L1=c*L2; a=p2+1; b=p1+1; =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= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (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)=tw2; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 tx=ones(a,b); for i=1:1:a for j=1:1:b y=(a-i)*deltaX; x=(j-1)*deltaX; m=sym('m'); g=(((-1)^(m+1)+1)/m)*sin(m*pi*x/L1)*sinh(m*pi*y/L1)/sinh(m*pi*L2/L1); h=symsum(g,m,1,100); tx(i,j)=2*h*(tw2-tw1)/pi+tw1; end end tx 迭代次数k = 695 数值解温度场 数值解每次迭代的最大误差max1 =9.8243e-07 解析温度场 tx = 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像： 数值温度场 图像图像 ：解析温度场tx图像：数值解与解析解的误差数值解与解析解的误差程序内容： 数值计算matlab程序内容：>> t0=90; =10; L=10; c=0.25; p2=20; p1=p2/c; B=c*L; d=0.5*B; h=10; a=p2+1; b=p1+1; deltaX=B/p2; lambda=160; Bi=h*deltaX/lambda; =ones(a,b)*10; m1=ones(a,b)*3; m1(2:a-1,1)=zeros(a-2,1); m1(a,2:b-1)=ones(1,b-2); m1(1,2:b-1)=ones(1,b-2)*6; m1(2:a-1,b)=ones(a-2,1)*2; m1(1,b)=ones(1,1)*4; m1(a,b)=ones(1,1)*5; m1(1,1)=7; m1(a,1)=8; tn= ; max1=1.0; k=0; while ( max1>1e-6) k=k+1; max1=0; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=t0; case 1 tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 2 tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4* )/(4+2*Bi)+ ; case 3 tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+1)+tn(i-1,j)+tn(i+1,j)); case 4 tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2* )/(2*Bi+2)+ ; case 5 tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2* )/(2*Bi+2)+ ; case 6 tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 7 tn(i,j)=t0; case 8 tn(i,j)=t0; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k ta=ones(a,b); Bi1=h*d/lambda; sbi=sqrt(Bi1); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end x=deltaX*(j-1); ta(i,j)=(cosh(sbi*(L-x)/d)+sbi*sinh(sbi*(L-x)/d))*(t0- )/(cosh(sbi*L/d)+sbi*sinh(sbi*L/d))+ ; end end ta 迭代次数k =1461 数值解温度场 解析温度场 ta 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像如下图像如下数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差程序内容：数值计算matlab程序内容：>> tw=10; L2=15; c=0.75; L1=L2/c; p2=24 ; p1=p2/c; deltaX=2*L2/p2; a=p2+1; b=p1+1; lambda=16; qv0=24; =ones(a,b)*5; m1=ones(a,b); m1(1,:)=zeros(1,b); m1(2:a,b)=zeros(a-1,1); m1(2:a,1)=zeros(a-1,1); m1(a,2:b-1)=zeros(1,b-2); tn= ; max1=1.0; k=0; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=tw; case 1 tn(i,j)=0.25*(tn(i-1,j)+tn(i+1,j)+tn(i,j-1)+tn(i,j+1)+qv0*(deltaX^2)/lambda); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k; tx=ones(a,b); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end if j>(b+1)/2 x=(j-(b+1)/2)*deltaX; else x=-((b+1)/2-j)*deltaX; end m=sym('m'); xi=(2*m-1)*pi/2; g=((-1)^m)/(xi^3)*(cosh(xi*y/L1)/cosh(xi*L2/L1))*cos(xi*x/L1); h=symsum(g,m,1,100); tx(i,j)=2*qv0*L1^2/lambda*h+qv0*(L1^2-x^2)/(2*lambda)+tw; end end tx 数值温度场 解析温度场tx 取第13行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点曲线为第13行的解析解的直线，散点为其数值解的点解析解 第13行的误差=[数值解(13行) –解析解(13行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差。
Matlab在热传递学课程中的应用
Matlab在热传递学课程中的应用热传递学是研究热能传递和传导的学科，广泛应用于工程、物理和环境领域。
在热传递学课程中，Matlab是一个常用的工具，可以帮助学生理解和分析热传递过程。
下面将介绍在热传递学课程中使用Matlab的步骤和应用。
第一步是建立热传递模型。
在研究热传递过程时，我们需要建立相应的数学模型。
可以使用Matlab来编写这些模型，并通过求解数学方程来分析热传递现象。
例如，我们可以使用Matlab编写热传导方程，并求解得到温度分布。
第二步是处理边界条件。
在热传递过程中，边界条件对结果有着重要的影响。
例如，我们可以设置材料的初始温度、表面的热通量或边界温度等。
Matlab提供了丰富的边界条件处理函数和图形界面，使得处理边界条件变得更加简便。
第三步是求解热传递问题。
在建立了合适的模型和边界条件后，我们可以使用Matlab的数值求解方法来求解热传递问题。
Matlab提供了许多数值求解算法，如有限差分法和有限元法，可以帮助我们得到准确的结果。
通过对求解结果的分析和可视化，我们可以更好地理解热传递过程。
第四步是进行参数敏感性分析。
在研究热传递过程时，我们通常需要考虑不同的参数对结果的影响。
Matlab提供了参数敏感性分析的工具，可以帮助我们理解不同参数对热传递问题的影响程度。
通过参数敏感性分析，我们可以选择最优的参数组合，并优化热传递系统的设计。
第五步是进行热传递实验和数据处理。
除了数值分析，实验也是研究热传递的重要手段。
Matlab 可以辅助我们进行热传递实验的数据处理和分析。
通过编写Matlab程序，我们可以快速地进行数据处理、绘图和拟合曲线，从而更好地理解实验数据和验证理论模型。
综上所述，Matlab在热传递学课程中具有广泛的应用。
它可以帮助学生建立热传递模型，处理边界条件，求解热传递问题，进行参数敏感性分析，并辅助实验数据处理。
通过使用Matlab，学生可以更好地理解和分析热传递过程，提高问题解决能力。
matlab差分法求解传热微分方程
matlab差分法求解传热微分方程matlab差分法求解传热微分方程是一种利用有限差分法对传热微分方程进行数值求解的方法。
有限差分法是一种广泛应用于偏微分方程求解的数值方法，它将连续的空间和时间离散化，通过离散的节点上的值来逼近连续的解。
1.matlab差分法求解传热微分方程求解步骤与规则：在使用matlab差分法求解传热微分方程时，通常遵循以下步骤和规则：（1）确定问题的物理模型：首先明确传热问题的类型，如稳态或非稳态传热，以及问题的边界条件和初始条件。
（2）离散化：将求解域划分为离散的网格节点，用于表示空间坐标。
同时，根据时间步长对时间进行离散化。
（3）建立差分方程：根据离散化的空间和时间，将传热微分方程转化为差分方程。
这通常包括对空间导数和时间导数的差分处理。
（4）选择适当的差分格式：根据问题的稳定性要求，选择合适的差分格式。
常见的有限差分格式有向前差分、向后差分、中心差分等。
（5）编写matlab程序：根据所选的差分格式和差分方程，编写matlab 程序进行求解。
程序中需要包含迭代算法，如三角剖分法、平方根法等，以收敛解为最终结果。
（6）验证和分析：通过与理论解或实验数据进行对比，验证求解结果的正确性和稳定性。
此外，分析结果以了解传热过程的规律和特点。
2.matlab差分法求解传热微分方程使用场景：matlab差分法求解传热微分方程在以下场景中具有广泛应用：（1）建筑节能分析：通过求解建筑墙体内部的温度场，评估建筑物的保温性能，为建筑设计提供依据。
（2）电子散热设计：分析电子设备内部的传热过程，优化散热结构，提高设备的工作效率和可靠性。
（3）能源工程：研究能源转换和传输过程中的热损失，为提高能源利用效率提供理论支持。
（4）材料科学：研究材料的热传导性能，为新材料的开发和应用提供参考。
（5）生物医学：分析生物组织内的传热过程，辅助诊断和治疗相关疾病。
总之，matlab差分法在求解传热微分方程方面具有较高的准确性和灵活性，为各种传热问题的分析和解决提供了有力工具。
传热学实验—-墙角matlab导热问题
二维导热物体温度场的数值模拟姓名小明学号 111111班级能动学院能动一、问题描述有一墙角模型，尺寸如图1所示，导热系数0.53W/(m·K)，墙角内外壁为第一类边界条件。
求解该模型的温度分布及导热量。
图1q=0二、计算原理根据热平衡法列出节点方程，各方向导入单元体的热量之和为零。
内节点和绝热边界点（图1点划线上的点）的方程形式不同。
图2 Array图2所示的内节点和绝热边界节点方程如下：内节点：)()()()(1,,1,,1,1,,1,=⎥⎦⎤⎢⎣⎡-+-+-+-••=+++-+-+x y t t x y t t y x t t y x t t j i j i j i j i j i j i j i j i W E S N ∆∆∆∆∆∆∆∆ΦΦΦΦλ绝热边界点：)(02)(2)(1,,1,1,,1,=⎥⎦⎤⎢⎣⎡-++-+-••=+++--+x y t t y xt t y x t t j i j i j i j i j i j i W E S N ∆∆∆∆∆∆ΦΦΦΦλ三、计算过程用Matlab7.1语言编写计算程序，初取网格步长m y x 1.0=∆=∆运行结果：图1：各个点的温度数值图2：分层设色等温线分布图3：等温线分布（每两条线间隔为三度）四、小结本次数值模拟是运用matlab程序用于数值计算。
小组成员共同讨论并复习了热传导问题的数学描述和热平衡法；从模拟过程中练习了不同节点迭代方程的建立；并简单学习了matlab语言的使用。
这次大作业对于我们以后的学习和可能的研究来说是一个很好的锻炼机会。
方腔内自然对流MATLAB程序数值传热学
(2)
其中， 是热扩散系数。 边界条件： 由无滑移边界条件得，四个壁面上的速度均为零，即： un us uw ue 0 。 在热壁面上 T 100 ，在冷壁面上 T 0 ，在上下绝热壁面上处
T 0。 y
数值处理：
区域离散化如图 2 所示。 对于动量守恒方程，在不考虑压力的情况下先计算出一个临时速度
%========================================================================== uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,1:ny+1)+u(1:nx+1,2:ny+2)); vv(1:nx+1,1:ny+1)=0.5*(v(1:nx+1,1:ny+1)+v(2:nx+2,1:ny+1));
%==========================================================================
un=0.0; us=0.0; vw=0.0; ve=0.0; Tw=100.0; Te=0.0; T(1,1:ny+1)=100; TT(1,1:ny+1)=100; %========================================================================== u(1:nx+1,1)=2*us-u(1:nx+1,2); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(1,1:ny+1)=2*vw-v(2,1:ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1);
- 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 平壁温度分布见下表