Multi-Physics Fields Simulations and Optimization of Solder Joints in Advanced Electronic Packaging

: The endurability of solder joints in the ball-grid array (BGA) packaging is crucial to the functioning of the microelectronic system. To improve electronic packaging reliability, this paper is dedicated to numerically optimize solder joint array conﬁguration and study the inﬂuence of multi-physical ﬁelds on solder joint reliability. The uniqueness of this study is that on the basis of temperature ﬁeld and stress ﬁeld, the electric ﬁeld is added to realize the coupling simulation of three physical ﬁelds. In addition, the “Open Angle” is mathematically deﬁned to describe the array conﬁguration, and it was used to reveal the inﬂuence factors of solder joint fatigue, including stress, temperature, and current density. In the single solder joint model, the impacts of geometric shape and working conditions on the maximum value and distribution of these evaluation factors (stress, temperature, and current density) were investigated. Overall, the numerical investigation gives the optimal conﬁguration, geometric shape, and working condition of solder joints, which beneﬁts the design of endurable and efﬁcient


Introduction
In recent decades, microelectronic chips were continuously shrinking in size [1,2] and increasing in numbers due to the urgent demand for high-speed data storage and operation. With the advancement of extreme ultraviolet lithography, the size of a single transistor reached 2 nm, which is almost the physical limit. To further improve chip performance, integrated circuit packaging became a potential way, which enables a higher density of transistors. The ball grid array (BGA) packaging is the most widely popular packaging method since the 1990s [3]. The BGA includes fine-pitch ball-grid array (FBGA), Wafer Level Chip Scale Packaging (WL-CSP) [4], and etc. Although these styles vary in material or solder joint arrangement, the basic structure is that the silicon-based chip and printed circuit board (PCB) are connected by solder joints [5].
The durability of BGA packaging heavily depends on the soldering, and the failure of a few solder joints can directly cause the malfunction of the system [6]. In order to elongate the lifespan of chips, main reasons for solder joint failure are investigated, including thermo-migration, electro-migration, and phase coarsening. Thermo-migration happens under the temperature gradient across the two joint ends, and it causes the solder alloy to become inhomogeneous [7][8][9]. Electro-migration happens under the electric current density gradient across the horizontal section, and it causes void nucleation on the cathode polar and mass accumulation on the anode polar [10][11][12]. Finally, the phase coarsening level is positively proportioned to the electric current density, causing the mechanical property degradation of the solder alloy [13,14]. The vacancies caused by the migration and the coarsening are detrimental to the microstructure of the solder joints. Intrinsically, the quality of the solder joint is directly affected by the temperature gradient and electric current gradient.
Several studies have investigated the fatigue of solder joints regarding different factors and geometries, using both theoretical analysis and simulation. Liu et al. [15] and Several studies have investigated the fatigue of solder joints regarding different factors and geometries, using both theoretical analysis and simulation. Liu et al. [15] and Peng et al. [16] derived the expression of force exerted on the solder joint, such function explains the fatigue mechanism. Wu et al. [17] applied cyclic thermal loading to a solder joint and derived its fatigue model. Furthermore, Tang et al. [18] investigated the temperature field distribution across the ball-grid array and in a single solder joint. In addition, Lau et al. [19] and Depriver et al. [20] found that simulation input parameters significantly influence the simulation results, such as convective heat transfer coefficient and solder joint material. Fan et al. [21] compared average peel stress with cumulated equivalent creep strain in the solder joint, a supplement to the von mises stress in previous studies. Furthermore, Ye et al. [10] proposed the surface marker movement technique to track electro-migration. However, few researchers took the electric current and current density distribution into consideration. Actually, the joule heat is the primary heat source, so the coupling mechanism of multi-physics fields of solder should be revealed, including the temperature field, mechanical field, and electrical field. In addition, the effect of geometric shape and working conditions on the solder joint performance was also less investigated.
This paper aims to find the most fragile joints in a solder joint array and then optimize the geometric shape and working condition of a single solder joint. The study uses software COMSOL Multiphysics as simulation platform. The performance of the solder joint is numerically evaluated by influence factors (stress, temperature, and current density). The solder failure was directly caused by stress, which was induced by temperature gradient. The simulation included coupling multi-physics modules of current density, solid mechanics, and heat transfer. Initially, the impacts of array configuration on the multiphysics fields in the solder joints were investigated. A theory of "open angle" was defined to compare solder joints within one or more configurations. Furthermore, to improve the endurability of central solder joints, single solder joint models were studied. The impacts of geometric shapes (solder ball curvature, interface area, and height) and working conditions (input current, chip temperature, and inner heat convection coefficient) on evaluation factors were studied, and the optimal working condition of the solder joints was concluded.

Physical Model
Firstly, a single solder joint model was established. The solder joint geometry had three components: a ball-shaped solder joint in the middle (the solder joint), a cuboid on the top (the chip), and a cuboid at the bottom (the PCB). The solder joint ( Figure 1a) was a ball cut horizontally on top and bottom to leave out the interfaces between the solder joint and the boards. The geometric shape was described by three controllable, independent factors: concave or convex extent, interface area, and height. Detailed definitions are discussed in Section 4. When one factor was modified, others remained unchanged.   The model of the entire ball-grid array (Figure 1b) was the expanded version of the single solder joint model. In order to describe the location of solder joint in a ball-grid array, the theory of "Open Angle" was proposed based on the array model by Edwards et al. [22]. In their simulation, a spatially periodic square array of circular cylinders was described by a set of position vectors. The open angle, in addition to square arrays, is an absolute standard applicable to the solder joint comparison between any array configuration.
As in Figure 2, a polar coordinate system is used to describe the open angle, and the solder joints are simplified as points on the grid. The coordinate system is centered in the targeted solder joint (Figure 2). In the angle aggregate θ, each element is defined as the degree of clockwise rotation of a line that connects the origin and an outermost point until the line encounters another point. The maximum angle in θ is the "Open Angle" (45 • in Figure 2). The number of elements of θ equals the number of solder joints in the outmost ring, which are defined as the solder joints that have open angles greater or equal to π. In a symmetric solder joint array, the degree of a point's open angle is equivalent to its distance to the outermost ring in the array. According to the definition of open angle, solder joints in the same ring have equal open angles.
Chips 2022, 2, FOR PEER REVIEW The model of the entire ball-grid array (Figure 1b) was the expanded version of th single solder joint model. In order to describe the location of solder joint in a ball-gri array, the theory of "Open Angle" was proposed based on the array model by Edward et al. [22]. In their simulation, a spatially periodic square array of circular cylinders wa described by a set of position vectors. The open angle, in addition to square arrays, is a absolute standard applicable to the solder joint comparison between any array configura tion.
As in Figure 2, a polar coordinate system is used to describe the open angle, and th solder joints are simplified as points on the grid. The coordinate system is centered in th targeted solder joint ( Figure 2). In the angle aggregate θ, each element is defined as th degree of clockwise rotation of a line that connects the origin and an outermost point unt the line encounters another point. The maximum angle in θ is the "Open Angle" (45° i  The solder joint is Sn62Pb36Ag2, the chip is silicon, and the PCB is FR-4. The param eters of these materials are listed in Table 1.  The solder joint is Sn62Pb36Ag2, the chip is silicon, and the PCB is FR-4. The parameters of these materials are listed in Table 1. Chips 2022, 1

194
The melting temperature of the solder material is approximately 456 K, which is higher than the usual chip work temperature of 333.15 K [19]. Therefore, the simulation does not take phase change into consideration. The residual stress is another factor that contributes to solder joint failure [23,24], usually caused by moisture cycling and coefficient of thermal expansion (CTE) mismatch on the interface. However, the residual stress is ignored in this study, because assumingly the solder joint has not gone through any thermal loading prior to the simulation.

Governing Equation
As shown in the following part, there are a series of governing equations to illustrate the multi-physics field heat transfer, electricity, and solid mechanics.
The theoretical equations are the basis of simulation mechanism. The heat conduction equation together relates the heat flow, the temperature gradient and the heat generated in the solder joint. This differential equation has the general form as the following: In this equation, λ is heat conduction coefficient, .
q v the heat generating power per unit volume, C p the specific heat capacity, ρ the density, and τ the time. The left-hand side is the net heat flow in three directions and the generated heat, which equals the change of heat with respect to time (the right-hand side). From the homogeneity of the solder joint material, we assume that C p and ρ are constants, so a simplified form of the equation will be: The heat convection equation describes the heat exchange between the solder joint and the surrounding air: In this equation, .
q is the change of heat with respect to time, h the heat convection coefficient, A the surface area, and T r the reference temperature (293.15 K in this paper). α = 1 in this paper, because the working temperature of chips is usually under 360 K, where the Newton's Law of Cooling best predicts the heat convection: For the solder-joint simulation, only the steady state needs to be studied, so the direct current assumption is made. Given the direct current, two controlling equations can be applied: the conservation of electric potential and the conservation of current. The forms of the equations applicable to solid are as the followings: In this equation, J is the electric current density, σ the electric conductivity, and ϕ the electric potential.
The expansion deformation can be simulated by tangent coefficient of thermal expansion, secant coefficient of thermal expansion, or thermal strain. The temperature-dependent thermal expansion coefficient is related to temperature by the following equation: where α m is the temperature-dependent secant thermal expansion coefficient, T m denotes the strain-free state of the material as far as the measured values of α m (T) is concerned, and α r is the redefined thermal expansion coefficient dependent on the reference temperature T ref .
Owing to the deformation of the solder joint is very small (approximately 10 −7 ), the stress and strain in the solder joint are directly proportional, described by the generalized Hooke's Law [25]. The equation can be written in the Voigt notation as the following: where σ αβ and ε αβ are the stress and strain tensors respectively for which α and β equal x, y or z; C ij is the stiffness matrix. The coefficients can be further reduced by the symmetric attributes of the plane perpendicular to the z axis and the lack of shear stress. The equation can be rewritten as: Young's modulus (E) and Poisson's ratio (v) are the input parameters in the simulation. Rewriting C ij with the two parameters, the equation then becomes:

Boundary Condition
To resemble real soldering, the two interfaces are set with fixed constraints. Because the temperature of chip will naturally increase when working, an initial temperature of 333.15 K is set on the top interface of the chip. Moreover, the heat convection coefficient of the inner surfaces (inner heat convection coefficient) ( Figure 3) follows forced convection [26], which is 10 W/(m 2 ·K) [27], while that of the other surfaces is natural convection, which is 5 W/(m 2 ·K) [28].
The interface between the solder joint and the PCB is set with a constant current, while the interface between the solder joint and the chip is set with U = U x = U y = U z = 0. The electric current density in the section of the solder joint has the initial condition of J x = J y = J z = I/(πR 2 ).

Coupling of CFD Modules and Structural Model
The commercial software for coupling COMSOL was applied for unifying and relating the physics field modules. In this approach, the heat transfer module (Equations (1)-(4)), the electric current module (Equations (5) and (6)), and the solid mechanics module (Equations (7)-(10)) run simultaneously. Coupling information is linearly transmitted along the three modules during the simulation. The joule heat generated by electric current is transferred within the solder joint and dissipated through air convection. Eventually, thermal stress and strain are induced by the thermal expansion. As temperature stabilizes, the generated and dissipated heat reaches an equilibrium, and the steady state provides all results used in this study.

Coupling of CFD Modules and Structural Model
The commercial software for coupling COMSOL was applied for unifying a ing the physics field modules. In this approach, the heat transfer module (Equati (4)), the electric current module (Equations (5)-(6)), and the solid mechanics (Equations (7)-(10)) run simultaneously. Coupling information is linearly tran along the three modules during the simulation. The joule heat generated by elec rent is transferred within the solder joint and dissipated through air convection. ally, thermal stress and strain are induced by the thermal expansion. As tempera bilizes, the generated and dissipated heat reaches an equilibrium, and the stea provides all results used in this study.

Model Validation
Solid mechanics, electric current, and heat transfer modules were used in th lation. The stress, distribution of current density, and temperature fluctuation w ulated by the modules, respectively. Concerning the efficiency and accuracy of sim two mesh sizes were applied to the model. The solder joint region used tetrahed fine mesh, while the two boards used tetrahedral finer mesh. The total mesh num over 2400+ for the solder joint region. As the mesh number increased from 500 to 4 mean temperature value remained constant, and the mean stress value started to at 2000 (Figure 4), suggesting that 2400 was an appropriate mesh number.

Model Validation
Solid mechanics, electric current, and heat transfer modules were used in the simulation. The stress, distribution of current density, and temperature fluctuation were simulated by the modules, respectively. Concerning the efficiency and accuracy of simulation, two mesh sizes were applied to the model. The solder joint region used tetrahedral extra fine mesh, while the two boards used tetrahedral finer mesh. The total mesh number was over 2400+ for the solder joint region. As the mesh number increased from 500 to 4000, the mean temperature value remained constant, and the mean stress value started to stabilize at 2000 (Figure 4), suggesting that 2400 was an appropriate mesh number.

Coupling of CFD Modules and Structural Model
The commercial software for coupling COMSOL was applied for unifying an ing the physics field modules. In this approach, the heat transfer module (Equati (4)), the electric current module (Equations (5)-(6)), and the solid mechanics (Equations (7)-(10)) run simultaneously. Coupling information is linearly tran along the three modules during the simulation. The joule heat generated by elec rent is transferred within the solder joint and dissipated through air convection. ally, thermal stress and strain are induced by the thermal expansion. As tempera bilizes, the generated and dissipated heat reaches an equilibrium, and the stea provides all results used in this study.

Model Validation
Solid mechanics, electric current, and heat transfer modules were used in th lation. The stress, distribution of current density, and temperature fluctuation w ulated by the modules, respectively. Concerning the efficiency and accuracy of sim two mesh sizes were applied to the model. The solder joint region used tetrahedr fine mesh, while the two boards used tetrahedral finer mesh. The total mesh num over 2400+ for the solder joint region. As the mesh number increased from 500 to 4 mean temperature value remained constant, and the mean stress value started to s at 2000 (Figure 4), suggesting that 2400 was an appropriate mesh number.    197 minimum temperature are at the edge. The decreasing slopes of stress and temperature increase with the solder joint deviating from the center. This same tendency of stress and temperature verified that the stress is in fact induced by temperature variation.

Open Angle Theory Validation
For ordered point symmetric configurations, the simulation proved a prediction made by the open angle theory. The isolines of stress and temperature are in a ring structure. The ring structure suggests that the variance of evaluation factors is little within one ring but huge between rings. To test the hypotheses, 3 × 3, 5 × 5 and 6 × 6 ball-grid arrays were simulated. The mean stress and temperature distributions in the simulations ( Figure 6) proved the hypotheses correct. The theory also predicts that the values evaluation factors should be equal within one ring. However, the numerical results showed that the solder joints at the corners have lower stress and temperature than those on the sides. To explain the difference, a Second (largest) Open Angle additional to the Open Angle can be involved to correct the theory.

General Distribution Trend
The 16 × 16 ball-grid array model was adopted. It can be seen from Figure 5 th maximum stress and maximum temperature are at the center and minimum stress a minimum temperature are at the edge. The decreasing slopes of stress and temperatu increase with the solder joint deviating from the center. This same tendency of stress a temperature verified that the stress is in fact induced by temperature variation.

Open Angle Theory Validation
For ordered point symmetric configurations, the simulation proved a predicti made by the open angle theory. The isolines of stress and temperature are in a ring stru ture. The ring structure suggests that the variance of evaluation factors is little within o ring but huge between rings. To test the hypotheses, 3 × 3, 5 × 5 and 6 × 6 ball-grid arra were simulated. The mean stress and temperature distributions in the simulations (Figu 6) proved the hypotheses correct. The theory also predicts that the values evaluation f tors should be equal within one ring. However, the numerical results showed that t  The open angle theory also reveals comparison between configurations, beyon single solder joint. A weighted mean of the open angle of all solder joints was use cause of the nonlinear increase of stress and temperature from outer rings to inner r Consequently, the small-open-angle solder joints in a configuration greatly diminis scoring coefficients in the weighted averaging process. Eventually, the configuration the greatest mean open angle is the most endurable.
Overall, when designing solder joint configurations, restrictions, such as num shape, dimensions and joint distance, are differently set in different cases, which m the arbitrary best configuration does not exist. However, the guiding rules that can to the best design in every occasion are concluded. First, in an ordered configur where rings are clear, the number of the ring should be minimized. Second, increasin edges of the polygonal configuration outperforms adding external solder joints. T each solder joint's contribution to the decrease of others' open angles should be m mized. In this way, the optimal configuration problem is changed into a 2D dynamic gramming problem. Overall, when designing solder joint configurations, restrictions, such as number, shape, dimensions and joint distance, are differently set in different cases, which means the arbitrary best configuration does not exist. However, the guiding rules that can lead to the best design in every occasion are concluded. First, in an ordered configuration where rings are clear, the number of the ring should be minimized. Second, increasing the edges of the polygonal configuration outperforms adding external solder joints. Third, each solder joint's contribution to the decrease of others' open angles should be minimized. In this way, the optimal configuration problem is changed into a 2D dynamic programming problem.

Geometric Shape-Curvature
The side of the solder ball can be either concave (inward contraction) or convex (outward extension), corresponding to negative or positive curvature. As the curvature of the side increases, the solder ball deviates the cylindrical shape. The evaluation factors show two distinguished clusters at positive curvatures and negative curvatures, so it is clearer to separately view the curves of convex and concave shapes in two graphs.
Overall, the value range of stress is evidently lower for concave curvatures, and the value ranges of temperature and current density are evidently lower for convex curvatures ( Figure 8).
Consider the convex curvatures. In the range of [7067.1 m −1 , 10,929.0 m −1 ], as curvature increases, maximum stress generally increases, maximum temperature decreases almost linearly, and maximum current density fluctuates with a maximum at approximately 9500 (m −1 ). Consider the concave curvatures. In the negative range, as the absolute value of the curvature decreases, maximum stress approximately linearly increases, and maximum temperature and current density decreases with a decreasing slope.

Geometric Shape-Curvature
The side of the solder ball can be either concave (inward contraction) or convex (outward extension), corresponding to negative or positive curvature. As the curvature of the side increases, the solder ball deviates the cylindrical shape. The evaluation factors show two distinguished clusters at positive curvatures and negative curvatures, so it is clearer to separately view the curves of convex and concave shapes in two graphs.
Overall, the value range of stress is evidently lower for concave curvatures, and the value ranges of temperature and current density are evidently lower for convex curvatures (Figure 8).
Consider the convex curvatures. In the range of [7067.1 m −1 , 10,929.0 m −1 ], as curvature increases, maximum stress generally increases, maximum temperature decreases almost linearly, and maximum current density fluctuates with a maximum at approximately 9500 (m −1 ). Consider the concave curvatures. In the negative range, as the absolute value of the curvature decreases, maximum stress approximately linearly increases, and maximum temperature and current density decreases with a decreasing slope. In each of the concave stress section figures (Figures 9-11), the absolute value of the curvature increases from the left to the right. In all figures, high-stress and high-currentdensity areas are at the corners (except for the current density figures of concave curvatures), and the temperature distribution stratifies. For the convex curvatures, the stress in the central area decreases, while the area of high-stress sections at the corners (or the interfacial rims in 3D perspective) increases as the absolute value of curvature increases. For the concave curvatures, the stress in the central area increases, while that at the corners decrease as the absolute value of curvature increases (Figure 9). For both curvature types, In each of the concave stress section figures (Figures 9-11), the absolute value of the curvature increases from the left to the right. In all figures, high-stress and highcurrent-density areas are at the corners (except for the current density figures of concave curvatures), and the temperature distribution stratifies. For the convex curvatures, the Chips 2022, 1 201 stress in the central area decreases, while the area of high-stress sections at the corners (or the interfacial rims in 3D perspective) increases as the absolute value of curvature increases. For the concave curvatures, the stress in the central area increases, while that at the corners decrease as the absolute value of curvature increases (Figure 9). For both curvature types, the temperature is evenly stratified, and the hotter layer extrudes the cooler layers as the absolute value of curvature increases ( Figure 10). The high-current-density areas are at the corners for convex curvatures and in the center for concave curvatures, which increase as the absolute value of curvature increases. The current density in the rest of the area, however, decreases as the absolute value of curvature increases (Figure 11).
Chips 2022, 2, FOR PEER REVIEW 11 absolute value of curvature increases ( Figure 10). The high-current-density areas are at the corners for convex curvatures and in the center for concave curvatures, which increase as the absolute value of curvature increases. The current density in the rest of the area, however, decreases as the absolute value of curvature increases ( Figure 11).   Chips 2022, 2, FOR PEER REVIEW 11 absolute value of curvature increases ( Figure 10). The high-current-density areas are at the corners for convex curvatures and in the center for concave curvatures, which increase as the absolute value of curvature increases. The current density in the rest of the area, however, decreases as the absolute value of curvature increases ( Figure 11).   Chips 2022, 2, FOR PEER REVIEW 12 Figure 11. Current density distribution with respect to concave and convex curvatures.

Geometric Shape-Interface Area
The interface includes the upper and lower surfaces of the solder joints. As shown in Figure 12, in the range of [0.04011 mm 2 , 0.111702 mm 2 ], maximum stress fluctuates and its increases with an increasing slope. Maximum temperature and maximum current density show a smooth, decelerated descent.

Geometric Shape-Interface Area
The interface includes the upper and lower surfaces of the solder joints. As shown in Figure 12, in the range of [0.04011 mm 2 , 0.111702 mm 2 ], maximum stress fluctuates and its increases with an increasing slope. Maximum temperature and maximum current density show a smooth, decelerated descent.
Chips 2022, 2, FOR PEER REVIEW 12 Figure 11. Current density distribution with respect to concave and convex curvatures.

Geometric Shape-Interface Area
The interface includes the upper and lower surfaces of the solder joints. As shown in Figure 12, in the range of [0.04011 mm 2 , 0.111702 mm 2 ], maximum stress fluctuates and its increases with an increasing slope. Maximum temperature and maximum current density show a smooth, decelerated descent. In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the interfacial area increases, section area of the solder joint increases, which causes the stress ( Figure 13) to decrease, the temperature to increase (Figure 14), and the current density to decrease ( Figure 15) in the central area. The area of the high-stress sections at the corners increases.    In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the interfacial area increases, section area of the solder joint increases, which causes the stress ( Figure 13) to decrease, the temperature to increase (Figure 14), and the current density to decrease ( Figure 15) in the central area. The area of the high-stress sections at the corners increases. In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the interfacial area increases, section area of the solder joint increases, which causes the stress ( Figure 13) to decrease, the temperature to increase (Figure 14), and the current density to decrease ( Figure 15) in the central area. The area of the high-stress sections at the corners increases.    In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the interfacial area increases, section area of the solder joint increases, which causes the stress ( Figure 13) to decrease, the temperature to increase (Figure 14), and the current density to decrease ( Figure 15) in the central area. The area of the high-stress sections at the corners increases.

Geometric Shape-Height
The height is defined as the distance between the two interfaces. As shown in Figure 16, in the range of [0.172 mm, 0.217 mm], maximum stress fluctuates and its increases with an increasing slope, maximum temperature slowly increases, and current density linearly increases.

Geometric Shape-Height
The height is defined as the distance between the two interfaces. As shown in Figure  16, in the range of [0.172 mm, 0.217 mm], maximum stress fluctuates and its increases with an increasing slope, maximum temperature slowly increases, and current density linearly increases.   In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the height increases, the stress in the central area increases, while that at the corners remains unchanged, and the rest of the section is filled with low-stress sections ( Figure 17). The temperature is barely affected by the change of the height (Figure 18). The current density across the section area uniforms as the height increases ( Figure 19).
Chips 2022, 2, FOR PEER REVIEW 15 In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the height increases, the stress in the central area increases, while that at the corners remains unchanged, and the rest of the section is filled with low-stress sections ( Figure 17). The temperature is barely affected by the change of the height (Figure 18). The current density across the section area uniforms as the height increases ( Figure 19).

Comparison of Impact
The impacts of the three geometric factors on the variation of the evaluation factors were compared. The impact on maximum stress decreases in the sequence of the interface area, the height, and finally within the convex or concave curvature range. The impact on maximum temperature is extremely high for concave curvature and similar for all the other geometric factors. The impact on maximum current density decreases in the sequence of the concave curvature, interface area, and the convex curvature.
Therefore, the optimal geometric shape can be concluded. In a given dimension, the optimal shape a solder joint would consist of a small and convex curvature, a medium or large interface area, and a low height. These results are intuitive to understand. As the curvature increases and the solder joint deviates from cylindrical shape, stress and current density concentrates at the interfacial rim and decreases in the center. As the interface area increases, the same amount of current and heat flow into a larger tunnel, and the average value naturally decreases. As the height increases and the distance between the chip and Chips 2022, 2, FOR PEER REVIEW 15 In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the height increases, the stress in the central area increases, while that at the corners remains unchanged, and the rest of the section is filled with low-stress sections ( Figure 17). The temperature is barely affected by the change of the height (Figure 18). The current density across the section area uniforms as the height increases ( Figure 19).

Comparison of Impact
The impacts of the three geometric factors on the variation of the evaluation factors were compared. The impact on maximum stress decreases in the sequence of the interface area, the height, and finally within the convex or concave curvature range. The impact on maximum temperature is extremely high for concave curvature and similar for all the other geometric factors. The impact on maximum current density decreases in the sequence of the concave curvature, interface area, and the convex curvature.
Therefore, the optimal geometric shape can be concluded. In a given dimension, the optimal shape a solder joint would consist of a small and convex curvature, a medium or large interface area, and a low height. These results are intuitive to understand. As the curvature increases and the solder joint deviates from cylindrical shape, stress and current density concentrates at the interfacial rim and decreases in the center. As the interface area increases, the same amount of current and heat flow into a larger tunnel, and the average value naturally decreases. As the height increases and the distance between the chip and Chips 2022, 2, FOR PEER REVIEW 15 In all figures, high-stress and high-current-density areas are at the corners, and the temperature distribution stratifies. As the height increases, the stress in the central area increases, while that at the corners remains unchanged, and the rest of the section is filled with low-stress sections ( Figure 17). The temperature is barely affected by the change of the height (Figure 18). The current density across the section area uniforms as the height increases ( Figure 19).

Comparison of Impact
The impacts of the three geometric factors on the variation of the evaluation factors were compared. The impact on maximum stress decreases in the sequence of the interface area, the height, and finally within the convex or concave curvature range. The impact on maximum temperature is extremely high for concave curvature and similar for all the other geometric factors. The impact on maximum current density decreases in the sequence of the concave curvature, interface area, and the convex curvature.
Therefore, the optimal geometric shape can be concluded. In a given dimension, the optimal shape a solder joint would consist of a small and convex curvature, a medium or large interface area, and a low height. These results are intuitive to understand. As the curvature increases and the solder joint deviates from cylindrical shape, stress and current density concentrates at the interfacial rim and decreases in the center. As the interface area increases, the same amount of current and heat flow into a larger tunnel, and the average value naturally decreases. As the height increases and the distance between the chip and

Comparison of Impact
The impacts of the three geometric factors on the variation of the evaluation factors were compared. The impact on maximum stress decreases in the sequence of the interface area, the height, and finally within the convex or concave curvature range. The impact on maximum temperature is extremely high for concave curvature and similar for all the other geometric factors. The impact on maximum current density decreases in the sequence of the concave curvature, interface area, and the convex curvature.
Therefore, the optimal geometric shape can be concluded. In a given dimension, the optimal shape a solder joint would consist of a small and convex curvature, a medium or large interface area, and a low height. These results are intuitive to understand. As the curvature increases and the solder joint deviates from cylindrical shape, stress and current density concentrates at the interfacial rim and decreases in the center. As the interface area increases, the same amount of current and heat flow into a larger tunnel, and the average value naturally decreases. As the height increases and the distance between the chip and the PCB increases, the medium-stress area in the center of the solder joint increases, causing greater reaction at the interfacial rim to counteract the increasing stress.

Working Condition
The impact of working conditions on evaluation factors (stress, temperature, and current density) were investigated in this part. In each graph, different curves represented the same geometric shape under different working conditions.
In the graphs of curvature, interface area, and height, as input current increases from 1 A to 5.5 A, the maximum stress, temperature, and current density curves all translate upward, confirming positive relations ( Figure 20); as the chip temperature increases from 328.15 K to 368.15 K, the maximum stress curves translate upward, confirming again a positive relation ( Figure 21); and as the inner heat convection coefficient increases from 10 to 60, the maximum stress curves translate downward, confirming a negative relation ( Figure 22). As curvature, interface area, or height increases, the impact of input current on stress decreases (Figure 20).
The optimal working condition for solder joints can be concluded. Within the restriction ranges, the boundary condition should have a low input current, a low chip temperature, and a high inner heat convection coefficient. All of the three conditions refer to less generated and more dissipated heat, which benefit the endurability of the solder material. Among the three boundary conditions, the chip temperature's impact on curve translation is the greatest, while the inner heat convection coefficient's impact is the least. This result suggests that the chip power and the chip heat-transfer ability deeply affect the durability of the solder joints beneath; an effective chip cooling system is the most effective improvement to the solder joints' working condition.
Chips 2022, 2, FOR PEER REVIEW 16 the PCB increases, the medium-stress area in the center of the solder joint increases, causing greater reaction at the interfacial rim to counteract the increasing stress.

Working Condition
The impact of working conditions on evaluation factors (stress, temperature, and current density) were investigated in this part. In each graph, different curves represented the same geometric shape under different working conditions.
In the graphs of curvature, interface area, and height, as input current increases from 1A to 5.5A, the maximum stress, temperature, and current density curves all translate upward, confirming positive relations ( Figure 20); as the chip temperature increases from 328.15 K to 368.15 K, the maximum stress curves translate upward, confirming again a positive relation ( Figure 21); and as the inner heat convection coefficient increases from 10 to 60, the maximum stress curves translate downward, confirming a negative relation ( Figure 22). As curvature, interface area, or height increases, the impact of input current on stress decreases ( Figure 20).
The optimal working condition for solder joints can be concluded. Within the restriction ranges, the boundary condition should have a low input current, a low chip temperature, and a high inner heat convection coefficient. All of the three conditions refer to less generated and more dissipated heat, which benefit the endurability of the solder material. Among the three boundary conditions, the chip temperature's impact on curve translation is the greatest, while the inner heat convection coefficient's impact is the least. This result suggests that the chip power and the chip heat-transfer ability deeply affect the durability of the solder joints beneath; an effective chip cooling system is the most effective improvement to the solder joints' working condition.

Conclusions
Based on previous studies, this paper improved the accuracy of solder joint simulation model by concerning the electric field. This paper investigated solder joint array configurations by proposing the Open Angle theory and quantified the influential factors on the solder joint performance by curves and sectional color-gradient figures. Single solder joints in different configurations and configurations as wholes were compared by open angles. Since the optimal configuration varies as condition changes, the guiding rules of finding it are reducing ring number, decreasing central solder joint number, and minimizing the weighed mean open angle of all solder joints. However, only square configurations of solder joint array are studied in this paper, and more configuration types are needed to prove the Open Angle theory more completely. The fluctuation of evaluation factors (stress, temperature, and current density) was investigated as geometric factors (curvature, interface area, and height) and boundary conditions (input current, chip temperature, inner heat convection coefficient) changed. Furthermore, the impact from geometric factors on evaluation factors' maximum value and high-value area were respectively negative, positive, and negative; the same impact from boundary conditions were negative, negative, and positive. Noticeably, only convex curvature was adopted because it has lower evaluation factor values than its concave counterpart, and the chip temperature is the most significant influencer among the boundary conditions. However, the study assumes the constant current, while cases of fluctuating current are yet to be discussed. Therefore, the guiding rules of the optimal configuration, optimal geometric shape, and optimal working condition were concluded, proposing a strategy that designs endurable and efficient ball-grid arrays.