# 熱擴散（Thermal Diffusion)

## 簡介


在這個宇宙中，能量的分佈最終會趨於平衡。而熱能在一個空間中的分佈情況可以由溫度隨著時間變化的過程來描述。於是我們可以推導出**熱方程式**以及其變化型態：


<br>

其中alpha為熱擴散率，為物體内部熱量擴散速率的量度，也就是物體内部溫度趨向均匀的能力。而當我們將空間設為一個一個的格點時，可以將上述方程式簡化為


代表著每一點的溫度變化會受到相鄰點溫度的影響，以此類推，完成熱擴散的過程。

## 例題
我們假設有一個邊長為0.1公尺的正方形銅板，其初始溫度值如下圖所示，格點的距離設為0.0025公尺，銅板的熱擴散率設為0.000111。目標是要觀察其溫度隨時間的擴散情況和中心點達到定溫的時間。

<br>

1. 參數(epsilon為收斂判定值、r_x\r_y為該點的溫度變化率）

%參數：熱擴散係數、邊界條件
alpha=1.11e-4;
L=0.05;
H=0.05;
dx=0.0025;
dy=0.0025;
tmax=50;
dt=0.01;
epsilon=1e-4;
r_x=alpha*dt/dx^2;
r_y=alpha*dt/dy^2;

<br>
<br>

2. 首先建立格點和銅板的中心點。接著建立一個對應格點數的溫度矩陣，並設定邊界的溫度條件。

% 建立 x,y軸網格點
nx=uint32(L/dx+1);
ny=uint32(H/dy+1);
[X,Y]=meshgrid(linspace(0,L,nx),linspace(0,H,ny));
%中心點
ic=uint32((nx-1)/2+1);
jc=uint32((ny-1)/2+1);   
% 建立初始溫度：左右25度、下50度、上0度
T=10*ones(ny,nx);
T(:,1)=25;
T(:,end)=25;
T(1,:)=50;
T(end,:)=0;

<br>
<br>

3. 最後來跑**while**迴圈，這邊是分成垂直方向和平行方向的溫度擴散來做計算。特別注意每一個點除了會受到周圍四個臨點的溫度影響，**自己也會貢獻熱量給周圍四個臨點**，因此要記得減去自己造成的影響。作圖的部分則是利用**contourf**函數來話出溫度擴散輪廓圖、**scatter**函數來畫出中心點溫度圖。記得要設定**pause(秒數)** 來更新圖片，才會有動圖的效果。最後當中心點溫度收斂時結束迴圈。

% 迴圈設定
n=0; 
nmax=uint32(tmax/dt);
while n < nmax
    n=n+1;
    T_n=T;
    for j=2:ny-1
        for i=2:nx-1
            T(j,i)=T_n(j,i)+r_x*(T_n(j,i+1)-2*T_n(j,i)+T_n(j,i-1))...
                +r_y*(T_n(j+1,i)-2*T_n(j,i)+T_n(j-1,i));
        end
    end
    if uint16(n/50) == n/50      
        % 溫度擴散輪廓圖
        subplot(1,2,1),contourf(X,Y,T,30);
        colormap(jet)
        title(sprintf('Time = %g s',n*dt)),
        colorbar(),
        xlabel('x (m)'),ylabel('y (m)')
        axis('equal','tight'),
        % 中心點溫度隨時間之變化
        subplot(1,2,2),scatter(n*dt,T(jc,ic),'r.');
        grid on;
        xlim([0 tmax]),xlabel('t (s)'),ylabel('Tcenter')
        hold('on')
        pause(0.01)
    end
    % 確認有無收斂
    err=max(max(abs((T-T_n))));
    if err<=epsilon
        break
    end  
end

<br>
<br>

4. 最後的圖片會長這樣，完整的動圖請點選下方的網址觀看。


https://www.youtube.com/watch?v=5hC7xpFavBA

## 練習
請下載附件的程式碼並加以修改來模擬不同物體、溫度、邊界條件的熱擴散情形並詳細比較。