二维肋片稳态导热问题的数值计算【DOC精选】.docVIP

例题4-5 二维肋片稳态导热问题的数值计算 （1）自主编程，编程语言自定，最后提交源程序 （2）提交电子报告（word格式），包括： （a）给出空间离散示意图（网格划分） （b）节点离散方程 （c） 图示温度等值线（可以利用origin或matlab） 解：（a）空间离散示意图(由origin的graph作图) （b）节点离散方程 由上图所示得各节点的节点离散方程 结点1：T(m,n)=0.25*(T(m+1,n)+T(m-1,n)+T(m,n+1)T(m,n-1)) 结点2：T(M,n)= 1/(4+2*Bi)*(T(M,n-1)+T(M,n+1)+2*T(M-1,n)) 结点3：T(m,N)=1/(4+2*Bi)*(T(m-1,N)+T(m+1,N)+2*T(m,N-1)) 结点4：T(M,N)= 1/(2+2*Bi)*(T(M-1,1)+T(M,2)) 结点5：T(m,1)=0.25*(T(m-1,1)+T(m+1,1)+2*T(m,2)) 结点6：T(M,1)= 1/(2+2*Bi)*(T(M,N-1)+T(M,2)) （c） 温度等值线 等温线图，工况1（Bi=0.01）η= 0.9670 温度与y轴分布图（工况1） 工况2（Bi=1）η=0.1910 温度与y轴分布图（工况2） 图像分析：四幅图的显示来看，结果是可信的。要是网格划分过松，就会发现，在肋板顶端的绝热边界上温度的分布是有问题的，温度的最高值并不是在半肋板顶端边界n=1处，而是在n1的不远处的离散点上，这是和我们的预期是相违背的，但是当网格划分到达一定的密度，就可以避免这个问题，虽然在图像上看不出来了，但是这个问题还是存在的，不过由于足够小的网格，是它可以忽略。 (D)用matlab编程，源程序 function example T0=input(T0=); Tf=input(Tf=); h=input(h=); k=input(k=); x=input(x=); H=input(H=); M=input(M=); st=H/(M-1); N=floor(x/st)+1; Bi=h*st/k; p=1; for m=1:(M) for n=1:(N) T(m,n)=0; end end for n=1:N T(1,n)=T0-Tf; end while p==1; p=0; for m=1:M;n=1:N; c(m,n)=T(m,n); end for m=2:M for n=1:N if (m=2mMn=2nN) T(m,n)=0.25*(T(m+1,n)+T(m-1,n)+T(m,n+1)+T(m,n-1)); elseif (m==Mn=2nN) T(M,n)=1/(4+2*Bi)*(T(M,n-1)+T(M,n+1)+2*T(M-1,n)); elseif (m==Mn==N) T(M,N)=1/(2+2*Bi)*(T(M,N-1)+T(M-1,N)); elseif (m=2mMn==N) T(m,N)=1/(4+2*Bi)*(T(m-1,N)+T(m+1,N)+2*T(m,N-1)); elseif (m=2mMn==1) T(m,1)=0.25*(T(m-1,1)+T(m+1,1)+2*T(m,2)); elseif (m==Mn==1); T(M,1)= 1/(2+2*Bi)*(T(M-1,N)+T(M,2)); end end end for m=1:M;n=1:N; if abs(c(m,n)-T(m,n))=1E-6; p=1; end end end T1=0; T2=0; for m=2:1:M T1=T1+T(m,N); end for n=2:1:(N-1) T2=T2+T(M,n); end Q=(0.5*(T(1,N)+T(M,1))+T1+T2)/(((M-1)+(N-1))*80); T=rot90(T+20); disp(Q) disp(Bi) disp(N) disp(T