二维肋片稳态导热问题的数值计算

传热学：第四章 导热问题的数值求解

例题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处，而是在n>1的不远处的离散点上，这是和我们的预期是相违背的，但是当网格划分到达一定的密度，就可以避免这个问题，虽然在图像上看不出来了，但是这个问题还是存在的，不过由于足够小的网格，是它可以忽略。

(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>=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

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) contour(T) end

总结：在编本题程序时，η值总是比书上的参考值小一半左右，但是反复核对和查找错误都没有确定症结所在，最终从对数据和图像的分析来看，我断定程序的主体是没有问题的，从而把计算η值的公式作为重点核对对象，终于发现是分子上的一个括号放错了位置，从而导致这个巨大的偏差，可见程序的完全正确和公式的正确需要非常仔细的核对，当然在结果的数据和图表的分析上，可以帮助锁定程序的错误发生的范围，从而更好的完善程序。

传热学：第四章 导热问题的数值求解

例题4-6 无限大平板的一维非稳态导热问题数值计算 （1）自主编程，编程语言自定，最后提交源程序 （2）提交电子报告（word格式），包括：

（a）给出空间离散示意图（网格划分） （b）节点离散方程（显示、隐式皆可）

（c） 图示温度分布（可以利用origin或matlab） （d） 分析空间步长和时间步长对计算结果的影响 解：（a）空间离散示意图

（b）节点离散方程

t(m)=Fo*(t(m+1)+t(m-1))+(1-2*Fo)*t(m) t(1)=2*Fo*t(2)+(1-2*Fo)*t(1)

t(M)=t(M)*(1-2*Fo*Bi-2*Fo)+2*Fo*t(M-1)+2*Fo*Bi*tf （c） 温度分布

传热学：第四章 导热问题的数值求解

平板中温度的瞬态分布Δτ=0.18，

M=40

平板中温度的瞬态分布Δτ=0.18，M=10

传热学：第四章 导热问题的数值求解

平板中温度的瞬态分布Δτ=0.018，

M=40

平板中温度的瞬态分布Δτ=0.018，M=10

传热学：第四章 导热问题的数值求解

平板中温度的瞬态分布Δτ=1.8的图由于Fo已经大于0.5，无法满足程序要求，故没有给出图像。

（d） 分析空间步长和时间步长对计算结果的影响

从两张图的比较看来，Δτ=0.018时的温度变化更大，猜想空间步长和时间步长的作用是一样的，空间步长大使温度变化不明显，时间步长大的话也使温度变化不明显。故若使结果精确就要使空间步长和时间步长都要足够的小，才能保证精度。

(E)用matlab编程，源程序

function example46 step2=0.18; h=1163; k=50;

a=1.39e-5; x1=0.1; tf=300; t0=80; M=40;

step1=x1/(M-1); x=1; min=30; while min>0

T=min*60/step2; Bi=h*step1/k;

Fo=a*step2/step1^2; if Fo<=0.5 for m=1:M t(m)=t0; end

for n=1:T

for m=2:(M-1)

t(m)=Fo*(t(m+1)+t(m-1))+(1-2*Fo)*t(m); end