基于微分方程的回焊爐焊接區(qū)域溫度分析

胡思哲，呂東育，姜高龍

（南京郵電大學(xué) 通信與信息工程學(xué)院， 南京 210023）

1 問題重述

2020 年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽（CUMCM）A 賽題的回焊爐設(shè)有11 個小溫區(qū)、爐前區(qū)域和爐后區(qū)域，每個小溫區(qū)長度為30.5 cm，相鄰小溫區(qū)之間有5 cm 的間隙，爐前區(qū)域和爐后區(qū)域長度均為25 cm[1]。附件給出某次實驗中爐溫曲線的數(shù)據(jù)，各溫區(qū)設(shè)定的溫度分別為175 ℃（小溫區(qū) 1—5）、195 ℃（小溫區(qū) 6）、235 ℃（小溫區(qū)7）、255 ℃（小溫區(qū) 8～9） 及 25 ℃（小溫區(qū) 10—11）。生產(chǎn)車間環(huán)境溫度保持25 ℃，傳送帶的過爐速度為70 cm/min，焊接區(qū)域的厚度為0.15 mm。溫度傳感器在焊接區(qū)域中心的溫度達(dá)到30 ℃時開始工作，電路板進(jìn)入回焊爐開始計時。

要求通過機(jī)理分析，建立數(shù)學(xué)模型，研究分析相關(guān)條件下的焊接區(qū)域中心溫度變化情況，列出小溫區(qū)域中心溫度，畫出相應(yīng)的爐溫曲線；確定允許的傳送帶最大過爐速度和最優(yōu)爐溫曲線；確定各溫區(qū)的設(shè)定溫度和傳送帶的過爐速度。

2 基于微分方程的熱傳導(dǎo)模型

2.1 確定熱傳導(dǎo)方程

熱力學(xué)定律——傅里葉定律指出[2]：單位時間通過給定截面的熱量與該截面垂直方向上的溫度變化率和橫截面積成正比，并且傳遞方向與溫度升高方向正好相反。假設(shè)回焊爐內(nèi)的傳熱方式以一維熱傳導(dǎo)為主，結(jié)合傅里葉熱傳導(dǎo)定律[3]和能量守恒定律，可得到熱傳導(dǎo)方程：

其中：u=u（t，x）為 t 時刻 x 坐標(biāo)處的溫度；c 為比熱容，ρ 為密度，k 為熱傳導(dǎo)率。

2.2 確定回焊爐溫度分布

根據(jù)回焊爐內(nèi)各區(qū)域的長度結(jié)構(gòu)和恒溫情況，可確定各恒溫區(qū)橫向端點的坐標(biāo)d0～d11依次為：0，15，25，197.5，202.5，233，238，268.5，273.5，339.5，344.5，435.5，單位為 cm。dr處的溫度記為ur（r=0，1，2，…，11）。已知恒溫情況為（℃）：

由穩(wěn)態(tài)導(dǎo)熱性質(zhì)[4]可知，在穩(wěn)態(tài)導(dǎo)熱過程中，溫度分布與時間無關(guān)，方程（1）成為：

由此解得回焊爐中溫度分布為分段線性函數(shù)：

3 焊接區(qū)域中心的溫度變化模型

3.1 偏微分方程模型

焊接區(qū)域中心的熱傳導(dǎo)方程仍為（1），其中x∈[l0，l1]是焊接區(qū)域中心的縱坐標(biāo)，l0、l1代表電路板的上下邊界，l1-l0=0.15 mm。

系統(tǒng)的初始條件為

設(shè)ke為邊界與高溫環(huán)境間的熱交換系數(shù)，結(jié)合牛頓冷卻定律和傅里葉熱傳導(dǎo)定律，可得系統(tǒng)的邊界條件為

其中ue是環(huán)境溫度，可由函數(shù)（2）給出ue=ue（t）=U（vt），v 為已知的傳送帶過爐速度。

綜上，式（1）—（5）構(gòu)成了焊接區(qū)域中心熱傳導(dǎo)的偏微分方程模型。

3.2 偏微分方程的差分格式

偏微分方程的數(shù)值求解方法通常有有限元素法[5]與有限差分法[6]兩種。上述焊接區(qū)域中心熱傳導(dǎo)模型函數(shù)均為一維偏微分方程，應(yīng)用有限元素法優(yōu)勢并不明顯，故采用有限差分法。

溫度函數(shù) u=u（t，x）定義在矩形區(qū)域[0，T]×[l0，l1]內(nèi)，其中 T 滿足 vT=d11，傳送帶過爐速度 v和回焊爐長度d11均已知。以時間步長Δt 和空間步長Δx 將定義區(qū)域分割成矩形網(wǎng)格，格點處的溫度記為表示第n 個時間點、第i 個空間位置的溫度。分別用一階差商和二階中心差商代替偏導(dǎo)數(shù)：

其中M 代表電路板最后層格點，滿足MΔx=l1-l0。

將這些差商代入式（1）、（4）、（5），再令

可得偏微分方程的差分格式

4 差分格式的數(shù)值求解

4.1 追趕法的求解原理

追趕法的基本原理是：將矩陣A 分解為A=LU，L 是對角線上元素全為1 的下二對角矩陣，U為上二對角矩陣；于是方程組Ax=B 化為兩個易解的方程組Ly=B 和Ux=y；先求得y=L-1B，再求解Ux=L-1B，即可得到各位置、全時間下的溫度分布。

4.2 基于最小二乘法的最佳參數(shù)確定

為了確定差分格式（6）中的參數(shù)λ 和μ，先取Δt=0.1 s，Δx=0.001 mm，然后建立最小二乘法的目標(biāo)規(guī)劃模型：

其中：u（t）為求解模型（6）獲得的 t 時刻焊接區(qū)域中心溫度，可稱之為理論溫度；u0（t）為附件[1]給出的實測溫度。

查閱資料表明[8]，在電路板到達(dá)冷卻區(qū)后，由于表面風(fēng)速變大，導(dǎo)致熱交換系數(shù)發(fā)生變化，因而參數(shù)μ 會發(fā)生變化，所以取

結(jié)合電路板的材質(zhì)等方面資料[9]，對λ、μ1、μ2的可能范圍分別進(jìn)行估計，得到：

具體做法是按一定的步長在上述范圍內(nèi)取定λ、μ1、μ2的值代入模型（6），用追趕法求出溫度函數(shù) u（t），以此搜索目標(biāo)規(guī)劃（7）的最小值。

為了提高搜索效率，按以下步驟實施搜索。

步驟1：在時間軸上對于0≤vt ＜d9（電路板位于水平方向冷卻區(qū)之前）的部分，搜索λ∈Λ 和μ1∈Ω1。

先取步長 Δλ =0.05× 106，Δμ1=0.1× 106，遍歷搜索得到 λ =0.5×106，μ1=4.1×106。

再在 λ∈[0.45 × 106，0.55 × 106]、μ1∈[4.05 ×106，4.15× 106]范圍內(nèi)，取步長 Δλ =0.01×106，Δμ1=0.01× 106，遍歷搜索得到 λ =0.5× 106，μ1=4.08×106。

步驟2：在時間軸上對于d9≤vt≤d11（電路板位于水平方向冷卻區(qū)之后）的部分，以及λ=0.5×106的情況下，搜索 μ2∈Ω2。

先取步長Δμ2=0.1×106，遍歷搜索得到μ2=12.4×106。

再取步長 Δμ2= 0.01 × 106，在 μ2∈[12.3 ×106，12.5 × 106]內(nèi)遍歷搜索得到 μ2=12.37 × 106。

4.3 模型優(yōu)化

第4.2 節(jié)得到參數(shù)的最佳值為λ=0.5×106，μ1= 4.08 × 106，μ2= 12.37 × 106。代入差分格式（6），求得焊接區(qū)域中心的溫度，在Matlab 軟件中繪制得到實測與理論爐溫曲線如圖1 所示。

圖1 爐溫理論曲線與實測曲線對比

觀察圖1 發(fā)現(xiàn)，爐溫理論曲線與實測曲線有三處分離程度過大，誤差不在可接受范圍內(nèi)，所以對模型進(jìn)行優(yōu)化。①②③三處誤差產(chǎn)生的原因分析如下：

（1）回顧回焊爐溫度分布函數(shù)（2），在爐前區(qū)域，從d=d1起溫度開始上升，設(shè)置d1=15 較小，溫度提升過早。

（2）查閱資料發(fā)現(xiàn)[10]，當(dāng)溫度高于200 ℃時，焊料開始熔化，其物理性質(zhì)發(fā)生較大變化，會影響參數(shù)λ 的值，所以需對λ 進(jìn)行分段處理。

（3）查閱資料表明[11]，在實際生產(chǎn)中，回焊爐中回流區(qū)與冷卻區(qū)間的過渡區(qū)較大（即不止一個間隙的寬度），所以要加大d10的取值。

根據(jù)上述分析，采取兩項措施予以改進(jìn)：一是在函數(shù)（2）中，將原先的 d1=15 和 d10=344.5，修改為d1=20，d10=380；二是 u≤200 ℃時，λ 取 λ1，反之λ 取λ2，分別進(jìn)行搜索確定。

進(jìn)行上述修改后，重復(fù)第4.1 和4.2 的過程，求解模型（6）和（7），可搜索得到參數(shù)的取值為λ1= 0.55 × 106，λ2= 0.71 × 106，μ1= 4.08 × 106，μ2=10.27 × 106。

將上述參數(shù)值代入，求解差分格式（6），可得到焊接區(qū)域的溫度分布。利用Matlab 軟件繪制溫度分布曲線如圖2 所示，繪制理論曲線與實測曲線并對比如圖3 所示。

圖2 焊接區(qū)域溫度分布曲線

圖3 模型優(yōu)化前后的爐溫曲線對比

由圖3 可見，兩條曲線基本重合，分離程度嚴(yán)重的三處消失，充分驗證了模型的準(zhǔn)確性。