传热学导热问题的数值计算实践报告.docx

传热学
导热问题的数值计算
实践报告

姓名：
学号：
班级：
完成日期：2019年12月12日

一、实践题目及要求
例题4-5二维肋片稳态导热问题的数值计算
（1）自主编程，编程语言自定，最后提交源程序
（2）提交电子报告（word格式），包括：
（a）给出空间离散示意图（网格划分）
（b）节点离散方程
（c）图示温度等值线（可以利用origin或matlab）
二、空间离散示意图

三、节点离散方程
由上图所示得各节点的节点离散方程
节点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)]
四、温度等值线
1、等温线图，工况1（Bi=0.01）η=0.9670

2、工况1，温度与y轴分布图

3、等温线图，工况2（Bi=1）η=0.1910

4、温度与y轴分布图（工况2）

图像分析：从四幅图的显示来看，结果是可信的。要是网格划分过松，就会出现在肋板顶端的绝热边界上温度的分布存在问题，温度的最高值并不是在半肋板顶端边界n=1处，而是在n>1的不远处的离散点上，这与预期是相违背的，但是当网格划分到达一定的密度，就可以避免这个问题，虽然在图像上看不出来，但此问题还是存在的，不过由于网格足够小，可忽略。
五、Matlab编程源程序
functionexample
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;
form=1:(M)
forn=1:(N)
T(m,n)=0;
end
end
forn=1:N
T(1,n)=T0-Tf;
end
whilep==1;
p=0;
form=1:M;n=1:N;
c(m,n)=T(m,n);
end
form=2:M
forn=1:N
if(m>=2&&m<M&&n>=2&&n<N)
T(m,n)=0.25*(T(m+1,n)+T(m-1,n)+T(m,n+1)+T(m,n-1));
elseif(m==M&&n>=2&&n<N)
T(M,n)=1/(4+2*Bi)*(T(M,n-1)+T(M,n+1)+2*T(M-1,n));
elseif(m==M&&n==N)
T(M,N)=1/(2+2*Bi)*(T(M,N-1)+T(M-1,N));
elseif(m>=2&&m<M&&n==N)
T(m,N)=1/(4+2*Bi)*(T(m-1,N)+T(m+1,N)+2*T(m,N-1));
elseif(m>=2&&m<M&&n==1)
T(m,1)=0.25*(T(m-1,1)+T(m+1,1)+2*T(m,2));
elseif(m==M&&n==1);
T(M,1)=1/(2+2*Bi)*(T(M-1,N)+T(M,2));
end
end
end
form=1:M;n=1:N;
ifabs(c(m,n)-T(m,n))>=1E-6;
p=1;
end
end
end
T1=0;
T2=0;
form=2:1:M
T1=T1+T(m,N);
end
forn=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)
contour(T)
end

六、个人总结与心得体会
在本次实践中，我取得了较大收获。不但重新巩固学习了已经许久未接触使用的MATLAB的基本操作，而且基本初步学习到了如何用编程的方法去求解实际导热数值计算问题，也加深了对课本导热理论知识的理解。而且在实践中也遇到问题并解决，让我明白编程时细心的重要性，程序的完全正确和公式的正确需要非常仔细的核对，当然在结果的数据和图表的分析上，可以帮助锁定程序的错误发生的范围，从而更好的完善程序。