當前位置:網站首頁>國賽助力:第三類邊界條件熱傳導方程及基於三對角矩陣的數值計算MATLAB實現(2020A)

國賽助力:第三類邊界條件熱傳導方程及基於三對角矩陣的數值計算MATLAB實現(2020A)

2021-08-20 06:35:16 slandarer

1第三類邊界條件的熱傳導方程

1.1 熱傳導方程
熱傳導在一維的各向同性介質裏的傳播可用以下方程錶達:
∂ u ∂ t = a ∂ 2 u ∂ x 2 (1) \frac{\partial u}{\partial t}=a \frac{\partial^{2} u}{\partial x^{2}} \tag{1} tu=ax22u(1)

其中, u = u ( x , t ) u=u(x,t) u=u(x,t) a = λ c ρ a=\frac{\lambda}{c\rho} a=cρλ λ \lambda λ錶示介質的熱傳導率, c c c錶示介質的比熱, ρ \rho ρ錶示介質的密度。
.
1.2 第三類邊界條件
考察介質放在另一種介質中的情形。外界介質的溫度 U U U與所考察介質錶面上的溫度 u u u往往並不相同,考慮流過所考察介質錶面的熱量,從所考察內部介質來看它應由 F o u r i e r Fourier Fourier定律確定,即:
d Q = − λ ∂ u ∂ n d S d t (2) d Q=-\lambda \frac{\partial u}{\partial n} d S d t \tag{2} dQ=λnudSdt(2)
其中 ∂ u ∂ n \frac{\partial u}{\partial n} nu錶示 u u u沿邊界 S S S上的單比特外法線方向 n n n的方向導數。從外部方面來看則應由牛頓冷卻定律决定,即:
d Q = h ( u − U ) d S d t (3) d Q=h\left(u-U\right) d S d t \tag{3} dQ=h(uU)dSdt(3)
結合 ( 2 ) ( 3 ) (2)(3) (2)(3)得到第三類邊界條件:
− λ ∂ u ∂ n = h ( u − U ) (4) -\lambda \frac{\partial u}{\partial n}=h\left(u-U\right) \tag{4} λnu=h(uU)(4)


2網格剖分

2.1 對符號更細致的說明
如下圖所示,以焊接區域中心的上側與爐內空氣接觸處為原點,指向電路板內部為正方向建立 x x x軸,熱量沿 x x x軸方向傳遞。
在這裏插入圖片描述
由於接觸面環境溫度 U U U是與時間 t t t和物件速度 v v v有關,則實際接觸面環境溫度寫作 U ( v , t ) U(v,t) U(v,t)較為合適,其中 v t vt vt為物件橫向移動距離:

在這裏插入圖片描述

因此我們可以將第一部分熱傳導方程進行如下整理:
.
2.2 方程整理
內部(熱傳導):
∂ u ( x , t ) ∂ t = a ∂ 2 u ( x , t ) ∂ x 2 \frac{\partial u(x, t)}{\partial t}=a \frac{\partial^{2} u(x, t)}{\partial x^{2}} tu(x,t)=ax22u(x,t)
上下兩邊界(第三邊界條件):
− λ ∂ u ( x , t ) ∂ t ∣ x = 0 + h u ( x , t ) ∣ x = 0 = h U ( v , t ) λ ∂ u ( x , t ) ∂ t ∣ x = d + h u ( x , t ) ∣ x = d = h U ( v , t ) \begin{aligned} &-\left.\lambda \frac{\partial u(x, t)}{\partial t}\right|_{x=0}+\left.h u(x, t)\right|_{x=0}=h U(v, t) \\ &\left.\quad\lambda \frac{\partial u(x, t)}{\partial t}\right|_{x=d}+\left.h u(x, t)\right|_{x=d}=h U(v, t) \end{aligned} λtu(x,t)x=0+hu(x,t)x=0=hU(v,t)λtu(x,t)x=d+hu(x,t)x=d=hU(v,t)
初值條件:
t = 0 t=0 t=0時,我們認為電路板溫度與生產車間的溫度 T 0 T_0 T0保持一致,故初值條件為:
u ( x , 0 ) = T 0 u(x,0)=T_0 u(x,0)=T0
整理:
{ ∂ u ( x , t ) ∂ t = a ∂ 2 u ( x , t ) ∂ x 2 − λ ∂ u ( x , t ) ∂ t ∣ x = 0 + h u ( x , t ) ∣ x = 0 = h U ( v , t ) λ ∂ u ( x , t ) ∂ t ∣ x = d + h u ( x , t ) ∣ x = d = h U ( v , t ) u ( x , 0 ) = T 0 \left\{\begin{array}{c} \frac{\partial u(x, t)}{\partial t}=a \frac{\partial^{2} u(x, t)}{\partial x^{2}} \\ -\left.\lambda \frac{\partial u(x, t)}{\partial t}\right|_{x=0}+\left.h u(x, t)\right|_{x=0}=h U(v, t) \\ \begin{array}{c} \left.\quad\lambda \frac{\partial u(x, t)}{\partial t}\right|_{x=d}+\left.h u(x, t)\right|_{x=d}=h U(v, t) \\ u(x, 0)=T_{0} \end{array} \end{array}\right. tu(x,t)=ax22u(x,t)λtu(x,t)x=0+hu(x,t)x=0=hU(v,t)λtu(x,t)x=d+hu(x,t)x=d=hU(v,t)u(x,0)=T0
.
2.3 網格拆分
我們對於方向 x x x及方向 t t t進行網格拆分,為敘述簡便起見我們記 u k , j = u ( x k , t j ) u_{k,j}=u(x_k,t_j) uk,j=u(xk,tj),其中: k = 0 , 1 , … , n , j = 0 , 1 , … , m , n = [ d Δ x ] , m = ∣ L v Δ t ⌋ \left.k=0,1, \ldots, n, j=0,1, \ldots, m, \quad n=\left[\frac{d}{\Delta x}\right], m=\mid \frac{L}{v \Delta t}\right\rfloor k=0,1,,n,j=0,1,,m,n=[Δxd],m=vΔtL
在這裏插入圖片描述
初始條件:
u k , 0 = u ( x k , 0 ) = T 0 ( k = 0 , 1 , … , n ) u_{k, 0}=u\left(x_{k}, 0\right)=T_{0} \quad(k=0,1, \ldots, n) uk,0=u(xk,0)=T0(k=0,1,,n)

內部(熱傳導):

∂ u ∂ t \frac{\partial u}{\partial t} tu采用向後差分公式:

∂ u ∂ t ∣ ( k , j ) = u k , j − u k , j − 1 Δ t + O ( Δ t ) \left.\frac{\partial u}{\partial t}\right|_{(k, j)}=\frac{u_{k, j}-u_{k, j-1}}{\Delta t}+O(\Delta t) tu(k,j)=Δtuk,juk,j1+O(Δt)

∂ 2 u ∂ x 2 \frac{\partial^{2} u}{\partial x^{2}} x22u采用二階中心差商公式:
∂ 2 u ∂ x 2 ∣ ( k , j ) = u k + 1 , j − 2 u k , j + u k − 1 , j Δ x 2 + O ( Δ x 2 ) \left.\frac{\partial^{2} u}{\partial x^{2}}\right|_{(k, j)}=\frac{u_{k+1, j}-2 u_{k, j}+u_{k-1, j}}{\Delta x^{2}}+O\left(\Delta x^{2}\right) x22u(k,j)=Δx2uk+1,j2uk,j+uk1,j+O(Δx2)
則上述一維熱傳導方程式可錶示為:

u k , j − u k , j − 1 Δ t − a u k + 1 , j − 2 u k , j + u k − 1 , j Δ x 2 = O ( Δ t + Δ x 2 ) \frac{u_{k, j}-u_{k, j-1}}{\Delta t}-a \frac{u_{k+1, j}-2 u_{k, j}+u_{k-1, j}}{\Delta x^{2}}=O\left(\Delta t+\Delta x^{2}\right) Δtuk,juk,j1aΔx2uk+1,j2uk,j+uk1,j=O(Δt+Δx2)
近似為:

u k , j − u k , j − 1 Δ t − a u k + 1 , j − 2 u k , j + u k − 1 , j Δ x 2 = 0 \frac{u_{k, j}-u_{k, j-1}}{\Delta t}-a \frac{u_{k+1, j}-2 u_{k, j}+u_{k-1, j}}{\Delta x^{2}}=0 Δtuk,juk,j1aΔx2uk+1,j2uk,j+uk1,j=0
上下兩邊界(第三邊界條件):
相似的我們可以獲得邊界處溫度變化方程:
{ − u 1 , j − u 0 , j Δ x + γ u 0 , j = γ U ( v , t j ) u n , j − u n − 1 , j Δ x + γ u n , j = γ U ( v , t j ) \left\{\begin{array}{l} -\frac{u_{1, j}-u_{0, j}}{\Delta x}+\gamma u_{0, j}=\gamma U\left(v, t_{j}\right) \\ \frac{u_{n, j}-u_{n-1, j}}{\Delta x}+\gamma u_{n, j}=\gamma U\left(v, t_{j}\right) \end{array}\right. { Δxu1,ju0,j+γu0,j=γU(v,tj)Δxun,jun1,j+γun,j=γU(v,tj)
其中 γ = h λ \gamma=\frac{h}{\lambda} γ=λh


3三對角矩陣

依據上述差分近似方程,我們可以列出形式如下的三對角遞推線性非齊次方程組:
A u j + 1 → = u j → + f j → A = ( 1 + 2 F 0 − F o 1 + B i − F o 0 ⋯ 0 0 − F 0 1 + 2 F o − F o ⋯ 0 0 ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 0 ⋯ 1 + 2 F o − F o 0 0 0 ⋯ − F o 1 + 2 F o − F o 1 + B i ) , u j ‾ = ( u 1 , j u 2 , j ⋮ ⋮ u n − 2 , j u n − 1 , j ) , f j → = ( U ( v , t j ) 0 ⋮ ⋮ 0 U ( v , t j ) ) , ( j = 0 , … , m ) A \overrightarrow{u_{j+1}}=\overrightarrow{u_{j}}+\overrightarrow{f_{j}} \\ A=\left(\begin{array}{cccccc} 1+2 F_{0}-\frac{F_{o}}{1+B i} & -F_{o} & 0 & \cdots & 0 & 0 \\ -F_{0} & 1+2 F_{o} & -F_{o} & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1+2 F_{o} & -F_{o} \\ 0 & 0 & 0 & \cdots & -F_{o} & 1+2 F_{o}-\frac{F_{o}}{1+B i} \end{array}\right), \overline{u_{j}}=\left(\begin{array}{c} u_{1, j} \\ u_{2, j} \\ \vdots \\ \vdots \\ u_{n-2, j} \\ u_{n-1, j} \end{array}\right), \overrightarrow{f_{j}}=\left(\begin{array}{c} U\left(v, t_{j}\right) \\ 0 \\ \vdots \\ \vdots \\ 0 \\ U\left(v, t_{j}\right) \end{array}\right), \\ (j=0, \ldots, m) Auj+1 =uj +fj A=1+2F01+BiFoF000Fo1+2Fo000Fo00001+2FoFo00Fo1+2Fo1+BiFo,uj=u1,ju2,jun2,jun1,j,fj =U(v,tj)00U(v,tj),(j=0,,m)
其中 F o = a Δ t Δ x 2 , B i = γ Δ x = h λ Δ x F_{o}=\frac{a \Delta t}{\Delta x^{2}}, \quad B i=\gamma \Delta x=\frac{h}{\lambda} \Delta x Fo=Δx2aΔt,Bi=γΔx=λhΔx分別為傳熱學中的網格傅裏葉數和網格畢奧數。


4MATLAB模擬

4.1 模擬問題再描述

某回焊爐內有11個小溫區及爐前區域和爐後區域,每個小溫區長度為30.5 cm,相鄰小溫區之間有5 cm的間隙,爐前區域和爐後區域長度均為25 cm:
在這裏插入圖片描述
各溫區設定的溫度分別為175ºC(小溫區1 ~ 5)、195ºC(小溫區6)、235ºC(小溫區7)、255ºC(小溫區8 ~ 9)及25ºC(小溫區10 ~ 11);傳送帶的過爐速度為70 cm/min;焊接區域的厚度為0.15 mm。溫度傳感器在焊接區域中心的溫度達到30ºC時開始工作,電路板進入回焊爐開始計時。

假設 a = 4.41 × 1 0 − 5   m 2 / s , γ = 3.53 × 1 0 − 2   m − 1 a=4.41 \times 10^{-5} \mathrm{~m}^{2} / \mathrm{s}, \gamma=3.53 \times 10^{-2} \mathrm{~m}^{-1} a=4.41×105 m2/s,γ=3.53×102 m1
F o = 196000 , B i = 5.3 e − 08 F_o=196000,Bi=5.3e-08 Fo=196000,Bi=5.3e08

以下使用MATLAB模擬在該條件下焊接元件中心區域溫度變化:

4.2 相關代碼

function reflowProfile
% @author : slandarer

% 參數定義及計算 ==========================================================
% 溫區相關數據 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
warmZone.Len=30.5;      % 溫區長度(cm)           
warmZone.SepLen=5;      % 溫區間隙長度(cm)           
warmZone.ForeLen=25;    % 爐前區域長度(cm)       
warmZone.BackLen=25;    % 爐後區域長度(cm)       
warmZone.Num=11;        % 溫區數量     
% 溫區總長=溫區長度*溫區數量+間隙長度*(溫區數量-1)+爐前長度+爐後長度
% warmZone.TotalLen=30.5*11+5*10+25+25;
warmZone.TotalLen=warmZone.Len*warmZone.Num+...
                  warmZone.SepLen*(warmZone.Num-1)+...
                  warmZone.ForeLen+...
                  warmZone.BackLen;

% 每個大溫區包含哪幾個小溫區
warmZone.Zone{
    1}=[1 2 3 4 5];               
warmZone.Zone{
    2}=6;                         
warmZone.Zone{
    3}=7;
warmZone.Zone{
    4}=[8 9];
warmZone.Zone{
    5}=[10 11];      

% 設置每個溫區溫度
warmZone.Temp(warmZone.Zone{
    1})=175;
warmZone.Temp(warmZone.Zone{
    2})=195;
warmZone.Temp(warmZone.Zone{
    3})=235;
warmZone.Temp(warmZone.Zone{
    4})=255;
warmZone.Temp(warmZone.Zone{
    5})=25;


% 電路板相關數據 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
ccBoard.v_cm_min=70;     % 電路板移動速度(cm/min)
ccBoard.v_cm_s=70/60;    % 電路板移動速度(cm/s)  
ccBoard.d=0.15;          % 焊接區域厚度(mm)      
ccBoard.Temp0=25;        % 電路板初始溫度(C)

% 以下屬性在該篇博文中並未用到
% ccBoard.Lim.ChangeRate=[-3 3];  % 溫度變化率上下限
% ccBoard.Lim.RiseTime=[60 120];  % 溫度上昇過程中在150ºC~190ºC的時間限制
% ccBoard.Lim.PeakTime=[40 90];   % 溫度大於217ºC的時間上下限
% ccBoard.Lim.PeakTemp=[240 250]; % 峰值溫度上下限


% 其他相關參數計算- - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
totalTime=warmZone.TotalLen./ccBoard.v_cm_s;
disp(['焊接區域比特於回焊爐內部時長:',num2str(totalTime)]) 

% 獲取各個溫區拐點和中點比特置(用於插值外界溫度曲線)
wzPosList(3:3:3*warmZone.Num)=(1:warmZone.Num).*(warmZone.Len+warmZone.SepLen);
wzPosList(2:3:3*warmZone.Num)=wzPosList(3:3:3*warmZone.Num)-warmZone.SepLen;
wzPosList(1:3:3*warmZone.Num)=wzPosList(3:3:3*warmZone.Num)-warmZone.SepLen-warmZone.Len/2;
wzPosList(end)=[];

% 獲取各個溫區拐點和中點比特置溫度(用於插值外界溫度曲線)
wzTempList(3:3:3*warmZone.Num)=warmZone.Temp;
wzTempList(2:3:3*warmZone.Num)=warmZone.Temp;
wzTempList(1:3:3*warmZone.Num)=warmZone.Temp;

% 這裏用end-6是因為依據題目所給圖像,最後10~11溫區並不是直接到25度,也需要插值
posNodes=[0 warmZone.ForeLen warmZone.ForeLen+wzPosList(1:end-6),...
            warmZone.ForeLen+warmZone.BackLen+wzPosList(end)]; % 比特置節點
% timeNodes=posNodes./ccBoard.v_cm_s;                            % 時間節點
tempNodes=[ccBoard.Temp0 wzTempList(1:end-6) ccBoard.Temp0];   % 溫度節點

% 用於進行溫度插值
% interp1(posNodes,tempNodes,pos);
timeSet=0:0.01:totalTime;                 % 將時間進行細分
posSet=timeSet.*ccBoard.v_cm_s;           % 元件中心比特置
U=interp1(posNodes,tempNodes,posSet);     % 元件中心比特置接觸面環境溫度


% 三對角矩陣構建 ==========================================================
N=101;  % 將元件細分的取的樣點數,取奇數是希望中間點恰巧被取到

u=25.*ones(N,1);  % 元件溫度分布,初始每一處都是25A=zeros(N,N);     % 初始化三對角矩陣

Fo=196000;        % 網格傅裏葉數
Bi=5.3e-08;       % 網格畢奧數

% 三對角矩陣賦值
A(diag(1:N)~=0)=1+2*Fo;
A(diag(1:N-1,1)~=0)=-Fo;
A(diag(1:N-1,-1)~=0)=-Fo;
A(1,1)=A(1,1)-Fo/(1+Bi);
A(end,end)=A(end,end)-Fo/(1+Bi);

invA=eye(N)/A;    % 三對角矩陣的逆矩陣



% 數據計算 ================================================================
for i=1:length(timeSet)
    f=zeros(N,1); % 由外界溫度决定的附加項
    f([1,N])=U(i)*Fo*Bi/(1+Bi);
    u(:,i+1)=invA*u(:,i)+invA*f;
end

% 獲取中間處溫度,這裏向上向下取整是應對N取偶數的情况
mid_u=(u(floor((N+1)/2),:)+u(ceil((N+1)/2),:))./2;



% 繪圖 ====================================================================
% 繪制爐溫曲線 
plot(timeSet,mid_u(1:end-1),'LineWidth',1.5)

% axes屬性設置
ax=gca;
hold(ax,'on');
box on
grid on
ax.LineWidth = 1;
ax.XLim=[0,373];
ax.GridLineStyle='--';
% X軸標簽
ax.XLabel.String='t(s)';
ax.XLabel.FontSize=13;
ax.XLabel.FontName='Cambria';
% Y軸標簽
ax.YLabel.String='T(^{\circ}C)';
ax.YLabel.FontSize=13;
ax.YLabel.FontName='Cambria';

% 繪制217ºC溫度線
plot(timeSet([1,end]),[217 217],'LineWidth',1.5,...
    'Color',[.6,.6,.6],'LineStyle','--')


end

4.3 模擬結果
在這裏插入圖片描述
在這裏插入圖片描述


5後言

本篇文章雖然只講解了如何基於三對角矩陣求解熱傳導方程,但實際上國賽題目2020A所有問題基本上都是在學會會了該方法的基礎上,在一定的限制條件下對部分參數進行更改和搜索以找出最優參數組,在此不做詳述。

版權聲明
本文為[slandarer]所創,轉載請帶上原文鏈接,感謝
https://cht.chowdera.com/2021/08/20210820063515164g.html

隨機推薦