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

前言：
上禮拜上課老師講解了10-2 範例五：簡易之熱流分析
並探討了2x2 和 3x3 的範例
我就想寫一個n x m 維通用Case的程式
因為探討在Steady state穩態熱流情況基本原理並不難

除了把探討的維度改成n x m 外，在這裡也改成可以自由選擇熱流輸入
和輸出的位置，也可以多點的輸入輸出。

注意的地方是Control Volume 控制體積仍然維持正方形，原因在上課講解的時候有解釋過了，而且這裡探討的是Steady state。


分析與流程：
這個問題是示範反矩陣解聯立的方法
所以重要的就是在於係數矩陣 (也就是聯立的係數) 怎麼寫
而C矩陣只和輸入輸出熱流項有關

首先介紹函式的用法
function [T]=heat_platen(n,m,ti,to,input,output)
輸入項 n 和m 代表熱流場的維度
Ti: 輸入溫度
To: 輸出溫度
Input： 熱流輸入的位置座標，以矩陣表示並可以多選。
Output：熱流輸出的位置座標

輸出項： T 溫度矩陣
並在輸入項加入了預設值的設計
也就是如果輸入沒有填熱流流出和流入的位置
程式將會預設成從( 1, 1 ) 傳入而傳到 ( n , m ) 熱流流出，也就是範例五的樣式

接下來簡略的說明想法
把2x2的係數矩陣和3x3的係數矩陣列出來

2x2 case :

3 -1 -1 0
-1 2 0 -1
-1 0 2 -1
-1 -1 0 3

3x3 case

3 -1 0 -1 0 0 0 0 0

-1 3 0 -1 -1 0 0 0 0

0 -1 2 0 0 -1 0 0 0

-1 0 0 3 -1 0 -1 0 0

0 -1 0 -1 4 -1 0 -1 0

0 0 -1 0 -1 3 0 0 -1

0 0 0 -1 0 0 2 -1 0

0 0 0 0 -1 0 -1 3 -1

0 0 0 0 0 -1 0 -1 3

可以觀察到有一定的規則性
在這裡分成矩陣的對角線元素的探討 A 和其餘部份的探討 B
然後最後再組合起來，最後在加上輸入輸出項

(1) 對角線元素為相鄰的區塊數目，如果是在四個角落為2，在矩形的四個邊
為3，而其他的區域為4，由這樣的規則成矩陣A
(程式 line 27 到 line 46 )

(2) 如果探討的控制體積是第k塊且是位在中間區域 (即不是四角及邊線)，
則在k+1、k-1、k+m、k-m的位置都要填 -1，而在邊線部份或四個角就只要填部份三個或兩的位置為 -1 ，並要特別注意i= 1 or n 或 j= 1 or m 的情況
由這樣的規則形成矩陣B
(程式 line 51 到 line 77)

(3) 加上輸入和輸出的部份，在輸入和輸出位置對應的對角線位置分別 +1 ，而在矩陣C 填上輸入溫度Ti 和輸出溫度To

(4) 最後利用反矩陣計算聯立方程式的解 T= inv(A)*C T為溫度矩陣
也是最終 Steady state 的溫度分佈 ，並繪製成3D圖


程式內容

function [T]=heat_platen(n,m,ti,to,input,output)
% 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,20,20,[1 1],[20 20])

if nargin==5
    outx=n;outy=m;
    inx=input(:,1);
    iny=input(:,2);
elseif nargin==4
    outx=n;outy=m;
    inx=1;iny=1;
else
inx=input(:,1);
iny=input(:,2);
outx=output(:,1);
outy=output(:,2);
end

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==1|i==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)

結果

Case 1 改變熱流場Size
圖一5x5 Case
圖二10x10 Case
圖三50x50 Case

(輸入溫度 Ti = 20度 在 (1,1)
輸出溫度 To = -10 度在 (n,m))

























Case2 熱流入流出位置的改變
圖四 Input = [ 1 1] Output = [1 20]
圖五 Input = [ 1 1] Output = [10 10]
圖六 Input = [ 10 10] Output = [1 1]

(輸入溫度 Ti = 20度 輸出溫度 To = -10 度
Size= 20 x 20 )
























Case3 多個熱流出入口的改變
圖七 Input = [ 1 1; 20 1] Output = [1 20]
圖八 Input = [ 1 1; 20 20] Output = [1 20 ; 20 1]
圖九 Input = [ 1 1 ; 1 20 ; 20 1 ; 20 20] Output = [10 10]
圖十 Input = [ 10 10] Output = [ 1 1 ; 1 20 ; 20 1 ; 20 20]