#### 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.

#### 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 is approximately linear with changes in combined parameters , 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

Figure 1 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.

For the development of the mathematical model, we take a reflow oven with 11 temperature zones as an example. As shown in Figure 1, 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 Figure 2. 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.

###### 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 values . 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]:where , is the specific heat at a constant pressure of influent gas, is the density of influent gas, is the thermal conductivity of influent gas, is the temperature, and 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:where , is the coordinates of the right boundary of the -th small temperature region, and is the coordinates of the left edge of the -th region. The analytic solution of the equation is as follows:where

Obviously, when the two boundary conditions have the same temperature , the solution of the equation is , 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:where is the natural convection heat transfer coefficient of influent gas, is the location coordinate of the right end of the front area, and is the ambient temperature. As a result, the temperature distribution of the front area is as follows:where

###### 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:

The analytic solution of this equation is as follows:where

##### 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:where , is the constant pressure specific heat capacity of the circuit board, is the density of the circuit board, and is the thermal conductivity of the circuit board. The initial condition is . 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, and . Then, according to Fourier’s law of heat conduction [29], the heat flowing in from the boundaries can be expressed as and . Since the quantity of heat presented in the two ways should be equal, we have the following equation:where represents the convective heat transfer coefficient between the circuit board and the gas in the furnace, is the position coordinates of the soldering area when the circuit board in its initial position, is the thickness of the PCB, stands for the passing speed of the conveyor belt, and denotes the temperature at 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, and 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. Figure 3 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. and 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.

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.where and ; and represent the temperature distribution of the liquid and solid phases, respectively; and denote the constant pressure heat capacity of liquid and solid solder, respectively; and signify the density of liquid and solid solder, respectively; and 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:where and represent the heat flux in the positive direction of the liquid and solid phases, respectively, and is the latent heat of solder melting. According to Fourier’s law [29], Equation (14) can be interpreted as follows:

The other issue is the cooling process, that is, the solder changes from liquid to solid phase as the temperature decreases. Figure 4 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.

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:

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 distribution of each zone calculated; the heat conduction problem can be resolved to obtain the temperature variation of the soldering area. Then, the temperature variation of the center of the soldering area is , 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:where , , and . The boundary conditions (Equation (12)) were discretized using backward difference for the upper and forward difference for the lower.where . The initial condition is . 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 . Substitute the boundary conditions, and the matrix form can be expressed as follows:where . The system of linear equations can be recorded as , and the nonzero elements in the coefficient matrix 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:where , , and . is the temperature of the liquid phase, and is the temperature of the solid phase. The cooling process is opposite, and equations are similar.where . 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:where and . Similarly, the interface conditions during cooling are as follows:

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:where , , , and when melting, and , , , and when cooling. Substitute the transformation boundary conditions (Equations (22) and (23)) that and during melting and and during cooling, and the following relation can be obtained:where and during the melting process, and and 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 . , and can be written in the following form:

Although the coefficient matrix changes with time, the system still satisfies the condition of the Crout decomposition method, that is, (i) , (ii) , , , and (iii) . Obviously, is a nonsingular matrix whose determinant . Thus, its sequential principal minor of any order .

A tridiagonal square matrix satisfying the above conditions has a unique Crout decomposition matrix as follows:where , , and . and are time-varying and piecewise functions, and the recurrence formulas can be expressed as follows:

Therefore, can be decomposed into two triangular equations and . The first step is to solve the system of equations as follows:

The recurrence formulas of solutions are as follows:

The second step is to solve the system of equations as follows:

The recurrence formulas of solutions are as follows:

where 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.

#### 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 Table 1.

Temperature profiles and their relative errors are calculated, as shown in Figures 6(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.

**(a)**

**(b)**

In the region marked by 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 . 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 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 . Figure 6(b) indicates that the relative error mostly fluctuates around 1.0%. However, the error becomes larger for , even up to 10%, which is due to the fact that the actual temperature is low for . 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 .

The overall temperature distributions of the soldering area can also be calculated, for both considering and ignoring phase transitions, as shown in Figures S1 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 Table 2, the lasting time that the temperature of the center of soldering area exceeds , , should not be too long, and the peak temperature should not be too high. We now define the integration: , where is the temperature function plotted in Figure 6. Note that should reach a minimum value for an ideal temperature profile.

Multiobjective 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.

**(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 (Figure 7(c)), although its optimization results of obj1. and obj3. 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 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 Figure 8, where both the objective function and decision variable 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.

##### 3.3. Verification of the Linear Relationship between and

is the heating factor; is optimal; and 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 , , , , , and , and with considering phase transitions, we have , , , , , and . Then, we establish the functional relationships according to the reference that , , , and for ignoring phase transformations, and , , , and for considering phase transitions. When we also make small modifications to , , and keep the other parameters unchanged, the relationship between and can be obtained, as shown in Figure 9.

Obviously, is linearly related to in the two cases for the segments of and , respectively, as shown in Figure 9 when . This result is consistent with the conclusion of previous work [4]. It needs to be emphasized that the conclusion, the linear relationship between and , 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 and , it can be found that although the variation trends of them tend to be nonlinear outside the range , the difference is that changes continuously, while has an obvious discontinuity point when the temperature is about or . Therefore, it can be concluded that phase transitions break up in some places, which is supposed to change continuously, and form a banded distribution. Hence, we further expand the range of , , and take more dense points to show the banded distribution of in Figure 10, calculated by the model of involving phase transformations.

Allowed bands are the places where can be continuously sampled, and forbidden bands are the areas where is discontinuous. In general, can always be determined by changing the value of , so is continuous disregarding the phase transitions. A new phenomenon has been found that under the condition of phase changes, cannot sample all the values continuously except the points in the range of allowed bands.

Based on our analysis, the discontinuity of the heating factor should be related to the constraint conditions of the process limitations, because the value of is further expanded on the basis of the results of the constrained multiobjective optimization. The influence of each constraint on can be revealed separately by calculating the range of allowed by process limitations shown in Table 2.

From Figures 11(b) and 11(c), it can be found that the first jump discontinuity, where is about , is caused by the violation of the restriction that the time when on the temperature profile must be , and the second jump discontinuity, where is about , is resulted from going against two limitations that the rate of temperature change must be , and the heating time must be when . However, there is no connection between the other two constraints and the phenomenon that there are discontinuities in function , as shown in Figures 11(a) and 11(d).

**(a)**

**(b)**

**(c)**

**(d)**

In summary, feedback of the heating factor 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 rows, while the original has 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 and 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 that has only two discontinuities, regardless of whether the process limitations are met. Let’s expand the range of even further, as shown in Figure 12.

**(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, parts and 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. and obj3. 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 and area 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 temperature , 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, , is discontinuous. Instead, in view of its definition, should be continuous, and you can always reasonably adjust and to make 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 factor is a continuous function of , but 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 , 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.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the National Natural Science Foundation of China (Grant no. 11901329), the Natural Science Foundation of Shandong Province (Grant no. ZR2019BA022), and Teaching Reform Projects of Qingdao University of Technology F2019–044 and F2020-079.

#### Supplementary Materials

Figures S1 and S2: temperature distribution of the soldering area calculated for both considering and ignoring phase transitions, respectively. Figure S3: the results that the temperature distribution in the case of phase transitions subtracts the other. Tables S1 and S2: three objective functions calculated by numerical simulations for both ignoring and considering phase transformations, respectively, and the results of five decision variables can also be seen in Tables S3 and S4.* (Supplementary Materials)*