星期六, 12月 23, 2006

星期四, 12月 14, 2006

吳子青:熱流分析之通用情況推廣 程式碼

function [T]=heat_platen(n,m,ti,to)
% Prog calculating heat transfer through a plate by 3x3.
% Inputs:
% ti,to : input & output temperature, C
% n,m : n x m mesh (Heat flow region)
% input : Heat flow input location [x,y] x=1~n,y=1~m
% output: Heat flow output location [x,y] x=1~n,y=1~m
% Outputs:
% temp:temperatures at each layer,C
% Example:
% T=heat_platen(20,-10)
inx=1;
iny=1;
outx=n;
outy=m;

A=zeros(n,m);
B=zeros(n,m);
C=zeros(n*m,1);
C((inx-1)*m+iny,1)=ti;
C((outx-1)*m+outy,1)=to;

for i=1:n
for j=1:m
if i==1i==n
A(i,j)=3;
elseif j==1j==m
A(i,j)=3;
else
A(i,j)=4;
end
if (i==1i==n)&(j==1j==m)
A(i,j)=2;
end
end
end
B(inx,iny)=1;
B(outx,outy)=1;
A=A+B;

TT=zeros(m*n,m*n);
for k=1:n*m
for i=1:n
for j=1:m
if k==(i-1)*m+j
TT(k,k)=A(i,j);
end
if i==1
TT((i-1)*m+j,i*m+j)=-1;
elseif i==n
TT((i-1)*m+j,(i-2)*m+j)=-1;
else
TT((i-1)*m+j,i*m+j)=-1;
TT((i-1)*m+j,(i-2)*m+j)=-1;
end
if j==1
TT((i-1)*m+j,(i-1)*m+j+1)=-1;
elseif j==m
TT((i-1)*m+j,(i-1)*m+j-1)=-1;
else
TT((i-1)*m+j,(i-1)*m+j+1)=-1;
TT((i-1)*m+j,(i-1)*m+j-1)=-1;
end

end
end
end
T=reshape(inv(TT)*C,n,m);
mesh (T); figure(gcf)

資料前處理