国赛助力：第三类边界条件热传导方程及基于三对角矩阵的数值计算MATLAB实现（2020A）

1第三类边界条件的热传导方程

1.1 热传导方程
热传导在一维的各向同性介质里的传播可用以下方程表达：
∂u∂t=a∂2u∂x2(1)\frac{\partial u}{\partial t}=a \frac{\partial^{2} u}{\partial x^{2}} \tag{1} ∂t∂u=a∂x2∂2u(1)

其中，u=u(x,t)u=u(x,t)u=u(x,t)，a=λcρa=\frac{\lambda}{c\rho}a=cρλ， λ\lambdaλ表示介质的热传导率， ccc表示介质的比热， ρ\rhoρ表示介质的密度。
.
1.2 第三类边界条件
考察介质放在另一种介质中的情形。外界介质的温度UUU与所考察介质表面上的温度uuu往往并不相同，考虑流过所考察介质表面的热量，从所考察内部介质来看它应由FourierFourierFourier定律确定，即：
dQ=−λ∂u∂ndSdt(2)d Q=-\lambda \frac{\partial u}{\partial n} d S d t \tag{2} dQ=−λ∂n∂udSdt(2)
其中∂u∂n\frac{\partial u}{\partial n}∂n∂u表示uuu沿边界SSS上的单位外法线方向nnn的方向导数。从外部方面来看则应由牛顿冷却定律决定，即：
dQ=h(u−U)dSdt(3)d Q=h\left(u-U\right) d S d t \tag{3} dQ=h(u−U)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} −λ∂n∂u=h(u−U)(4)

2网格剖分

2.1 对符号更细致的说明
如下图所示，以焊接区域中心的上侧与炉内空气接触处为原点，指向电路板内部为正方向建立xxx轴，热量沿xxx轴方向传递。

由于接触面环境温度UUU是与时间ttt和物件速度vvv有关，则实际接触面环境温度写作U(v,t)U(v,t)U(v,t)较为合适，其中vtvtvt为物件横向移动距离：

因此我们可以将第一部分热传导方程进行如下整理：
.
2.2 方程整理
内部（热传导）：
∂u(x,t)∂t=a∂2u(x,t)∂x2\frac{\partial u(x, t)}{\partial t}=a \frac{\partial^{2} u(x, t)}{\partial x^{2}} ∂t∂u(x,t)=a∂x2∂2u(x,t)
上下两边界（第三边界条件）：
−λ∂u(x,t)∂t∣x=0+hu(x,t)∣x=0=hU(v,t)λ∂u(x,t)∂t∣x=d+hu(x,t)∣x=d=hU(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} −λ∂t∂u(x,t)∣∣∣∣x=0+hu(x,t)∣x=0=hU(v,t)λ∂t∂u(x,t)∣∣∣∣x=d+hu(x,t)∣x=d=hU(v,t)
初值条件：
在t=0t=0t=0时，我们认为电路板温度与生产车间的温度T0T_0T0保持一致，故初值条件为：
u(x,0)=T0u(x,0)=T_0 u(x,0)=T0
整理：
{∂u(x,t)∂t=a∂2u(x,t)∂x2−λ∂u(x,t)∂t∣x=0+hu(x,t)∣x=0=hU(v,t)λ∂u(x,t)∂t∣x=d+hu(x,t)∣x=d=hU(v,t)u(x,0)=T0\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. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧∂t∂u(x,t)=a∂x2∂2u(x,t)−λ∂t∂u(x,t)∣∣∣x=0+hu(x,t)∣x=0=hU(v,t)λ∂t∂u(x,t)∣∣∣x=d+hu(x,t)∣x=d=hU(v,t)u(x,0)=T0
.
2.3 网格拆分
我们对于方向xxx及方向ttt进行网格拆分，为叙述简便起见我们记 uk,j=u(xk,tj)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=∣LvΔ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\rfloork=0,1,…,n,j=0,1,…,m,n=[Δxd],m=∣vΔtL⌋：

初始条件：
uk,0=u(xk,0)=T0(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}∂t∂u采用向后差分公式：

∂u∂t∣(k,j)=uk,j−uk,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) ∂t∂u∣∣∣∣(k,j)=Δtuk,j−uk,j−1+O(Δt)

对∂2u∂x2\frac{\partial^{2} u}{\partial x^{2}}∂x2∂2u采用二阶中心差商公式：
∂2u∂x2∣(k,j)=uk+1,j−2uk,j+uk−1,jΔx2+O(Δx2)\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) ∂x2∂2u∣∣∣∣(k,j)=Δx2uk+1,j−2uk,j+uk−1,j+O(Δx2)
则上述一维热传导方程式可表示为：

uk,j−uk,j−1Δt−auk+1,j−2uk,j+uk−1,jΔx2=O(Δt+Δx2)\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,j−uk,j−1−aΔx2uk+1,j−2uk,j+uk−1,j=O(Δt+Δx2)
近似为：

uk,j−uk,j−1Δt−auk+1,j−2uk,j+uk−1,jΔx2=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,j−uk,j−1−aΔx2uk+1,j−2uk,j+uk−1,j=0
上下两边界（第三边界条件）：
相似的我们可以获得边界处温度变化方程：
{−u1,j−u0,jΔx+γu0,j=γU(v,tj)un,j−un−1,jΔx+γun,j=γU(v,tj)\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,j−u0,j+γu0,j=γU(v,tj)Δxun,j−un−1,j+γun,j=γU(v,tj)
其中γ=hλ\gamma=\frac{h}{\lambda}γ=λh

3三对角矩阵

依据上述差分近似方程，我们可以列出形式如下的三对角递推线性非齐次方程组：
Auj+1→=uj→+fj→A=(1+2F0−Fo1+Bi−Fo0⋯00−F01+2Fo−Fo⋯00⋮⋮⋱⋮⋮000⋯1+2Fo−Fo000⋯−Fo1+2Fo−Fo1+Bi),uj‾=(u1,ju2,j⋮⋮un−2,jun−1,j),fj→=(U(v,tj)0⋮⋮0U(v,tj)),(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+fjA=⎝⎜⎜⎜⎜⎜⎛1+2F0−1+BiFo−F0⋮00−Fo1+2Fo⋮000−Fo⋱00⋯⋯⋯⋯00⋮1+2Fo−Fo00⋮−Fo1+2Fo−1+BiFo⎠⎟⎟⎟⎟⎟⎞,uj=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛u1,ju2,j⋮⋮un−2,jun−1,j⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞,fj=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛U(v,tj)0⋮⋮0U(v,tj)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞,(j=0,…,m)
其中Fo=aΔtΔx2,Bi=γΔx=hλΔxF_{o}=\frac{a \Delta t}{\Delta x^{2}}, \quad B i=\gamma \Delta x=\frac{h}{\lambda} \Delta xFo=Δ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×10−5m2/s,γ=3.53×10−2m−1a=4.41 \times 10^{-5} \mathrm{~m}^{2} / \mathrm{s}, \gamma=3.53 \times 10^{-2} \mathrm{~m}^{-1}a=4.41×10−5 m2/s,γ=3.53×10−2 m−1
即Fo=196000,Bi=5.3e−08F_o=196000,Bi=5.3e-08Fo=196000,Bi=5.3e−08

以下使用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);  % 元件温度分布，初始每一处都是25度
A=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所有问题基本上都是在学会会了该方法的基础上，在一定的限制条件下对部分参数进行更改和搜索以找出最优参数组，在此不做详述。