photo

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

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)    

資料前處理