大家好，又见面了，我是你们的朋友全栈君。
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} ∂t∂u=a∂x2∂2u(1)
其中， u = u ( x , t ) u=u(x,t) u=u(x,t)， a = λ c ρ a=frac{lambda}{crho} 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=−λ∂n∂udSdt(2) 其中 ∂ u ∂ n frac{partial u}{partial n} ∂n∂u表示 u u u沿边界 S S S上的单位外法线方向 n n n的方向导数。从外部方面来看则应由牛顿冷却定律决定，即： d Q = h ( u − U ) d S d t (3) d Q=hleft(u-Uright) 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}=hleft(u-Uright) tag{4} −λ∂n∂u=h(u−U)(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}} ∂t∂u(x,t)=a∂x2∂2u(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.quadlambda 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 = 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.quadlambda 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 网格拆分 我们对于方向 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}rightrfloor 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}=uleft(x_{k}, 0right)=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 ) = 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) ∂t∂u∣∣∣∣(k,j)=Δtuk,j−uk,j−1 O(Δt)
对 ∂ 2 u ∂ x 2 frac{partial^{2} u}{partial x^{2}} ∂x2∂2u采用二阶中心差商公式： ∂ 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}} Oleft(Delta x^{2}right) ∂x2∂2u∣∣∣∣(k,j)=Δx2uk 1,j−2uk,j uk−1,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}}=Oleft(Delta t Delta x^{2}right) Δtuk,j−uk,j−1−aΔx2uk 1,j−2uk,j uk−1,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,j−uk,j−1−aΔx2uk 1,j−2uk,j uk−1,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 Uleft(v, t_{j}right) \ frac{u_{n, j}-u_{n-1, j}}{Delta x} gamma u_{n, j}=gamma Uleft(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三对角矩阵
依据上述差分近似方程，我们可以列出形式如下的三对角递推线性非齐次方程组： 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} Uleft(v, t_{j}right) \ 0 \ vdots \ vdots \ 0 \ Uleft(v, t_{j}right) end{array}right), \ (j=0, ldots, m) Auj 1 =uj fj A=⎝⎜⎜⎜⎜⎜⎛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) 其中 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×10−5 m2/s,γ=3.53×10−2 m−1 即 F o = 196000 , B i = 5.3 e − 08 F_o=196000,Bi=5.3e-08 Fo=196000,Bi=5.3e−08
以下使用MATLAB模拟在该条件下焊接元件中心区域温度变化：
4.2 相关代码
代码语言：javascript复制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','--')
endJetbrains全家桶1年46，售后保障稳定
4.3 模拟结果
5后言
本篇文章虽然只讲解了如何基于三对角矩阵求解热传导方程，但实际上国赛题目2020A所有问题基本上都是在学会会了该方法的基础上，在一定的限制条件下对部分参数进行更改和搜索以找出最优参数组，在此不做详述。
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。
发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/209924.html原文链接：https://javaforall.cn
