# Effects of the Solder Phase Transformation on the Optimization of Reflow Soldering Parameters and Temperature Profiles

Authors: Yongze Xu; Yu Liu; Lei Zhang
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this study, a heat convection model of the reflow oven and a heat conduction model of the soldering area are proposed based on heat transfer theory, and a dynamic Thomas algorithm is developed for solving linear equations with coefficient matrix evolving over time in the tridiagonal system, which is derived from a heat transfer problem with moving boundaries in the solder phase transition process. We have also carried out numerical simulations for investigating the accuracy of the mathematical model, in which the temperature profiles are calculated and compared for different cases with considering or ignoring phase transformations, respectively. Parameters of reflow soldering, such as the conveyor speed, the set temperature in each zone, and a part of the heating factor, are optimized by the use of the nondominated sorting genetic algorithm II. By comparing the temperature profile and optimal solutions in the two cases, numerical results show that phase transitions of the solder have great impacts on optimal parameters and the slope of temperature profiles. Moreover, the phenomenon that the heating factor varies with the maximum set temperature in a banded distribution is investigated and analyzed, which is an important part of this work.


---

## Body

## 1. Introduction

In electronic manufacturing industries such as integrated circuit boards, the soldering process of printed circuit boards (PCB) is of particular importance. Currently, the PCB soldering equipment and technology are diverse, represented by reflow ovens and surface mount technology. The main factor of reflow soldering is the temperature profile, that is, a curve of the temperature in the center of the soldering area [1]. Moreover, the set temperature of each zone and the conveyor speed play a crucial role in the quality of products. Therefore, how to adjust and control them reasonably, especially the optimization of the temperature profile, has been a research focus [2–6]. At present, most work in this field is checked and adjusted by experimental tests, for lack of a complete theoretical model.There have been some studies on how to optimize the reflow soldering temperature profile. Whalley used a simplified representation of products and processes to propose a simplified reflow soldering model [7]. A finite-difference method (FDM) based on an alternating-direction implicit scheme was a good way to calculate the heat transfer in vapor phase soldering [8]. In fact, the mainstream was the finite-element method (FEM) with advantages that there is no need to consider the complex differential equations and unintelligible physical concepts. While there was a fact that its accuracy was incomparable, the computation was particularly enormous [5, 9–13]. The surface temperature of the motherboard on flexible printed circuit boards was predicted using the simulation software ANSYS, though the reliability of the results still needs further experimental verification [14]. By applying simulation software, the computational fluid dynamics (CFD) models of reflow soldering were also tried [10–12, 15]. On the other hand, as the reflow soldering model becomes more accurate, on the basis of the temperature profile, the evaluation of soldering quality and optimization of relevant parameters, such as the heat transfer coefficient, heat flow, and reflux cooling time, have gradually been developed [5, 16]. For example, a linear regression model derived from the least-squares estimation methodology was utilized for describing the discovery that the heating factor Qη is approximately linear with changes in combined parameters Ht, providing a simple way for optimization based on the heating factor [4], which was a significant factor of the evaluation of solder joint quality and the bismuth/tin solder intermetallic layers too [6, 17]. The relationship between temperature, heat transfer coefficient and heat flux distribution, and cooling mass was investigated numerically; the optimal cooling period 11 s was obtained [16]. The application of the gray Taguchi method figured out the optimization problem of thermal stress and cooling rate of the solder joints of the ball net array package [18]. In addition, traditional artificial intelligence methods have also been applied to search for the thermal parameters of the reflow soldering process, and three alternative optimization methods were discussed [3].It is one of the difficulties in modeling owing to diverse ways of heat transfer. There was a literature using CFD software to model the infrared convection reflow oven and simulate the structure heating ball grid array (BGA) packages. The unexpected results indicated that the convection mode in the infrared convection oven had little effect on the heat transfer of the PCB. In other words, the outcomes of natural convection and forced convection were similar in an infrared-convection reflow oven [10]. Due to the huge temperature changes in the soldering process, the thermal conductivity and convection coefficients are inconstant. Therefore, research on the difference in the heat transfer coefficient of the printed circuit board in the reflow soldering process was meaningful [16, 19, 20]. Dziurdzia et al. performed statistical methods to evaluate convection reflow and vacuum vapor phase reflow, while the randomness of the distribution of voids in solder joints and vacuum convection reflow was a major disadvantage [21]. Furthermore, in the light of the former models, many intelligent algorithms, such as the genetic algorithm (GA) and the back-propagation neural network (BPNN), were utilized to seek the optimal parameters of the reflow oven [12]. For instance, Pan et al. put forward a parameter optimization model combined with the BPNN and GA to determine the best parameter setting for reflow soldering contours, contributing to reduce the trial-and-error time in practical applications [2]. Illés et al. provided a good summary of the previous models, as well as detailed discussions on the operating principle, heat transfer mechanisms, and numerical simulation of different reflow ovens [22]. Recently, the study on the nanocomposite solder paste attracted widespread attention again, since the discovery prevailed that adding nanoparticles enables to increase the solder properties in many terms, and it was proved to be true by numerical simulation [23].Note that micromorphological dynamics in solder alloys have been extensively studied through experiments and numerical methods in the past decades. In 2002, a model for phase separation driven by mechanical effects was proposed by Bonetti et al. [24]. Anders et al. employed an extended Cahn–Hilliard phase-field model to capture the essence of the microstructural evolution in solder alloys, and the numerical simulations were presented by the innovative isogeometric finite-element approach [25, 26]. We noted that a thermal coupling method for the BGA components in the forced convection reflow soldering process was put forth, showing that phase changes of the solder joints caused a wide range of temperature changes [27].Obviously, most of the previous studies in this area have been limited with costly experiments and simulation models with massive computation; especially, little research has been conducted on the effect of phase changes of solder on reflow soldering. In this work, we will focus on the rational control scheme through analyzing the heat transfer mechanism of reflow soldering, for studying the method of adjusting and controlling the temperature profile based on the theory of heat transfer and multiobjective optimization and also investigating the influence of phase changes on the temperature profile.


## 2. Mathematical Model and Computing Method

Mathematical models of reflow soldering and related computing methods are analyzed and developed based on the rigorous theory of heat transfer, which is the basis of the development of the dynamic Thomas algorithm in this work.


### 2.1. Heat Convection Model of Reflow Oven

Figure1 shows the internal structure of the reflow oven and the method of temperature measurement of the soldering area. The temperature distribution is mainly determined by the heating method, the temperature of each zone, and the size of the gap.Figure 1
Structure of reflow oven with 11 temperature zones.For the development of the mathematical model, we take a reflow oven with 11 temperature zones as an example. As shown in Figure1, zones 1∼5 are the preheating sections; zones 6∼7 are the constant temperature sections; zones 8∼9 are the recirculation sections; and zones 10∼11 are the cooling sections. Since the heat transfer characteristics of the heating and the cooling processes are different, they should be considered separately.Under normal circumstances of temperature measurement, some typical soldering areas should be selected to ensure the reliability of results, as shown in Figure2. The temperature-measurement equipment must be resistant to high temperatures since it has to pass through the reflow oven along with the PCB during the process of measurement. In addition, it should be noted that the probe must be small enough, generally less than 1 mm, to ensure that it can be embedded in the solder. Generally, the temperature measured in this way is basically the same as the solder itself.Figure 2
Temperature measurement of the soldering area.


#### 2.1.1. Heating Process

The heating process is mainly conducted in the range of zones 1∼9, while heat radiation and convection are the major heating methods. Since the height difference between two heating devices is small, each zone can be heated to its set temperature quickly, and this temperature can be kept unchanged. Therefore, we assume that the temperatures of zones 1∼9 are their set valuesux=unn=1,2,…,9. For the gap between two small zones, mass and heat transfer are realized by the forced convection of the influent gas in the furnace, so the heat transfer mode should contain the processes of convection and conduction. Generally, the soldering of the circuit board will be carried out only after the reflow oven is stable, and the temperature distribution in the reflow oven should be stable at this point. Therefore, it is basically a problem of steady-state heat transfer with no internal heat source, and the problem couples heat conduction and convection. Thus, the temperature distribution can be given by the Fourier–Kirchhoff equation as follows [28]:(1)a∇2u−vx∂u∂x+vy∂u∂y+vz∂u∂z=∂u∂t,where a=k/cρ, c is the specific heat at a constant pressure of influent gas, ρ is the density of influent gas, k is the thermal conductivity of influent gas, u is the temperature, and v is the flow velocity of influent gas.By considering the symmetry of the temperature field, the temperature distribution in the reflow oven is only related to the horizontal coordinates, instead of the vertical coordinates. Since the medium in the reflow oven is usually an inert gas, and its physical and chemical properties are relatively stable, the related thermal parameters of the medium can be considered as a constant. Thus, the issue of steady-state heat transfer in the reflow oven is as follows:(2)a∂2u∂2x−vx∂u∂x=0,u|x=xn,R=un,u|x=xn+1,L=un+1,where n=1,2,…,9, xn,R is the coordinates of the right boundary of the n-th small temperature region, and xn+1,L is the coordinates of the left edge of the n+1-th region. The analytic solution of the equation is as follows:(3)u=A1n+A2nevxx/a,where(4)A1n=unevxxn+1−xn/a−un+1evxxn+1−xn/a−1,A2n=un+1−unevxxn+1/a−evxxn/a.Obviously, when the two boundary conditions have the same temperatureun=un+1, the solution of the equation is u=un, which means the temperature of the gap is equal to that of its adjacent zones.Generally, there is natural convection instead of forced convection in the front area, so the boundary condition of free cooling is adopted on the left side. Hence, there is:(5)a∂2u∂2x−vx∂u∂x=0,k∂u∂x|x=0=hu|x=0−ua,u|x=x0,R=u1,where h is the natural convection heat transfer coefficient of influent gas, x0,R is the location coordinate of the right end of the front area, and ua is the ambient temperature. As a result, the temperature distribution of the front area is as follows:(6)u=B1n+B2nevvx/a,where(7)B1n=u1−u1−uaevxx0,L/aevxx0,L/a−1+kvx/ha,B2n=u1−uaevxx0,L/a−1+kvx/ha.


#### 2.1.2. Cooling Process

In the course of cooling, the high temperature in the recirculating section will restrict the cooling effect, so that the temperature in the cooling section is incapable of reaching its set temperature. As a result, the cooling zones and the back area of the furnace should be considered as a whole. The heat transfer problem of the cooling process is then as follows:(8)a∂2u∂2x−vx∂u∂x=0,u|x=x9,R=u9,u|x=xend=u11.The analytic solution of this equation is as follows:(9)u=C1n+C2nevxx/a,where(10)C1n=u11−u9evxxb/a−evxx9/a,C2n=u9evxxend−x9/a−u11evxxend−x9/a−1.


### 2.2. Heat Conduction Model of Soldering Area

The soldering area is generally distributed on the surface of the PCB, with few or no solder joints on the boundary. Its internal heat transfer mode should be based on the process of heat conduction, since the PCB and soldering area are both solid objects. In view of its small thickness, the temperature change of the soldering area is mainly affected by the convection heat transfer of the upside and downside gas media but with little effects resulting from the horizontal heat conduction.


#### 2.2.1. Without Phase Transformations

Usually, the solder material of the circuit board is Sn96.5Ag3Cu0.5, and its melting point is about 217°C. Before reaching the melting point, the heat transfer of the soldering area is a nonsteady heat conduction problem without an internal heat source, and its temperature can be calculated using the following equation:(11)∂T∂t−λ∂2T∂2y=0,where λ=k′/c′ρ′, c′ is the constant pressure specific heat capacity of the circuit board, ρ′ is the density of the circuit board, and k′ is the thermal conductivity of the circuit board. The initial condition is T|t=0=25∘C. As far as the boundary condition is concerned, since the heat transfer between the soldering area and the reflow oven is dominated by the process of heat convection, the third kind of boundary condition is adopted. The amount of heat convection can then be given by Newton’s law of cooling [29]. Hence, Q1=kcux0+vt−T|y=b and Q2=kcux0+vt−T|y=0. Then, according to Fourier’s law of heat conduction [29], the heat flowing in from the boundaries can be expressed as Q1=k∂T/∂y|y=b and Q2=−k∂T/∂y|y=0. Since the quantity of heat presented in the two ways should be equal, we have the following equation:(12)k∂T∂y|y=b=kcux0+vt−T|y=b,−k∂T∂y|y=0=kcux0+vt−T|y=0,where kc represents the convective heat transfer coefficient between the circuit board and the gas in the furnace, x0 is the position coordinates of the soldering area when the circuit board in its initial position, b is the thickness of the PCB, v stands for the passing speed of the conveyor belt, and ux0+vt denotes the temperature at x0+vt in reflow oven.There’s a huge difference in temperature distribution. Generally, the temperature of the recirculating zone is about 255°C, while the cooling zone is about 25°C. As a result,k and kc are functions of temperature instead of constants approximately. Different simplification methods lead to different results slightly. Theoretically, the more the coefficient, the higher the degree of fitting between the results and the reference data, but multifarious coefficients will lead to overfitting and huge computation [30].


#### 2.2.2. With Phase Transformations

One of the two processes for the solder to change phases is melting, where it changes from solid to liquid phase, which requires heat absorption. Figure3 displays the heat transfer process in the soldering area when the solder melts, in which heat flows into the solid phase passing through the liquid phase, from the upper and lower surfaces. Meanwhile, the latent heat of phase transitions is absorbed in the interface. YLt and YHt represent the two interfaces of the solid and the liquid phases, which change with time from both ends to the middle gradually, and finally, the soldering area becomes liquid phase.Figure 3
Phase transformations in thawing process.There is no forced convection in molten solder, and natural convection has little effect on heat conduction. Therefore, it can be assumed that the heat transfer model of liquid phase is based on the process of heat conduction. Thus, the temperature distribution of the two phases should meet the heat conduction equations, respectively.(13)∂Tl∂t−λl∂2Tl∂2y=0,0<y<YLtorYHt<y<b,∂Ts∂t−λs∂2Ts∂2y=0,YLt<y<YHt,where λl=kl/clρl and λl=ks/csρs; Tl and Ts represent the temperature distribution of the liquid and solid phases, respectively; cl and cs denote the constant pressure heat capacity of liquid and solid solder, respectively; ρl and ρs signify the density of liquid and solid solder, respectively; kl and ks indicate the thermal conductivity of liquid and solid solder, respectively. The temperature at the interface is the phase-transition temperature. According to the law of conservation of energy, the heat flux at the inflow interface minus that at the outflow interface is equal to the latent heat of phase changes absorbed per unit time [31]. Therefore:(14)y=YLt:ql−qs=ρsLdYLtdt,y=YHt:−ql+qs=−ρsLdYHtdt,where ql and qs represent the heat flux in the positive direction of the liquid and solid phases, respectively, and L is the latent heat of solder melting. According to Fourier’s law [29], Equation (14) can be interpreted as follows:(15)y=YLt:−kl∂Tl∂y+ks∂Ts∂y=ρsLdYLtdt,y=YHt:kl∂Tl∂y−ks∂Ts∂y=−ρsLdYHtdt.The other issue is the cooling process, that is, the solder changes from liquid to solid phase as the temperature decreases. Figure4 shows the heat transfer process in the soldering area of the cooling process. Heat flows out from the liquid phase to the outside, passing through the solid phase and two surfaces. At the same time, the latent heat of phase changes is released at the interface. Interfaces between the solid and the liquid phases get closer gradually from the two sides towards the middle, and the soldering area becomes the solid phase eventually.Figure 4
Phase transformations in the cooling process.Similarly, the temperature distribution of the two phases ought to meet the heat conduction equation (Equation (13)), but the solid and the liquid phases are interchanged. Then, the temperature at the interfaces can also be deduced from the law of conservation of energy and Fourier’s law [29] as follows:(16)y=YLt:kl∂Tl∂y−ks∂Ts∂y=−ρlLdYLtdt,y=YHt:−kl∂Tl∂y+ks∂Ts∂y=ρlLdYHtdt.The third boundary condition of heat conduction is still adopted for the upper and lower boundaries of the soldering area in the phase transition process (Equation (12)).The temperature distributionux of each zone calculated; the heat conduction problem can be resolved to obtain the temperature variation Ty,t of the soldering area. Then, the temperature variation of the center of the soldering area is Tb/2,t, namely, the temperature profile.


### 2.3. A Dynamic Thomas Algorithm

FDM and FEM are the two most popular approaches to solve differential equations. The advantage of the FDM is its simple principle, but it can only get the numerical solution in the regular region. Instead, the FEM is inclined to solve the irregular region problems [32]. Since the boundary of the mathematical model in our study is a regular rectangular region, the finite-difference method has been adopted.


#### 2.3.1. Ignoring Phase Transformations

Firstly, the heat conduction equation of the soldering area was solved by the implicit difference method, according to Equation (11), as follows:(17)Ti,j−1=αTi−1,j+1−2αTi,j+αTi+1,j,where α=kΔt/cρΔy2=λΔt/Δy2, i=y/Δy=0,1,2,…,m, and j=t/Δt=0,1,2,…,n. The boundary conditions (Equation (12)) were discretized using backward difference for the upper and forward difference for the lower.(18)1+μTm,j−Tm−1,j=μux0+vt,1+μT0,j−T1,j=μux0+vt,where μ=kcΔy/k. The initial condition is Ti,0=25. Equation (17) is an implicit difference scheme for a heat conduction model of a soldering region without phase changes, which can be treated as a tridiagonal system of linear equations with respect to Ti,j. Substitute the boundary conditions, and the matrix form can be expressed as follows:(19)1+2α−α1+μ−α0⋯0−α1+2α−α⋯0⋮⋱⋱⋱⋮0⋯−α1+2α−α0⋯0−α1+2α−α1+μT1,jT2,j⋮Tm−2,jTm−1,j=T1,j−1+UjT2,j⋮Tm−2,j−1Tm−1,j−1+Uj,where Uj=αμ/1+μux0+vjΔt. The system of linear equations can be recorded as AX=B, and the nonzero elements in the coefficient matrix A are distributed on the leading diagonal and the two pairs of adjacent diagonals above and below. This form of linear equations can be solved quickly by the forward elimination and backward substitution method, namely, the Thomas algorithm [32].


#### 2.3.2. Considering Phase Transitions

The heat conduction equation and boundary conditions are similar to the model without phase transitions, but the difference is that the initial condition should be the melting point of the solder. Consequently, the following are the equations for the melting process:(20)Ti,j−1l=αlTi−1,jl+1−2αlTi,jl+αlTi+1,jl,Ti,j−1s=αsTi−1,js+1−2αsTi,js+αsTi+1,js,Ti,0=Tmelting point,1+μlTm,j−Tm−1,j=μlux0+vt,1+μlT0,j−T1,j=μlux0+vt,where μl=kclΔy/kl, αl=klΔt/clρlΔy2=λlΔt/Δy2, and αs=ksΔt/csρsΔy2=λsΔt/Δy2. Tl is the temperature of the liquid phase, and Ts is the temperature of the solid phase. The cooling process is opposite, and equations are similar.(21)Ti,j−1s=αsTi−1,js+1−2αsTi,js+αsTi+1,js,Ti,j−1l=αlTi−1,jl+1−2αlTi,jl+αlTi+1,jl,Ti,0=Tfreezing point,1+μsTm,j−Tm−1,j=μsux0+vt,1+μsT0,j−T1,j=μsux0+vt,where μs=kcsΔy/kl. When phases changing, the initial temperature of the solid and the liquid phases are supposed to satisfy the heat transfer equation with the latent heat of phase transitions. For the boundaries (Equation (15)), the liquid phase adopted forward difference, and the solid phase adopted backward difference at the upper boundary, and it is just the reverse at the lower boundary. Hence, the interface conditions of phase transitions for the melting process are as follows:(22)−klTYL,j−TYL−1,j+ksTYL+1,j−TYL,j=ρsLvLΔy,klTYH+1,j−TYH,j−ksTYH,j−TYH−1,j=−ρsLvHΔy,where vL=ΔYL/Δt and vH=ΔYH/Δt. Similarly, the interface conditions during cooling are as follows:(23)klTYL,j−TYL−1,j−ksTYL+1,j−TYL,j=−ρlLvLΔy,−klTYH+1,j−TYH,j+ksTYH,j−TYH−1,j=ρlLvHΔy.For the heat conduction model of the phase transition process, the solid and the liquid regions are the same as the model without phase transitions, yet there are two more dynamic boundary conditions, so the soldering area should be divided into three sections for calculations as follows:(24)1+2α1−α11+μ−α10⋯0−α11+2α1−α1⋯0⋮⋱⋱⋱⋮0⋯−α11+2α1−α10⋯0−α11+2α1T1,jT2,j⋮TYL−2,jTYL−1,j=T1,j−1+UjT2,j⋮TYL−2,j−1TYL−1,j−1+αTYL,j,1+2α2−α20⋯0−α21+2α2−α2⋯0⋮⋱⋱⋱⋮0⋯−α21+2α2−α20⋯0−α21+2α2TYL+1,jTYL+2,j⋮TYH−2,jTYH−1,j=TYL+1,j−1+αTYL,jTYL+2,j⋮TYH−2,j−1TYH−1,j−1+αTYH,j,1+2α1−α10⋯0−α11+2α1−α⋯0⋮⋱⋱⋱⋮0⋯−α1+2α1−α10⋯0−α11+2α1−α11+μTYH+1,jTYH+2,j⋮Tm−2,jTm−1,j=TYH+1,j−1+αTYL,jTYH+2,j⋮Tm−2,j−1Tm−1,j−1+Uj,where μ=μl, α1=αl, α2=αs, and Uj=αlμl/1+μlux0+vjΔt when melting, and μ=μs, α1=αs, α2=αl, and Uj=αsμs/1+μsux0+vjΔt when cooling. Substitute the transformation boundary conditions (Equations (22) and (23)) that TYL,j=−ρlLvLΔy+klTL−1,j+ksTL+1,j/kl+ks and TYH,j=−ρlLvHΔy+klTH+1,j+ksTH−1,j/kl+ks during melting and TYL,j=−ρsLvLΔy+klTL−1,j+ksTL+1,j/kl+ks and TYH,j=−ρsLvHΔy+klTH+1,j+ksTH−1,j/kl+ks during cooling, and the following relation can be obtained:(25)1+2α1−α11+μ1−α1−α11+2α1−α1⋱⋱⋱−α11+2α1klks+kl−α1ksks+kl−α2klks+kl1+2α2ksks+kl−α2⋱⋱⋱−α21+2α2−α2⋱⋱⋱−α21+2α2ksks+kl−α2klks+kl−α1ksks+kl1+2α1klks+kl−α1⋱⋱⋱−α11+2α1−α1−α11+2α1−α11+μT1,jT2,j⋮TYL−1,jTYL+1,j⋮Ti,j⋮TYH−1,jTYH+1,j⋮Tm−2,jTm−1,j=T1,j−1+UjT2,j⋮TYL−1,j−1−LjTYL+1,j−1−Lj⋮Ti,j−1⋮TYH−1,j−1−HjTYH+1,j−1−Hj⋮Tm−2,j−1Tm−1,j−1+Uj,where Lj=αlρlLvLΔy/kskl and Hj=αlρlLvHΔy/kskl during the melting process, and Lj=αlρlLvLΔy/kskl and Hj=αlρlLvHΔy/kskl during the cooling process. The coefficient matrix of equations without phase transformations is independent of time and can be solved quickly by the Thomas algorithm. As for the phase transition process, the coefficient matrix changes with time because interfaces of the solid and the liquid phases approach gradually. Consequently, Equation (25) can be expressed as AtXt=Yt. At, Xt and Yt can be written in the following form:(26)At=b1tc1t0⋯0a2tb2tc2t⋯0⋮⋱⋱⋱⋮0⋯am−4tbm−4tcm−4t0⋯0am−3tbm−3t,Xt=x1tx2t⋮xm−4txm−3t,Yt=y1ty2t⋮ym−4tym−3t.Although the coefficient matrix changes with time, the system still satisfies the condition of the Crout decomposition method, that is, (i)b1t≥c1t≥0, (ii) bit≥ait+cit, aitcit≠0, i=2,3,…,m−3, and (iii) bm−3t≥am−3t≥0. Obviously, At is a nonsingular matrix whose determinant detAt≠0. Thus, its sequential principal minor of any order Dkt≠0.A tridiagonal square matrix satisfying the above conditions has a unique Crout decomposition matrix as follows:(27)At=PtQt=p1t00⋯0a2tp2t0⋯0⋮⋱⋱⋱⋮0⋯am−4tpm−4t00⋯0am−3tpm−3t1q1t0⋯001q2t⋯0⋮⋱⋱⋱⋮0⋯01qm−4t0⋯001,where p1t=b1t, qit=cit/piti=1,2,…,m−4, and pit=bit−aitqi−1ti=2,3,…,m−3. pit and qit are time-varying and piecewise functions, and the recurrence formulas can be expressed as follows:(28)qit=−α11+2α1−α1/1+μ,i=1;−α11+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;−α1ks/ks+kl1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;−α21+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;−α21+2α2+α2qi−1t,YLt≤i≤YHt−3;−α2kl/ks+kl1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;−α11+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1.pit=1+2α1−α11+μ,i=1;1+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;1+2α1klks+kl+α1qi−1t,i=YLt−1;1+2α2ksks+kl+α2klks+klqi−1t,i=YLt;1+2α2+α2qi−1t,YLt≤i≤YHt−3;1+2α2kSks+kl+α2qi−1t,i=YHt−2;1+2α2klks+kl+α1ksks+klqi−1t,i=YHt−1;1+2α1−α11+μ+α1qi−1t,i=m−3.Therefore,AtXt=PtQtXt=Yt can be decomposed into two triangular equations PtDt=Yt and QtXt=Dt. The first step is to solve the system of equations PtDt=Yt as follows:(29)p1t00⋯0a2tp2t0⋯0⋮⋱⋱⋱⋮0⋯am−4tpm−4t00⋯0am−3tpm−3td1td2t⋮dm−4tdm−3t=y1ty2t⋮ym−4tym−3t.The recurrence formulas of solutions are as follows:(30)dit=y1t1+2α1−α1/1+μ,i=1;yit+α1di−1t1+2α1+α1qi−1t2≤i≤YLt−2orYHt≤i≤m−4;yit+α1di−1t1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;yit+α2kl/ks+kldi−1t1+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;yit+α2di−1t1+2α2+α2qi−1t,YLt≤i≤YHt−3;yit+α2di−1t1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;yit+α1ks/ks+kldi−1t1+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1;yit+α1di−1t1+2α1−α1/1+μ+α1qi−1t,i=m−3.The second step is to solve the system of equationsQtXt=Dt as follows:(31)1q1t0⋯001q2t⋯0⋮⋱⋱⋱⋮0⋯01qm−4t0⋯001x1tx2t⋮xm−4txm−3t=d1td2t⋮dm−4tdm−3t.The recurrence formulas of solutions are as follows:(32)xit=yit+α1di−1t1+2α1−α1/1+μ+α1qi−1t,i=m−3;yit+α1ks/ks+kldi−1t+α1xi+1t1+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1;yit+α2di−1t+α1kl/ks+klxi+1t1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;yit+α2di−1t+α2xi+1t1+2α2+α2qi−1t,YLt≤i≤YHt−3;yit+α2kl/ks+kldi−1t+α2xi+1t1+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;yit+α1di−1t+α1ks/ks+klxi+1t1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;yit+α1di−1t+α1xi+1t1+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;y1t+α1xi+1t1+2α1−α1/1+μ,i=1.wherexit is the discrete solution of the heat transfer equation in the case of considering phase transitions. A dynamic Thomas algorithm is proposed and developed in this study, in order to address the time-dependent solution of the tridiagonal linear system. Flowchart of the developed algorithm is indicated in Figure 5. By cyclic iteration, the temperature of the soldering area during the entire process can be obtained.Figure 5
Flowchart of the dynamic Thomas algorithm, that is, forward elimination and backward substitution with time.


## 2.1. Heat Convection Model of Reflow Oven

Figure1 shows the internal structure of the reflow oven and the method of temperature measurement of the soldering area. The temperature distribution is mainly determined by the heating method, the temperature of each zone, and the size of the gap.Figure 1
Structure of reflow oven with 11 temperature zones.For the development of the mathematical model, we take a reflow oven with 11 temperature zones as an example. As shown in Figure1, zones 1∼5 are the preheating sections; zones 6∼7 are the constant temperature sections; zones 8∼9 are the recirculation sections; and zones 10∼11 are the cooling sections. Since the heat transfer characteristics of the heating and the cooling processes are different, they should be considered separately.Under normal circumstances of temperature measurement, some typical soldering areas should be selected to ensure the reliability of results, as shown in Figure2. The temperature-measurement equipment must be resistant to high temperatures since it has to pass through the reflow oven along with the PCB during the process of measurement. In addition, it should be noted that the probe must be small enough, generally less than 1 mm, to ensure that it can be embedded in the solder. Generally, the temperature measured in this way is basically the same as the solder itself.Figure 2
Temperature measurement of the soldering area.


### 2.1.1. Heating Process

The heating process is mainly conducted in the range of zones 1∼9, while heat radiation and convection are the major heating methods. Since the height difference between two heating devices is small, each zone can be heated to its set temperature quickly, and this temperature can be kept unchanged. Therefore, we assume that the temperatures of zones 1∼9 are their set valuesux=unn=1,2,…,9. For the gap between two small zones, mass and heat transfer are realized by the forced convection of the influent gas in the furnace, so the heat transfer mode should contain the processes of convection and conduction. Generally, the soldering of the circuit board will be carried out only after the reflow oven is stable, and the temperature distribution in the reflow oven should be stable at this point. Therefore, it is basically a problem of steady-state heat transfer with no internal heat source, and the problem couples heat conduction and convection. Thus, the temperature distribution can be given by the Fourier–Kirchhoff equation as follows [28]:(1)a∇2u−vx∂u∂x+vy∂u∂y+vz∂u∂z=∂u∂t,where a=k/cρ, c is the specific heat at a constant pressure of influent gas, ρ is the density of influent gas, k is the thermal conductivity of influent gas, u is the temperature, and v is the flow velocity of influent gas.By considering the symmetry of the temperature field, the temperature distribution in the reflow oven is only related to the horizontal coordinates, instead of the vertical coordinates. Since the medium in the reflow oven is usually an inert gas, and its physical and chemical properties are relatively stable, the related thermal parameters of the medium can be considered as a constant. Thus, the issue of steady-state heat transfer in the reflow oven is as follows:(2)a∂2u∂2x−vx∂u∂x=0,u|x=xn,R=un,u|x=xn+1,L=un+1,where n=1,2,…,9, xn,R is the coordinates of the right boundary of the n-th small temperature region, and xn+1,L is the coordinates of the left edge of the n+1-th region. The analytic solution of the equation is as follows:(3)u=A1n+A2nevxx/a,where(4)A1n=unevxxn+1−xn/a−un+1evxxn+1−xn/a−1,A2n=un+1−unevxxn+1/a−evxxn/a.Obviously, when the two boundary conditions have the same temperatureun=un+1, the solution of the equation is u=un, which means the temperature of the gap is equal to that of its adjacent zones.Generally, there is natural convection instead of forced convection in the front area, so the boundary condition of free cooling is adopted on the left side. Hence, there is:(5)a∂2u∂2x−vx∂u∂x=0,k∂u∂x|x=0=hu|x=0−ua,u|x=x0,R=u1,where h is the natural convection heat transfer coefficient of influent gas, x0,R is the location coordinate of the right end of the front area, and ua is the ambient temperature. As a result, the temperature distribution of the front area is as follows:(6)u=B1n+B2nevvx/a,where(7)B1n=u1−u1−uaevxx0,L/aevxx0,L/a−1+kvx/ha,B2n=u1−uaevxx0,L/a−1+kvx/ha.


### 2.1.2. Cooling Process

In the course of cooling, the high temperature in the recirculating section will restrict the cooling effect, so that the temperature in the cooling section is incapable of reaching its set temperature. As a result, the cooling zones and the back area of the furnace should be considered as a whole. The heat transfer problem of the cooling process is then as follows:(8)a∂2u∂2x−vx∂u∂x=0,u|x=x9,R=u9,u|x=xend=u11.The analytic solution of this equation is as follows:(9)u=C1n+C2nevxx/a,where(10)C1n=u11−u9evxxb/a−evxx9/a,C2n=u9evxxend−x9/a−u11evxxend−x9/a−1.


## 2.1.1. Heating Process

The heating process is mainly conducted in the range of zones 1∼9, while heat radiation and convection are the major heating methods. Since the height difference between two heating devices is small, each zone can be heated to its set temperature quickly, and this temperature can be kept unchanged. Therefore, we assume that the temperatures of zones 1∼9 are their set valuesux=unn=1,2,…,9. For the gap between two small zones, mass and heat transfer are realized by the forced convection of the influent gas in the furnace, so the heat transfer mode should contain the processes of convection and conduction. Generally, the soldering of the circuit board will be carried out only after the reflow oven is stable, and the temperature distribution in the reflow oven should be stable at this point. Therefore, it is basically a problem of steady-state heat transfer with no internal heat source, and the problem couples heat conduction and convection. Thus, the temperature distribution can be given by the Fourier–Kirchhoff equation as follows [28]:(1)a∇2u−vx∂u∂x+vy∂u∂y+vz∂u∂z=∂u∂t,where a=k/cρ, c is the specific heat at a constant pressure of influent gas, ρ is the density of influent gas, k is the thermal conductivity of influent gas, u is the temperature, and v is the flow velocity of influent gas.By considering the symmetry of the temperature field, the temperature distribution in the reflow oven is only related to the horizontal coordinates, instead of the vertical coordinates. Since the medium in the reflow oven is usually an inert gas, and its physical and chemical properties are relatively stable, the related thermal parameters of the medium can be considered as a constant. Thus, the issue of steady-state heat transfer in the reflow oven is as follows:(2)a∂2u∂2x−vx∂u∂x=0,u|x=xn,R=un,u|x=xn+1,L=un+1,where n=1,2,…,9, xn,R is the coordinates of the right boundary of the n-th small temperature region, and xn+1,L is the coordinates of the left edge of the n+1-th region. The analytic solution of the equation is as follows:(3)u=A1n+A2nevxx/a,where(4)A1n=unevxxn+1−xn/a−un+1evxxn+1−xn/a−1,A2n=un+1−unevxxn+1/a−evxxn/a.Obviously, when the two boundary conditions have the same temperatureun=un+1, the solution of the equation is u=un, which means the temperature of the gap is equal to that of its adjacent zones.Generally, there is natural convection instead of forced convection in the front area, so the boundary condition of free cooling is adopted on the left side. Hence, there is:(5)a∂2u∂2x−vx∂u∂x=0,k∂u∂x|x=0=hu|x=0−ua,u|x=x0,R=u1,where h is the natural convection heat transfer coefficient of influent gas, x0,R is the location coordinate of the right end of the front area, and ua is the ambient temperature. As a result, the temperature distribution of the front area is as follows:(6)u=B1n+B2nevvx/a,where(7)B1n=u1−u1−uaevxx0,L/aevxx0,L/a−1+kvx/ha,B2n=u1−uaevxx0,L/a−1+kvx/ha.


## 2.1.2. Cooling Process

In the course of cooling, the high temperature in the recirculating section will restrict the cooling effect, so that the temperature in the cooling section is incapable of reaching its set temperature. As a result, the cooling zones and the back area of the furnace should be considered as a whole. The heat transfer problem of the cooling process is then as follows:(8)a∂2u∂2x−vx∂u∂x=0,u|x=x9,R=u9,u|x=xend=u11.The analytic solution of this equation is as follows:(9)u=C1n+C2nevxx/a,where(10)C1n=u11−u9evxxb/a−evxx9/a,C2n=u9evxxend−x9/a−u11evxxend−x9/a−1.


## 2.2. Heat Conduction Model of Soldering Area

The soldering area is generally distributed on the surface of the PCB, with few or no solder joints on the boundary. Its internal heat transfer mode should be based on the process of heat conduction, since the PCB and soldering area are both solid objects. In view of its small thickness, the temperature change of the soldering area is mainly affected by the convection heat transfer of the upside and downside gas media but with little effects resulting from the horizontal heat conduction.


### 2.2.1. Without Phase Transformations

Usually, the solder material of the circuit board is Sn96.5Ag3Cu0.5, and its melting point is about 217°C. Before reaching the melting point, the heat transfer of the soldering area is a nonsteady heat conduction problem without an internal heat source, and its temperature can be calculated using the following equation:(11)∂T∂t−λ∂2T∂2y=0,where λ=k′/c′ρ′, c′ is the constant pressure specific heat capacity of the circuit board, ρ′ is the density of the circuit board, and k′ is the thermal conductivity of the circuit board. The initial condition is T|t=0=25∘C. As far as the boundary condition is concerned, since the heat transfer between the soldering area and the reflow oven is dominated by the process of heat convection, the third kind of boundary condition is adopted. The amount of heat convection can then be given by Newton’s law of cooling [29]. Hence, Q1=kcux0+vt−T|y=b and Q2=kcux0+vt−T|y=0. Then, according to Fourier’s law of heat conduction [29], the heat flowing in from the boundaries can be expressed as Q1=k∂T/∂y|y=b and Q2=−k∂T/∂y|y=0. Since the quantity of heat presented in the two ways should be equal, we have the following equation:(12)k∂T∂y|y=b=kcux0+vt−T|y=b,−k∂T∂y|y=0=kcux0+vt−T|y=0,where kc represents the convective heat transfer coefficient between the circuit board and the gas in the furnace, x0 is the position coordinates of the soldering area when the circuit board in its initial position, b is the thickness of the PCB, v stands for the passing speed of the conveyor belt, and ux0+vt denotes the temperature at x0+vt in reflow oven.There’s a huge difference in temperature distribution. Generally, the temperature of the recirculating zone is about 255°C, while the cooling zone is about 25°C. As a result,k and kc are functions of temperature instead of constants approximately. Different simplification methods lead to different results slightly. Theoretically, the more the coefficient, the higher the degree of fitting between the results and the reference data, but multifarious coefficients will lead to overfitting and huge computation [30].


### 2.2.2. With Phase Transformations

One of the two processes for the solder to change phases is melting, where it changes from solid to liquid phase, which requires heat absorption. Figure3 displays the heat transfer process in the soldering area when the solder melts, in which heat flows into the solid phase passing through the liquid phase, from the upper and lower surfaces. Meanwhile, the latent heat of phase transitions is absorbed in the interface. YLt and YHt represent the two interfaces of the solid and the liquid phases, which change with time from both ends to the middle gradually, and finally, the soldering area becomes liquid phase.Figure 3
Phase transformations in thawing process.There is no forced convection in molten solder, and natural convection has little effect on heat conduction. Therefore, it can be assumed that the heat transfer model of liquid phase is based on the process of heat conduction. Thus, the temperature distribution of the two phases should meet the heat conduction equations, respectively.(13)∂Tl∂t−λl∂2Tl∂2y=0,0<y<YLtorYHt<y<b,∂Ts∂t−λs∂2Ts∂2y=0,YLt<y<YHt,where λl=kl/clρl and λl=ks/csρs; Tl and Ts represent the temperature distribution of the liquid and solid phases, respectively; cl and cs denote the constant pressure heat capacity of liquid and solid solder, respectively; ρl and ρs signify the density of liquid and solid solder, respectively; kl and ks indicate the thermal conductivity of liquid and solid solder, respectively. The temperature at the interface is the phase-transition temperature. According to the law of conservation of energy, the heat flux at the inflow interface minus that at the outflow interface is equal to the latent heat of phase changes absorbed per unit time [31]. Therefore:(14)y=YLt:ql−qs=ρsLdYLtdt,y=YHt:−ql+qs=−ρsLdYHtdt,where ql and qs represent the heat flux in the positive direction of the liquid and solid phases, respectively, and L is the latent heat of solder melting. According to Fourier’s law [29], Equation (14) can be interpreted as follows:(15)y=YLt:−kl∂Tl∂y+ks∂Ts∂y=ρsLdYLtdt,y=YHt:kl∂Tl∂y−ks∂Ts∂y=−ρsLdYHtdt.The other issue is the cooling process, that is, the solder changes from liquid to solid phase as the temperature decreases. Figure4 shows the heat transfer process in the soldering area of the cooling process. Heat flows out from the liquid phase to the outside, passing through the solid phase and two surfaces. At the same time, the latent heat of phase changes is released at the interface. Interfaces between the solid and the liquid phases get closer gradually from the two sides towards the middle, and the soldering area becomes the solid phase eventually.Figure 4
Phase transformations in the cooling process.Similarly, the temperature distribution of the two phases ought to meet the heat conduction equation (Equation (13)), but the solid and the liquid phases are interchanged. Then, the temperature at the interfaces can also be deduced from the law of conservation of energy and Fourier’s law [29] as follows:(16)y=YLt:kl∂Tl∂y−ks∂Ts∂y=−ρlLdYLtdt,y=YHt:−kl∂Tl∂y+ks∂Ts∂y=ρlLdYHtdt.The third boundary condition of heat conduction is still adopted for the upper and lower boundaries of the soldering area in the phase transition process (Equation (12)).The temperature distributionux of each zone calculated; the heat conduction problem can be resolved to obtain the temperature variation Ty,t of the soldering area. Then, the temperature variation of the center of the soldering area is Tb/2,t, namely, the temperature profile.


## 2.2.1. Without Phase Transformations

Usually, the solder material of the circuit board is Sn96.5Ag3Cu0.5, and its melting point is about 217°C. Before reaching the melting point, the heat transfer of the soldering area is a nonsteady heat conduction problem without an internal heat source, and its temperature can be calculated using the following equation:(11)∂T∂t−λ∂2T∂2y=0,where λ=k′/c′ρ′, c′ is the constant pressure specific heat capacity of the circuit board, ρ′ is the density of the circuit board, and k′ is the thermal conductivity of the circuit board. The initial condition is T|t=0=25∘C. As far as the boundary condition is concerned, since the heat transfer between the soldering area and the reflow oven is dominated by the process of heat convection, the third kind of boundary condition is adopted. The amount of heat convection can then be given by Newton’s law of cooling [29]. Hence, Q1=kcux0+vt−T|y=b and Q2=kcux0+vt−T|y=0. Then, according to Fourier’s law of heat conduction [29], the heat flowing in from the boundaries can be expressed as Q1=k∂T/∂y|y=b and Q2=−k∂T/∂y|y=0. Since the quantity of heat presented in the two ways should be equal, we have the following equation:(12)k∂T∂y|y=b=kcux0+vt−T|y=b,−k∂T∂y|y=0=kcux0+vt−T|y=0,where kc represents the convective heat transfer coefficient between the circuit board and the gas in the furnace, x0 is the position coordinates of the soldering area when the circuit board in its initial position, b is the thickness of the PCB, v stands for the passing speed of the conveyor belt, and ux0+vt denotes the temperature at x0+vt in reflow oven.There’s a huge difference in temperature distribution. Generally, the temperature of the recirculating zone is about 255°C, while the cooling zone is about 25°C. As a result,k and kc are functions of temperature instead of constants approximately. Different simplification methods lead to different results slightly. Theoretically, the more the coefficient, the higher the degree of fitting between the results and the reference data, but multifarious coefficients will lead to overfitting and huge computation [30].


## 2.2.2. With Phase Transformations

One of the two processes for the solder to change phases is melting, where it changes from solid to liquid phase, which requires heat absorption. Figure3 displays the heat transfer process in the soldering area when the solder melts, in which heat flows into the solid phase passing through the liquid phase, from the upper and lower surfaces. Meanwhile, the latent heat of phase transitions is absorbed in the interface. YLt and YHt represent the two interfaces of the solid and the liquid phases, which change with time from both ends to the middle gradually, and finally, the soldering area becomes liquid phase.Figure 3
Phase transformations in thawing process.There is no forced convection in molten solder, and natural convection has little effect on heat conduction. Therefore, it can be assumed that the heat transfer model of liquid phase is based on the process of heat conduction. Thus, the temperature distribution of the two phases should meet the heat conduction equations, respectively.(13)∂Tl∂t−λl∂2Tl∂2y=0,0<y<YLtorYHt<y<b,∂Ts∂t−λs∂2Ts∂2y=0,YLt<y<YHt,where λl=kl/clρl and λl=ks/csρs; Tl and Ts represent the temperature distribution of the liquid and solid phases, respectively; cl and cs denote the constant pressure heat capacity of liquid and solid solder, respectively; ρl and ρs signify the density of liquid and solid solder, respectively; kl and ks indicate the thermal conductivity of liquid and solid solder, respectively. The temperature at the interface is the phase-transition temperature. According to the law of conservation of energy, the heat flux at the inflow interface minus that at the outflow interface is equal to the latent heat of phase changes absorbed per unit time [31]. Therefore:(14)y=YLt:ql−qs=ρsLdYLtdt,y=YHt:−ql+qs=−ρsLdYHtdt,where ql and qs represent the heat flux in the positive direction of the liquid and solid phases, respectively, and L is the latent heat of solder melting. According to Fourier’s law [29], Equation (14) can be interpreted as follows:(15)y=YLt:−kl∂Tl∂y+ks∂Ts∂y=ρsLdYLtdt,y=YHt:kl∂Tl∂y−ks∂Ts∂y=−ρsLdYHtdt.The other issue is the cooling process, that is, the solder changes from liquid to solid phase as the temperature decreases. Figure4 shows the heat transfer process in the soldering area of the cooling process. Heat flows out from the liquid phase to the outside, passing through the solid phase and two surfaces. At the same time, the latent heat of phase changes is released at the interface. Interfaces between the solid and the liquid phases get closer gradually from the two sides towards the middle, and the soldering area becomes the solid phase eventually.Figure 4
Phase transformations in the cooling process.Similarly, the temperature distribution of the two phases ought to meet the heat conduction equation (Equation (13)), but the solid and the liquid phases are interchanged. Then, the temperature at the interfaces can also be deduced from the law of conservation of energy and Fourier’s law [29] as follows:(16)y=YLt:kl∂Tl∂y−ks∂Ts∂y=−ρlLdYLtdt,y=YHt:−kl∂Tl∂y+ks∂Ts∂y=ρlLdYHtdt.The third boundary condition of heat conduction is still adopted for the upper and lower boundaries of the soldering area in the phase transition process (Equation (12)).The temperature distributionux of each zone calculated; the heat conduction problem can be resolved to obtain the temperature variation Ty,t of the soldering area. Then, the temperature variation of the center of the soldering area is Tb/2,t, namely, the temperature profile.


## 2.3. A Dynamic Thomas Algorithm

FDM and FEM are the two most popular approaches to solve differential equations. The advantage of the FDM is its simple principle, but it can only get the numerical solution in the regular region. Instead, the FEM is inclined to solve the irregular region problems [32]. Since the boundary of the mathematical model in our study is a regular rectangular region, the finite-difference method has been adopted.


### 2.3.1. Ignoring Phase Transformations

Firstly, the heat conduction equation of the soldering area was solved by the implicit difference method, according to Equation (11), as follows:(17)Ti,j−1=αTi−1,j+1−2αTi,j+αTi+1,j,where α=kΔt/cρΔy2=λΔt/Δy2, i=y/Δy=0,1,2,…,m, and j=t/Δt=0,1,2,…,n. The boundary conditions (Equation (12)) were discretized using backward difference for the upper and forward difference for the lower.(18)1+μTm,j−Tm−1,j=μux0+vt,1+μT0,j−T1,j=μux0+vt,where μ=kcΔy/k. The initial condition is Ti,0=25. Equation (17) is an implicit difference scheme for a heat conduction model of a soldering region without phase changes, which can be treated as a tridiagonal system of linear equations with respect to Ti,j. Substitute the boundary conditions, and the matrix form can be expressed as follows:(19)1+2α−α1+μ−α0⋯0−α1+2α−α⋯0⋮⋱⋱⋱⋮0⋯−α1+2α−α0⋯0−α1+2α−α1+μT1,jT2,j⋮Tm−2,jTm−1,j=T1,j−1+UjT2,j⋮Tm−2,j−1Tm−1,j−1+Uj,where Uj=αμ/1+μux0+vjΔt. The system of linear equations can be recorded as AX=B, and the nonzero elements in the coefficient matrix A are distributed on the leading diagonal and the two pairs of adjacent diagonals above and below. This form of linear equations can be solved quickly by the forward elimination and backward substitution method, namely, the Thomas algorithm [32].


### 2.3.2. Considering Phase Transitions

The heat conduction equation and boundary conditions are similar to the model without phase transitions, but the difference is that the initial condition should be the melting point of the solder. Consequently, the following are the equations for the melting process:(20)Ti,j−1l=αlTi−1,jl+1−2αlTi,jl+αlTi+1,jl,Ti,j−1s=αsTi−1,js+1−2αsTi,js+αsTi+1,js,Ti,0=Tmelting point,1+μlTm,j−Tm−1,j=μlux0+vt,1+μlT0,j−T1,j=μlux0+vt,where μl=kclΔy/kl, αl=klΔt/clρlΔy2=λlΔt/Δy2, and αs=ksΔt/csρsΔy2=λsΔt/Δy2. Tl is the temperature of the liquid phase, and Ts is the temperature of the solid phase. The cooling process is opposite, and equations are similar.(21)Ti,j−1s=αsTi−1,js+1−2αsTi,js+αsTi+1,js,Ti,j−1l=αlTi−1,jl+1−2αlTi,jl+αlTi+1,jl,Ti,0=Tfreezing point,1+μsTm,j−Tm−1,j=μsux0+vt,1+μsT0,j−T1,j=μsux0+vt,where μs=kcsΔy/kl. When phases changing, the initial temperature of the solid and the liquid phases are supposed to satisfy the heat transfer equation with the latent heat of phase transitions. For the boundaries (Equation (15)), the liquid phase adopted forward difference, and the solid phase adopted backward difference at the upper boundary, and it is just the reverse at the lower boundary. Hence, the interface conditions of phase transitions for the melting process are as follows:(22)−klTYL,j−TYL−1,j+ksTYL+1,j−TYL,j=ρsLvLΔy,klTYH+1,j−TYH,j−ksTYH,j−TYH−1,j=−ρsLvHΔy,where vL=ΔYL/Δt and vH=ΔYH/Δt. Similarly, the interface conditions during cooling are as follows:(23)klTYL,j−TYL−1,j−ksTYL+1,j−TYL,j=−ρlLvLΔy,−klTYH+1,j−TYH,j+ksTYH,j−TYH−1,j=ρlLvHΔy.For the heat conduction model of the phase transition process, the solid and the liquid regions are the same as the model without phase transitions, yet there are two more dynamic boundary conditions, so the soldering area should be divided into three sections for calculations as follows:(24)1+2α1−α11+μ−α10⋯0−α11+2α1−α1⋯0⋮⋱⋱⋱⋮0⋯−α11+2α1−α10⋯0−α11+2α1T1,jT2,j⋮TYL−2,jTYL−1,j=T1,j−1+UjT2,j⋮TYL−2,j−1TYL−1,j−1+αTYL,j,1+2α2−α20⋯0−α21+2α2−α2⋯0⋮⋱⋱⋱⋮0⋯−α21+2α2−α20⋯0−α21+2α2TYL+1,jTYL+2,j⋮TYH−2,jTYH−1,j=TYL+1,j−1+αTYL,jTYL+2,j⋮TYH−2,j−1TYH−1,j−1+αTYH,j,1+2α1−α10⋯0−α11+2α1−α⋯0⋮⋱⋱⋱⋮0⋯−α1+2α1−α10⋯0−α11+2α1−α11+μTYH+1,jTYH+2,j⋮Tm−2,jTm−1,j=TYH+1,j−1+αTYL,jTYH+2,j⋮Tm−2,j−1Tm−1,j−1+Uj,where μ=μl, α1=αl, α2=αs, and Uj=αlμl/1+μlux0+vjΔt when melting, and μ=μs, α1=αs, α2=αl, and Uj=αsμs/1+μsux0+vjΔt when cooling. Substitute the transformation boundary conditions (Equations (22) and (23)) that TYL,j=−ρlLvLΔy+klTL−1,j+ksTL+1,j/kl+ks and TYH,j=−ρlLvHΔy+klTH+1,j+ksTH−1,j/kl+ks during melting and TYL,j=−ρsLvLΔy+klTL−1,j+ksTL+1,j/kl+ks and TYH,j=−ρsLvHΔy+klTH+1,j+ksTH−1,j/kl+ks during cooling, and the following relation can be obtained:(25)1+2α1−α11+μ1−α1−α11+2α1−α1⋱⋱⋱−α11+2α1klks+kl−α1ksks+kl−α2klks+kl1+2α2ksks+kl−α2⋱⋱⋱−α21+2α2−α2⋱⋱⋱−α21+2α2ksks+kl−α2klks+kl−α1ksks+kl1+2α1klks+kl−α1⋱⋱⋱−α11+2α1−α1−α11+2α1−α11+μT1,jT2,j⋮TYL−1,jTYL+1,j⋮Ti,j⋮TYH−1,jTYH+1,j⋮Tm−2,jTm−1,j=T1,j−1+UjT2,j⋮TYL−1,j−1−LjTYL+1,j−1−Lj⋮Ti,j−1⋮TYH−1,j−1−HjTYH+1,j−1−Hj⋮Tm−2,j−1Tm−1,j−1+Uj,where Lj=αlρlLvLΔy/kskl and Hj=αlρlLvHΔy/kskl during the melting process, and Lj=αlρlLvLΔy/kskl and Hj=αlρlLvHΔy/kskl during the cooling process. The coefficient matrix of equations without phase transformations is independent of time and can be solved quickly by the Thomas algorithm. As for the phase transition process, the coefficient matrix changes with time because interfaces of the solid and the liquid phases approach gradually. Consequently, Equation (25) can be expressed as AtXt=Yt. At, Xt and Yt can be written in the following form:(26)At=b1tc1t0⋯0a2tb2tc2t⋯0⋮⋱⋱⋱⋮0⋯am−4tbm−4tcm−4t0⋯0am−3tbm−3t,Xt=x1tx2t⋮xm−4txm−3t,Yt=y1ty2t⋮ym−4tym−3t.Although the coefficient matrix changes with time, the system still satisfies the condition of the Crout decomposition method, that is, (i)b1t≥c1t≥0, (ii) bit≥ait+cit, aitcit≠0, i=2,3,…,m−3, and (iii) bm−3t≥am−3t≥0. Obviously, At is a nonsingular matrix whose determinant detAt≠0. Thus, its sequential principal minor of any order Dkt≠0.A tridiagonal square matrix satisfying the above conditions has a unique Crout decomposition matrix as follows:(27)At=PtQt=p1t00⋯0a2tp2t0⋯0⋮⋱⋱⋱⋮0⋯am−4tpm−4t00⋯0am−3tpm−3t1q1t0⋯001q2t⋯0⋮⋱⋱⋱⋮0⋯01qm−4t0⋯001,where p1t=b1t, qit=cit/piti=1,2,…,m−4, and pit=bit−aitqi−1ti=2,3,…,m−3. pit and qit are time-varying and piecewise functions, and the recurrence formulas can be expressed as follows:(28)qit=−α11+2α1−α1/1+μ,i=1;−α11+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;−α1ks/ks+kl1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;−α21+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;−α21+2α2+α2qi−1t,YLt≤i≤YHt−3;−α2kl/ks+kl1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;−α11+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1.pit=1+2α1−α11+μ,i=1;1+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;1+2α1klks+kl+α1qi−1t,i=YLt−1;1+2α2ksks+kl+α2klks+klqi−1t,i=YLt;1+2α2+α2qi−1t,YLt≤i≤YHt−3;1+2α2kSks+kl+α2qi−1t,i=YHt−2;1+2α2klks+kl+α1ksks+klqi−1t,i=YHt−1;1+2α1−α11+μ+α1qi−1t,i=m−3.Therefore,AtXt=PtQtXt=Yt can be decomposed into two triangular equations PtDt=Yt and QtXt=Dt. The first step is to solve the system of equations PtDt=Yt as follows:(29)p1t00⋯0a2tp2t0⋯0⋮⋱⋱⋱⋮0⋯am−4tpm−4t00⋯0am−3tpm−3td1td2t⋮dm−4tdm−3t=y1ty2t⋮ym−4tym−3t.The recurrence formulas of solutions are as follows:(30)dit=y1t1+2α1−α1/1+μ,i=1;yit+α1di−1t1+2α1+α1qi−1t2≤i≤YLt−2orYHt≤i≤m−4;yit+α1di−1t1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;yit+α2kl/ks+kldi−1t1+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;yit+α2di−1t1+2α2+α2qi−1t,YLt≤i≤YHt−3;yit+α2di−1t1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;yit+α1ks/ks+kldi−1t1+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1;yit+α1di−1t1+2α1−α1/1+μ+α1qi−1t,i=m−3.The second step is to solve the system of equationsQtXt=Dt as follows:(31)1q1t0⋯001q2t⋯0⋮⋱⋱⋱⋮0⋯01qm−4t0⋯001x1tx2t⋮xm−4txm−3t=d1td2t⋮dm−4tdm−3t.The recurrence formulas of solutions are as follows:(32)xit=yit+α1di−1t1+2α1−α1/1+μ+α1qi−1t,i=m−3;yit+α1ks/ks+kldi−1t+α1xi+1t1+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1;yit+α2di−1t+α1kl/ks+klxi+1t1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;yit+α2di−1t+α2xi+1t1+2α2+α2qi−1t,YLt≤i≤YHt−3;yit+α2kl/ks+kldi−1t+α2xi+1t1+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;yit+α1di−1t+α1ks/ks+klxi+1t1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;yit+α1di−1t+α1xi+1t1+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;y1t+α1xi+1t1+2α1−α1/1+μ,i=1.wherexit is the discrete solution of the heat transfer equation in the case of considering phase transitions. A dynamic Thomas algorithm is proposed and developed in this study, in order to address the time-dependent solution of the tridiagonal linear system. Flowchart of the developed algorithm is indicated in Figure 5. By cyclic iteration, the temperature of the soldering area during the entire process can be obtained.Figure 5
Flowchart of the dynamic Thomas algorithm, that is, forward elimination and backward substitution with time.


## 2.3.1. Ignoring Phase Transformations

Firstly, the heat conduction equation of the soldering area was solved by the implicit difference method, according to Equation (11), as follows:(17)Ti,j−1=αTi−1,j+1−2αTi,j+αTi+1,j,where α=kΔt/cρΔy2=λΔt/Δy2, i=y/Δy=0,1,2,…,m, and j=t/Δt=0,1,2,…,n. The boundary conditions (Equation (12)) were discretized using backward difference for the upper and forward difference for the lower.(18)1+μTm,j−Tm−1,j=μux0+vt,1+μT0,j−T1,j=μux0+vt,where μ=kcΔy/k. The initial condition is Ti,0=25. Equation (17) is an implicit difference scheme for a heat conduction model of a soldering region without phase changes, which can be treated as a tridiagonal system of linear equations with respect to Ti,j. Substitute the boundary conditions, and the matrix form can be expressed as follows:(19)1+2α−α1+μ−α0⋯0−α1+2α−α⋯0⋮⋱⋱⋱⋮0⋯−α1+2α−α0⋯0−α1+2α−α1+μT1,jT2,j⋮Tm−2,jTm−1,j=T1,j−1+UjT2,j⋮Tm−2,j−1Tm−1,j−1+Uj,where Uj=αμ/1+μux0+vjΔt. The system of linear equations can be recorded as AX=B, and the nonzero elements in the coefficient matrix A are distributed on the leading diagonal and the two pairs of adjacent diagonals above and below. This form of linear equations can be solved quickly by the forward elimination and backward substitution method, namely, the Thomas algorithm [32].


## 2.3.2. Considering Phase Transitions

The heat conduction equation and boundary conditions are similar to the model without phase transitions, but the difference is that the initial condition should be the melting point of the solder. Consequently, the following are the equations for the melting process:(20)Ti,j−1l=αlTi−1,jl+1−2αlTi,jl+αlTi+1,jl,Ti,j−1s=αsTi−1,js+1−2αsTi,js+αsTi+1,js,Ti,0=Tmelting point,1+μlTm,j−Tm−1,j=μlux0+vt,1+μlT0,j−T1,j=μlux0+vt,where μl=kclΔy/kl, αl=klΔt/clρlΔy2=λlΔt/Δy2, and αs=ksΔt/csρsΔy2=λsΔt/Δy2. Tl is the temperature of the liquid phase, and Ts is the temperature of the solid phase. The cooling process is opposite, and equations are similar.(21)Ti,j−1s=αsTi−1,js+1−2αsTi,js+αsTi+1,js,Ti,j−1l=αlTi−1,jl+1−2αlTi,jl+αlTi+1,jl,Ti,0=Tfreezing point,1+μsTm,j−Tm−1,j=μsux0+vt,1+μsT0,j−T1,j=μsux0+vt,where μs=kcsΔy/kl. When phases changing, the initial temperature of the solid and the liquid phases are supposed to satisfy the heat transfer equation with the latent heat of phase transitions. For the boundaries (Equation (15)), the liquid phase adopted forward difference, and the solid phase adopted backward difference at the upper boundary, and it is just the reverse at the lower boundary. Hence, the interface conditions of phase transitions for the melting process are as follows:(22)−klTYL,j−TYL−1,j+ksTYL+1,j−TYL,j=ρsLvLΔy,klTYH+1,j−TYH,j−ksTYH,j−TYH−1,j=−ρsLvHΔy,where vL=ΔYL/Δt and vH=ΔYH/Δt. Similarly, the interface conditions during cooling are as follows:(23)klTYL,j−TYL−1,j−ksTYL+1,j−TYL,j=−ρlLvLΔy,−klTYH+1,j−TYH,j+ksTYH,j−TYH−1,j=ρlLvHΔy.For the heat conduction model of the phase transition process, the solid and the liquid regions are the same as the model without phase transitions, yet there are two more dynamic boundary conditions, so the soldering area should be divided into three sections for calculations as follows:(24)1+2α1−α11+μ−α10⋯0−α11+2α1−α1⋯0⋮⋱⋱⋱⋮0⋯−α11+2α1−α10⋯0−α11+2α1T1,jT2,j⋮TYL−2,jTYL−1,j=T1,j−1+UjT2,j⋮TYL−2,j−1TYL−1,j−1+αTYL,j,1+2α2−α20⋯0−α21+2α2−α2⋯0⋮⋱⋱⋱⋮0⋯−α21+2α2−α20⋯0−α21+2α2TYL+1,jTYL+2,j⋮TYH−2,jTYH−1,j=TYL+1,j−1+αTYL,jTYL+2,j⋮TYH−2,j−1TYH−1,j−1+αTYH,j,1+2α1−α10⋯0−α11+2α1−α⋯0⋮⋱⋱⋱⋮0⋯−α1+2α1−α10⋯0−α11+2α1−α11+μTYH+1,jTYH+2,j⋮Tm−2,jTm−1,j=TYH+1,j−1+αTYL,jTYH+2,j⋮Tm−2,j−1Tm−1,j−1+Uj,where μ=μl, α1=αl, α2=αs, and Uj=αlμl/1+μlux0+vjΔt when melting, and μ=μs, α1=αs, α2=αl, and Uj=αsμs/1+μsux0+vjΔt when cooling. Substitute the transformation boundary conditions (Equations (22) and (23)) that TYL,j=−ρlLvLΔy+klTL−1,j+ksTL+1,j/kl+ks and TYH,j=−ρlLvHΔy+klTH+1,j+ksTH−1,j/kl+ks during melting and TYL,j=−ρsLvLΔy+klTL−1,j+ksTL+1,j/kl+ks and TYH,j=−ρsLvHΔy+klTH+1,j+ksTH−1,j/kl+ks during cooling, and the following relation can be obtained:(25)1+2α1−α11+μ1−α1−α11+2α1−α1⋱⋱⋱−α11+2α1klks+kl−α1ksks+kl−α2klks+kl1+2α2ksks+kl−α2⋱⋱⋱−α21+2α2−α2⋱⋱⋱−α21+2α2ksks+kl−α2klks+kl−α1ksks+kl1+2α1klks+kl−α1⋱⋱⋱−α11+2α1−α1−α11+2α1−α11+μT1,jT2,j⋮TYL−1,jTYL+1,j⋮Ti,j⋮TYH−1,jTYH+1,j⋮Tm−2,jTm−1,j=T1,j−1+UjT2,j⋮TYL−1,j−1−LjTYL+1,j−1−Lj⋮Ti,j−1⋮TYH−1,j−1−HjTYH+1,j−1−Hj⋮Tm−2,j−1Tm−1,j−1+Uj,where Lj=αlρlLvLΔy/kskl and Hj=αlρlLvHΔy/kskl during the melting process, and Lj=αlρlLvLΔy/kskl and Hj=αlρlLvHΔy/kskl during the cooling process. The coefficient matrix of equations without phase transformations is independent of time and can be solved quickly by the Thomas algorithm. As for the phase transition process, the coefficient matrix changes with time because interfaces of the solid and the liquid phases approach gradually. Consequently, Equation (25) can be expressed as AtXt=Yt. At, Xt and Yt can be written in the following form:(26)At=b1tc1t0⋯0a2tb2tc2t⋯0⋮⋱⋱⋱⋮0⋯am−4tbm−4tcm−4t0⋯0am−3tbm−3t,Xt=x1tx2t⋮xm−4txm−3t,Yt=y1ty2t⋮ym−4tym−3t.Although the coefficient matrix changes with time, the system still satisfies the condition of the Crout decomposition method, that is, (i)b1t≥c1t≥0, (ii) bit≥ait+cit, aitcit≠0, i=2,3,…,m−3, and (iii) bm−3t≥am−3t≥0. Obviously, At is a nonsingular matrix whose determinant detAt≠0. Thus, its sequential principal minor of any order Dkt≠0.A tridiagonal square matrix satisfying the above conditions has a unique Crout decomposition matrix as follows:(27)At=PtQt=p1t00⋯0a2tp2t0⋯0⋮⋱⋱⋱⋮0⋯am−4tpm−4t00⋯0am−3tpm−3t1q1t0⋯001q2t⋯0⋮⋱⋱⋱⋮0⋯01qm−4t0⋯001,where p1t=b1t, qit=cit/piti=1,2,…,m−4, and pit=bit−aitqi−1ti=2,3,…,m−3. pit and qit are time-varying and piecewise functions, and the recurrence formulas can be expressed as follows:(28)qit=−α11+2α1−α1/1+μ,i=1;−α11+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;−α1ks/ks+kl1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;−α21+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;−α21+2α2+α2qi−1t,YLt≤i≤YHt−3;−α2kl/ks+kl1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;−α11+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1.pit=1+2α1−α11+μ,i=1;1+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;1+2α1klks+kl+α1qi−1t,i=YLt−1;1+2α2ksks+kl+α2klks+klqi−1t,i=YLt;1+2α2+α2qi−1t,YLt≤i≤YHt−3;1+2α2kSks+kl+α2qi−1t,i=YHt−2;1+2α2klks+kl+α1ksks+klqi−1t,i=YHt−1;1+2α1−α11+μ+α1qi−1t,i=m−3.Therefore,AtXt=PtQtXt=Yt can be decomposed into two triangular equations PtDt=Yt and QtXt=Dt. The first step is to solve the system of equations PtDt=Yt as follows:(29)p1t00⋯0a2tp2t0⋯0⋮⋱⋱⋱⋮0⋯am−4tpm−4t00⋯0am−3tpm−3td1td2t⋮dm−4tdm−3t=y1ty2t⋮ym−4tym−3t.The recurrence formulas of solutions are as follows:(30)dit=y1t1+2α1−α1/1+μ,i=1;yit+α1di−1t1+2α1+α1qi−1t2≤i≤YLt−2orYHt≤i≤m−4;yit+α1di−1t1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;yit+α2kl/ks+kldi−1t1+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;yit+α2di−1t1+2α2+α2qi−1t,YLt≤i≤YHt−3;yit+α2di−1t1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;yit+α1ks/ks+kldi−1t1+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1;yit+α1di−1t1+2α1−α1/1+μ+α1qi−1t,i=m−3.The second step is to solve the system of equationsQtXt=Dt as follows:(31)1q1t0⋯001q2t⋯0⋮⋱⋱⋱⋮0⋯01qm−4t0⋯001x1tx2t⋮xm−4txm−3t=d1td2t⋮dm−4tdm−3t.The recurrence formulas of solutions are as follows:(32)xit=yit+α1di−1t1+2α1−α1/1+μ+α1qi−1t,i=m−3;yit+α1ks/ks+kldi−1t+α1xi+1t1+2α2kl/ks+kl+α1ks/ks+klqi−1t,i=YHt−1;yit+α2di−1t+α1kl/ks+klxi+1t1+2α2kS/ks+kl+α2qi−1t,i=YHt−2;yit+α2di−1t+α2xi+1t1+2α2+α2qi−1t,YLt≤i≤YHt−3;yit+α2kl/ks+kldi−1t+α2xi+1t1+2α2ks/ks+kl+α2kl/ks+klqi−1t,i=YLt;yit+α1di−1t+α1ks/ks+klxi+1t1+2α1kl/ks+kl+α1qi−1t,i=YLt−1;yit+α1di−1t+α1xi+1t1+2α1+α1qi−1t,2≤i≤YLt−2orYHt≤i≤m−4;y1t+α1xi+1t1+2α1−α1/1+μ,i=1.wherexit is the discrete solution of the heat transfer equation in the case of considering phase transitions. A dynamic Thomas algorithm is proposed and developed in this study, in order to address the time-dependent solution of the tridiagonal linear system. Flowchart of the developed algorithm is indicated in Figure 5. By cyclic iteration, the temperature of the soldering area during the entire process can be obtained.Figure 5
Flowchart of the dynamic Thomas algorithm, that is, forward elimination and backward substitution with time.


## 3. Results of Numerical Simulation

Numerical simulations are carried out for investigating the accuracy of the mathematical model and the validity of the dynamic Thomas algorithm. The temperature profiles are obtained by the calculations with numerical simulations, while the best experimental parameters are determined by a multiobjective optimization model.


### 3.1. Comparison of Temperature Profiles

The computing method and results are given by taking the reflow oven with 11 small temperature zones as an example. In an experiment, the set temperature in each zone is presented in Table1.Table 1
Set temperature of each zone.
Area1–5678 and 910 and 11Temperature175°C195°C235°C255°C25°CTemperature profiles and their relative errors are calculated, as shown in Figures6(a) and 6(b), respectively. Figure 6(a) shows the temperature profiles calculated for considering phase changes, ignoring phase transitions, and experimental measurement. Figure 6(b) indicates the relative errors between the numerical results and experimental values.Figure 6
Results of numerical simulations: (a) temperature profiles and (b) relative errors.
(a)(b)In the region marked byPa of Figure 6(a), one can see that the deviation of the temperature-profile curve obtained through considering phase transitions from the experimental values is larger than the result of ignoring phase transitions, which can also be seen from the error curve of Figure 6(b) in the time range 240s≤t≤300s. The reason for this result could be that the solder density has changed when the phase transition of the solder happens in the high-temperature zone (Equation (24)). In the region marked by Pb of Figure 6(a), it is obvious that the freezing rate of the red curve is smaller than that of the yellow one, which reveals that the heat is absorbed when melting and released during solidification, leading to the temperature change more slowly. As shown in Figure 6(b), the values of the blue curve, which represents the relative errors between the theoretical and the experimental values for the case of ignoring phase transitions, are much larger than that of the purple curve when t≥300s. Figure 6(b) indicates that the relative error mostly fluctuates around 1.0%. However, the error becomes larger for t≤60s, even up to 10%, which is due to the fact that the actual temperature is low for t≤60s. As expected, phase transitions have little effect on the section where the temperature is lower than the melting point. As shown in Figures 6(a) and 6(b), the two curves almost coincide when t<175s.The overall temperature distributions of the soldering area can also be calculated, for both considering and ignoring phase transitions, as shown in FiguresS1 and S2 in the Supplementary Materials, respectively. Besides, in order to highlight the difference, Figure S3 in the Supplementary Materials shows the results that the temperature distribution in the case of phase transitions subtracts the other. In general, either the temperature profile or the temperature distribution of the soldering area calculated with phase changes is preferable.


### 3.2. Related Parameters Optimization

After verifying the validity of the mathematical model of the soldering furnace, the relevant experimental parameters that affect the temperature profile can be optimized, which can decrease the cost and increase the productivity. Actually, the conveyor speed and the temperature profile are supposed to meet certain requirements, for the purpose of ensuring the quality of the soldering. With the limitations indicated in Table2, the lasting time that the temperature of the center of soldering area exceeds 217∘C, Δt=t217,R−t217,L, should not be too long, and the peak temperature Tpeak should not be too high. We now define the integration: Qη=∫t217,Ltpeakft,vp−217dt, where ft,vp is the temperature function plotted in Figure 6. Note that Qη should reach a minimum value for an ideal temperature profile.Table 2
Process limitations.
ProcessMinimumMaximumSlope−3°C/s3°C/sHeating time when150∘C<T<190∘C60 s120 sThe time whenT>217∘C40 s90 sTemperature peak240∘C250∘CMultiobjective optimization is the feature of this problem, which contains three objective functions, five decision variables, and four constraint conditions. Such a problem usually needs a huge amount of computation and yields a set of solutions that are incomparable in terms of the overall objective functions. This kind of solution is called the nondominated solution or the Pareto optimal solution, whose characteristic is that it is impossible to improve any objective functions without weakening at least one of the other objective functions [33]. For all Pareto optimal solutions that constitute the Pareto optimal solution set, the surface, formed by the target vector in the corresponding target space, is called the Pareto optimality, which is the key to solve the multiobjective optimization problem [34]. An improved genetic algorithm, nondominated sorting genetic algorithm II (NSGA-II), was developed to solve the problem [35–37]. Through numerical experiments where the number of individuals and iterations are 80 and 500, respectively, the scatter diagrams of objective functions are shown in Figure 7, and the specific results of 10 individuals, which were calculated in the case of ignoring and considering phase transitions are shown in the Supplementary Materials, where Tables S1 and S2 show the values of three objective functions, and Tables S3 and S4 show the values of five decision variables.Figure 7
Numerical results of optimization of three objective functions: (a) ignoring phase transitions, (b) considering phase changes, and (c) relative position of results in the two cases.
(a)(b)(c)The results reveal that the values of the objective functions in the case of considering phase transitions are more concentrated than the other (Figure7(c)), although its optimization results of obj1. Qη and obj3. Δt are higher than that under the condition of ignoring phase transformations (Tables S1 and S2 in the Supplementary Materials). There is a similarity worth noting that the results of obj2. Peak all tend to the minimum 240∘C in the two cases, as can be seen in Tables S1 and S2 in the Supplementary Materials. For decision variables, the set temperature of zones 1∼5, zone 6, and zone 7 and the belt speed are lower than the other, comparing Tables S3 and S4. Thus, the conclusion is that either 3 objective functions or 5 decision variables in the two cases are dissimilar enough to influence the actual production.The optimal temperature profiles are shown in Figure8, where both the objective function Qη and decision variable v vary greatly. Therefore, the time taken for the entire soldering process is different too. The red curve that is the optimal temperature profile taking phase transformations into account ends at 275 s approximately, but the blue one, calculated in the case of leaving phase changes out, ends at about 290 s. As expected, their peaks are the same. In short, phase transitions have a significant influence on the optimization of the soldering temperature profile.Figure 8
The optimal temperature profiles in the two cases, considering and ignoring phase changes.


### 3.3. Verification of the Linear Relationship betweenQη′min and Ht

Qη′=∫t217,Lt217,Rft,vp−217dt is the heating factor; Qη′min is optimal; and Ht stands for the maximum set temperature. Results of the multiobjective optimization, shown in Figure 8, are taken as a reference. For taking no account of phase changes, the results are u8=u9=259.68∘C, u7=227.20∘C, u6=202.32∘C, u1∼5=179.32∘C, v=93.89cm/min, and Qη′min=573.03∘C⋅s, and with considering phase transitions, we have u8=u9=259.51∘C, u7=232.78∘C, u6=190.55∘C, u1∼5=173.37∘C, v=90.39cm/min, and Qη′min=652.87∘C⋅s. Then, we establish the functional relationships according to the reference that u8=u9=Ht, u7=H7−32.48, u6=Ht−57.36, and u1∼5=Ht−80.36 for ignoring phase transformations, and u8=u9=Ht, u7=H7−26.73, u6=Ht−68.96, and u1∼5=Ht−86.14 for considering phase transitions. When we also make small modifications to Ht, ΔHt≤10∘C, and keep the other parameters unchanged, the relationship between Qη′min and Ht can be obtained, as shown in Figure 9.Figure 9
The relationship between the optimal heating factorQη′min and the maximum set temperature Ht in the range of ΔHt≤10∘C.Obviously,Qη′min is linearly related to Ht in the two cases for the segments of BaBc and RaRc, respectively, as shown in Figure 9 when ΔHt≤5∘C. This result is consistent with the conclusion of previous work [4]. It needs to be emphasized that the conclusion, the linear relationship between Qη′min and Ht, was previously obtained from a large number of experiments [4], while we have had the same conclusion based on our mathematical model of heat transfer and multiobjective optimization. Moreover, by comparing the curve of Qηb′min and Qηr′min, it can be found that although the variation trends of them tend to be nonlinear outside the range ΔHt≤5∘C, the difference is that Qη′minb changes continuously, while Qηr′min has an obvious discontinuity point when the temperature is about 252∘C or 269∘C. Therefore, it can be concluded that phase transitions break up Qη′min in some places, which is supposed to change continuously, and form a banded distribution. Hence, we further expand the range of Ht, ΔHt≤15∘C, and take more dense points to show the banded distribution of Qηr′min in Figure 10, calculated by the model of involving phase transformations.Figure 10
The banded distribution of the heating factorQη′min in the range of ΔHt≤15∘C.Allowed bands are the places whereQη′min can be continuously sampled, and forbidden bands are the areas where Qηr′min is discontinuous. In general, Qη′min can always be determined by changing the value of Ht, so Qηb′min is continuous disregarding the phase transitions. A new phenomenon has been found that under the condition of phase changes, Qηr′min cannot sample all the values continuously except the points in the range of allowed bands.Based on our analysis, the discontinuity of the heating factorQηr′min should be related to the constraint conditions of the process limitations, because the value of Ht is further expanded on the basis of the results of the constrained multiobjective optimization. The influence of each constraint on Qηr′min can be revealed separately by calculating the range of Ht allowed by process limitations shown in Table 2.From Figures11(b) and 11(c), it can be found that the first jump discontinuity, where Ht is about 252∘C, is caused by the violation of the restriction that the time when T>217∘C on the temperature profile must be 40s∼90s, and the second jump discontinuity, where Ht is about 269∘C, is resulted from going against two limitations that the rate of temperature change must be −3∘C/s∼3∘C/s, and the heating time must be 60s<T<120s when 150∘C<T<190∘C. However, there is no connection between the other two constraints and the phenomenon that there are discontinuities in function Qηr′min=FHt, as shown in Figures 11(a) and 11(d).Figure 11
Effects of four process constraints in the range ofΔHt≤15∘C: (a) peak limiting, (b) slope limiting, (c) time limiting when T>217∘C, and (d) heating time limiting when 150∘C<T<190∘C.
(a)(b)(c)(d)In summary, feedback of the heating factorQη′min to process limitations with or without phase changes is completely different, especially the interaction between phase transformations, the heating factor, and slope of temperature profiles. However, the main reason that this phenomenon only occurs when considering the phase transition may include the following views. There are two interfaces between the solid and the liquid phases in the soldering area when considering phase changes. The temperature on the interface, like the boundary, cannot be calculated by the heat conduction equation but can only be figured out by the boundary condition (Equations (22) and (23)). This is why the matrix in the dynamic Thomas algorithm is m−3 rows, while the original has m−1 rows (Equations (26) and (19)). Therefore, when considering phase transitions, the temperature calculated by the dynamic Thomas algorithm is broken at the two interfaces, because TYL,j and TYH,j are removed (Equation (25)). This explanation is difficult to be verified directly in theory, but it can be inferred that if the addition of two interfaces during phase changes does result in the discontinuities in the heating factor, the number of jump discontinuities must be consistent with the number of interfaces. That is the function Qηr′min=FHt that has only two discontinuities, regardless of whether the process limitations are met. Let’s expand the range of Ht even further, as shown in Figure 12.Figure 12
Effects of four process constraints in the range ofΔHt≤60∘C: (a) peak limiting, (b) slope limiting, (c) time limiting when T>217∘C, and (d) heating time limiting when 150∘C<T<190∘C.
(a)(b)(c)(d)Obviously, this inference is correct. It can be concluded that the jump discontinuities appear due to the increase of two interfaces between the solid and the liquid phases during phase transitions and the process limitations, but the number of discontinuity points is only related to that of interfaces between two phases.


## 3.1. Comparison of Temperature Profiles

The computing method and results are given by taking the reflow oven with 11 small temperature zones as an example. In an experiment, the set temperature in each zone is presented in Table1.Table 1
Set temperature of each zone.
Area1–5678 and 910 and 11Temperature175°C195°C235°C255°C25°CTemperature profiles and their relative errors are calculated, as shown in Figures6(a) and 6(b), respectively. Figure 6(a) shows the temperature profiles calculated for considering phase changes, ignoring phase transitions, and experimental measurement. Figure 6(b) indicates the relative errors between the numerical results and experimental values.Figure 6
Results of numerical simulations: (a) temperature profiles and (b) relative errors.
(a)(b)In the region marked byPa of Figure 6(a), one can see that the deviation of the temperature-profile curve obtained through considering phase transitions from the experimental values is larger than the result of ignoring phase transitions, which can also be seen from the error curve of Figure 6(b) in the time range 240s≤t≤300s. The reason for this result could be that the solder density has changed when the phase transition of the solder happens in the high-temperature zone (Equation (24)). In the region marked by Pb of Figure 6(a), it is obvious that the freezing rate of the red curve is smaller than that of the yellow one, which reveals that the heat is absorbed when melting and released during solidification, leading to the temperature change more slowly. As shown in Figure 6(b), the values of the blue curve, which represents the relative errors between the theoretical and the experimental values for the case of ignoring phase transitions, are much larger than that of the purple curve when t≥300s. Figure 6(b) indicates that the relative error mostly fluctuates around 1.0%. However, the error becomes larger for t≤60s, even up to 10%, which is due to the fact that the actual temperature is low for t≤60s. As expected, phase transitions have little effect on the section where the temperature is lower than the melting point. As shown in Figures 6(a) and 6(b), the two curves almost coincide when t<175s.The overall temperature distributions of the soldering area can also be calculated, for both considering and ignoring phase transitions, as shown in FiguresS1 and S2 in the Supplementary Materials, respectively. Besides, in order to highlight the difference, Figure S3 in the Supplementary Materials shows the results that the temperature distribution in the case of phase transitions subtracts the other. In general, either the temperature profile or the temperature distribution of the soldering area calculated with phase changes is preferable.


## 3.2. Related Parameters Optimization

After verifying the validity of the mathematical model of the soldering furnace, the relevant experimental parameters that affect the temperature profile can be optimized, which can decrease the cost and increase the productivity. Actually, the conveyor speed and the temperature profile are supposed to meet certain requirements, for the purpose of ensuring the quality of the soldering. With the limitations indicated in Table2, the lasting time that the temperature of the center of soldering area exceeds 217∘C, Δt=t217,R−t217,L, should not be too long, and the peak temperature Tpeak should not be too high. We now define the integration: Qη=∫t217,Ltpeakft,vp−217dt, where ft,vp is the temperature function plotted in Figure 6. Note that Qη should reach a minimum value for an ideal temperature profile.Table 2
Process limitations.
ProcessMinimumMaximumSlope−3°C/s3°C/sHeating time when150∘C<T<190∘C60 s120 sThe time whenT>217∘C40 s90 sTemperature peak240∘C250∘CMultiobjective optimization is the feature of this problem, which contains three objective functions, five decision variables, and four constraint conditions. Such a problem usually needs a huge amount of computation and yields a set of solutions that are incomparable in terms of the overall objective functions. This kind of solution is called the nondominated solution or the Pareto optimal solution, whose characteristic is that it is impossible to improve any objective functions without weakening at least one of the other objective functions [33]. For all Pareto optimal solutions that constitute the Pareto optimal solution set, the surface, formed by the target vector in the corresponding target space, is called the Pareto optimality, which is the key to solve the multiobjective optimization problem [34]. An improved genetic algorithm, nondominated sorting genetic algorithm II (NSGA-II), was developed to solve the problem [35–37]. Through numerical experiments where the number of individuals and iterations are 80 and 500, respectively, the scatter diagrams of objective functions are shown in Figure 7, and the specific results of 10 individuals, which were calculated in the case of ignoring and considering phase transitions are shown in the Supplementary Materials, where Tables S1 and S2 show the values of three objective functions, and Tables S3 and S4 show the values of five decision variables.Figure 7
Numerical results of optimization of three objective functions: (a) ignoring phase transitions, (b) considering phase changes, and (c) relative position of results in the two cases.
(a)(b)(c)The results reveal that the values of the objective functions in the case of considering phase transitions are more concentrated than the other (Figure7(c)), although its optimization results of obj1. Qη and obj3. Δt are higher than that under the condition of ignoring phase transformations (Tables S1 and S2 in the Supplementary Materials). There is a similarity worth noting that the results of obj2. Peak all tend to the minimum 240∘C in the two cases, as can be seen in Tables S1 and S2 in the Supplementary Materials. For decision variables, the set temperature of zones 1∼5, zone 6, and zone 7 and the belt speed are lower than the other, comparing Tables S3 and S4. Thus, the conclusion is that either 3 objective functions or 5 decision variables in the two cases are dissimilar enough to influence the actual production.The optimal temperature profiles are shown in Figure8, where both the objective function Qη and decision variable v vary greatly. Therefore, the time taken for the entire soldering process is different too. The red curve that is the optimal temperature profile taking phase transformations into account ends at 275 s approximately, but the blue one, calculated in the case of leaving phase changes out, ends at about 290 s. As expected, their peaks are the same. In short, phase transitions have a significant influence on the optimization of the soldering temperature profile.Figure 8
The optimal temperature profiles in the two cases, considering and ignoring phase changes.


## 3.3. Verification of the Linear Relationship betweenQη′min and Ht

Qη′=∫t217,Lt217,Rft,vp−217dt is the heating factor; Qη′min is optimal; and Ht stands for the maximum set temperature. Results of the multiobjective optimization, shown in Figure 8, are taken as a reference. For taking no account of phase changes, the results are u8=u9=259.68∘C, u7=227.20∘C, u6=202.32∘C, u1∼5=179.32∘C, v=93.89cm/min, and Qη′min=573.03∘C⋅s, and with considering phase transitions, we have u8=u9=259.51∘C, u7=232.78∘C, u6=190.55∘C, u1∼5=173.37∘C, v=90.39cm/min, and Qη′min=652.87∘C⋅s. Then, we establish the functional relationships according to the reference that u8=u9=Ht, u7=H7−32.48, u6=Ht−57.36, and u1∼5=Ht−80.36 for ignoring phase transformations, and u8=u9=Ht, u7=H7−26.73, u6=Ht−68.96, and u1∼5=Ht−86.14 for considering phase transitions. When we also make small modifications to Ht, ΔHt≤10∘C, and keep the other parameters unchanged, the relationship between Qη′min and Ht can be obtained, as shown in Figure 9.Figure 9
The relationship between the optimal heating factorQη′min and the maximum set temperature Ht in the range of ΔHt≤10∘C.Obviously,Qη′min is linearly related to Ht in the two cases for the segments of BaBc and RaRc, respectively, as shown in Figure 9 when ΔHt≤5∘C. This result is consistent with the conclusion of previous work [4]. It needs to be emphasized that the conclusion, the linear relationship between Qη′min and Ht, was previously obtained from a large number of experiments [4], while we have had the same conclusion based on our mathematical model of heat transfer and multiobjective optimization. Moreover, by comparing the curve of Qηb′min and Qηr′min, it can be found that although the variation trends of them tend to be nonlinear outside the range ΔHt≤5∘C, the difference is that Qη′minb changes continuously, while Qηr′min has an obvious discontinuity point when the temperature is about 252∘C or 269∘C. Therefore, it can be concluded that phase transitions break up Qη′min in some places, which is supposed to change continuously, and form a banded distribution. Hence, we further expand the range of Ht, ΔHt≤15∘C, and take more dense points to show the banded distribution of Qηr′min in Figure 10, calculated by the model of involving phase transformations.Figure 10
The banded distribution of the heating factorQη′min in the range of ΔHt≤15∘C.Allowed bands are the places whereQη′min can be continuously sampled, and forbidden bands are the areas where Qηr′min is discontinuous. In general, Qη′min can always be determined by changing the value of Ht, so Qηb′min is continuous disregarding the phase transitions. A new phenomenon has been found that under the condition of phase changes, Qηr′min cannot sample all the values continuously except the points in the range of allowed bands.Based on our analysis, the discontinuity of the heating factorQηr′min should be related to the constraint conditions of the process limitations, because the value of Ht is further expanded on the basis of the results of the constrained multiobjective optimization. The influence of each constraint on Qηr′min can be revealed separately by calculating the range of Ht allowed by process limitations shown in Table 2.From Figures11(b) and 11(c), it can be found that the first jump discontinuity, where Ht is about 252∘C, is caused by the violation of the restriction that the time when T>217∘C on the temperature profile must be 40s∼90s, and the second jump discontinuity, where Ht is about 269∘C, is resulted from going against two limitations that the rate of temperature change must be −3∘C/s∼3∘C/s, and the heating time must be 60s<T<120s when 150∘C<T<190∘C. However, there is no connection between the other two constraints and the phenomenon that there are discontinuities in function Qηr′min=FHt, as shown in Figures 11(a) and 11(d).Figure 11
Effects of four process constraints in the range ofΔHt≤15∘C: (a) peak limiting, (b) slope limiting, (c) time limiting when T>217∘C, and (d) heating time limiting when 150∘C<T<190∘C.
(a)(b)(c)(d)In summary, feedback of the heating factorQη′min to process limitations with or without phase changes is completely different, especially the interaction between phase transformations, the heating factor, and slope of temperature profiles. However, the main reason that this phenomenon only occurs when considering the phase transition may include the following views. There are two interfaces between the solid and the liquid phases in the soldering area when considering phase changes. The temperature on the interface, like the boundary, cannot be calculated by the heat conduction equation but can only be figured out by the boundary condition (Equations (22) and (23)). This is why the matrix in the dynamic Thomas algorithm is m−3 rows, while the original has m−1 rows (Equations (26) and (19)). Therefore, when considering phase transitions, the temperature calculated by the dynamic Thomas algorithm is broken at the two interfaces, because TYL,j and TYH,j are removed (Equation (25)). This explanation is difficult to be verified directly in theory, but it can be inferred that if the addition of two interfaces during phase changes does result in the discontinuities in the heating factor, the number of jump discontinuities must be consistent with the number of interfaces. That is the function Qηr′min=FHt that has only two discontinuities, regardless of whether the process limitations are met. Let’s expand the range of Ht even further, as shown in Figure 12.Figure 12
Effects of four process constraints in the range ofΔHt≤60∘C: (a) peak limiting, (b) slope limiting, (c) time limiting when T>217∘C, and (d) heating time limiting when 150∘C<T<190∘C.
(a)(b)(c)(d)Obviously, this inference is correct. It can be concluded that the jump discontinuities appear due to the increase of two interfaces between the solid and the liquid phases during phase transitions and the process limitations, but the number of discontinuity points is only related to that of interfaces between two phases.


## 4. Discussions

Results of the temperature profile and distribution of the soldering area are capable of demonstrating that our approximation and symmetry treatment are feasible. Furthermore, the dynamic Thomas algorithm used to figure out the tridiagonal system of linear equations with coefficient matrix evolving over time is accurate and efficient. As expected, there are two big differences in the temperature profile, partsPa and Pb in Figure 6(a). This phenomenon may be explained by the fact that the phase transformation process has the effect of slowing down the temperature change, so the temperature decreases faster in the case of ignoring phase transitions than the other. Moreover, the optimization results have indicated that the values of obj1. Qη and obj3. Δt leaving phase changes out are smaller than the other, though its distribution is more scattered. Again, the culprit is the phase transition that makes the slope of the temperature profile smaller when considering phase transformation in the cooling process. Therefore, the time Δt=t217,R−t217,L and area Qη=∫t217,Ltpeakft,vp−217dt become larger when taking phase changes into account. The reason why the results of obj2. Peak are similar is that the peak of the temperature profile is independent of the slope. In other words, the peaks are in accordance no matter you consider phase transitions or not.Finally, the change of the heating factor with the highest set temperatureQηr′min=FHt, in the case of considering phase changes, is not continuous but in a band-like distribution. This conclusion does not mean that the heating factor itself, Qηr′min=Gun,vp, is discontinuous. Instead, in view of its definition, Qηr′min=Gun,vp should be continuous, and you can always reasonably adjust un and vp to make Qηr′min sample all the possible values.


## 5. Conclusions

A mathematical model is developed in this study, including the heat transfer mechanism of reflow soldering and the optimization of experimental parameters and temperature profiles. The results calculated in the two cases, that is, considering and ignoring the solder phase transformations, can be summarized as follows. Firstly, temperature profiles are quite different after the start of the phase transition. Secondly, the results of multiobjective optimization are obviously dissimilar. The optimal solution with ignoring phase transition is relatively dispersed, while the results with considering phase transformation are very concentrated. Thirdly, the heating factorQηb′min is a continuous function of Ht, but Qηr′min has two jump discontinuities and shows a banded distribution. Thus, we can make the following concluding points. Phase transitions cause a tremendous impact on the cooling rates of soldering and the optimal recipe for control parameters, but it is independent of the peak temperature. This finding provides a deeper insight into the optimization and adjustment of temperature profiles. Moreover, it is the limitations related to the slope of temperature profiles and two additional interfaces that result in discontinuities of the heating factor that is a function of the maximum set temperature Qηr′min=FHt, when phase transitions are taken into account. The universality of this conclusion and the deeper reasons deserve further study. Especially, the effect of phase transitions is reflected not only in the temperature profile but also in the continuity of some control parameters. From the computing method point of view, the dynamic Thomas algorithm proposed for the first time was proved to work out the problem accurately and effectively that the coefficient matrix of the tridiagonal system of linear equations evolves with time. In numerical calculation, this algorithm is helpful to resolve more complex diagonal matrix problems, such as the five-diagonal and seven-diagonal matrices related to time. On the other hand, analyses and evaluations regarding the multiobjective optimization results of experimental parameters worth further exploration, since evaluating the convergence, uniformity, and spread has invariably been the hot topic in multiobjective optimization for result analysis.



---

*Source: 9955967-2021-06-17.xml*