Decreasing Shear Stresses of the Solder Joints for Mechanical and Thermal Loads by Topological Optimization.

A methodology for obtaining the optimal structure and distribution for the gradient properties of a material in order to reduce the stress level in a soldered joint was constructed. The developed methodology was based on a combination of topological optimization methods (the moving asymptotes method) and the finite elements method; it was first implemented to solve problems of optimizing soldered joints. Using the proposed methodology, a number of problems were solved, allowing one to obtain optimal structural characteristics, in which a decrease in stress is revealed. Designing compounds using this technique will provide more robust designs. The proposed technique can be applied to a wide class of practical problems.


Introduction
Most structures in the aviation, shipbuilding and rocket industries are aimed at achieving a high resistance against various loads while being simultaneously lightweight. When designing and manufacturing mechanical (or electromechanical) structures, the choice of bonding components plays a crucial role in achieving the industry-required characteristics. Mechanical bonding, including beating and twisting, is widely used in design, although a technology of adhesion joints either treated separately or matched with mechanical fastening can essentially improve mechanical characteristics from the point of view of stiffness, strength and durability [1][2][3]. Solder and adhesive joints stand are frequently used advantageous ways of connecting different materials.
Solder and adhesive joints have advantages over other alternative ways of fastening due to stress distribution over a wider area than that based on adhesive joints or beating fastening. In contrast to the non-uniform load distribution under fastening by bonding elements, the load transmission through solder or adhesive components is continuous along the entire connecting layer. It enables the production of simple and lightweight joints. In other words, such a joint makes it possible to decrease the weight of the structure while maintaining the required mechanical loading capacity. That is why solder and adhesive joints are so often used in the design of mechanical systems [4][5][6][7][8].
Although it is well known that adhesive is one of the best methods for fastening composites with metals, the non-uniform distributions of stress/strain generated in bond-lines under shear loads often cause damage to bonded joints. However, the solder joints, suffer from temperature changes in the exploitation process due to the different linear temperature expansion coefficients of the bonded materials.
Most failures of solder joints that occur in electronic systems are problems caused by thermal mismatch between various materials, including those involved in the creation of electronic systems. During manufacture and operation, the structure goes through various temperature cycles that cause thermal expansion. Materials cannot expand freely because they are limited by packaged assembly. Therefore, significant stresses arise in soldered joints. In [9], these stresses causing the destruction of soldered joints were mathematically modeled using software and a comparison was made with experimental results. The thermomechanical reliability of the solder connection has been extensively investigated using finite element method [10,11].
Thermal stresses generated by solder are changed in wide intervals depending on the nature of the temperature fields, configurations and material properties. High stresses can cause cracks, thermoplastic deformations and other harmful effects which reduce the load carrying ability of the skip solders. It is expected that the optimization of the solder layer topology can eliminate problems associated with the occurrence of peak stresses in the intervals of the imposed structure and technological constraints. In the case of adhesive joints composed of functionally graded material/composite, stress fields are also often affected by the lack of symmetry. The maximum stress-strain characteristics are near the ends of the joints and suggest their rupture [12].
Groth and Nordlund [13] proposed the optimization of the shape of bonded joints to create strong and light joints under static loading by introducing optimal profiling of adherents. A significant decrease in the level of stress of the adhesive layer was obtained considering single/double laps and a double strap.
Hildebrand [14] analyzed single lap joints between fiber-reinforced plastics and metals using the non-linear finite element methods. Different shapes of adhesive fillet, rounded edges, reverse tapering and denting were used to increase the joint strength.
Rispeler et al. [15] employed the evolutionary structural optimization approach in order to optimize the shape of adhesive fillets in tabs of tensile test specimens. They showed a reduction in the maximum principal stresses in the adhesive for all considered cases.
Taib et al. [16] described the effects of joint configuration, defects, humidity, adhesive stiffness of glass-fiber-reinforced vinyl-ester composite laminates. In particular, they demonstrated a decrease of failure load and displacement associated with the increase of the adhesive layer thickness and aging of the joint in a hot and humid environment.
Hints and design methods to increase the load capacity of various joints and the state-of-the-art approaches aiming at the modification of the geometry of adhesive joints as well as their influence on the magnitude of stresses and effective bonding forces were considered in [3]. The properties of adhesive material, the thickness of adhesive layer, length and width of the adhesive contact layer as well as residual stresses were taken into account.
Silva and Adams [17] combined two adhesives (one for strength at high temperatures and one for strength at low temperatures) to design joints required for the fuselage of a supersonic aircraft. A numerical analysis based on finite element models made it possible to reduce the stress distribution and design the best titanium/titanium and titanium/composite double lap joints.
Haghani et al. [18] developed a stress reduction method to modify the geometry of the adhesive joint by either tapering the laminate end or by adding an adhesive fillet. An optical measurement system was used to monitor the strain field. It was shown that normal tapering of the laminate did not affect the shear and principal stress components.
Lang and Malick [19] investigated how the geometry of adhesive spew affects peak stresses and stress distribution on adhesively bonded single lap joints. A linear 2D strain analysis was carried out using the ANSYS finite element software. It was illustrated how by shaping the spew, a smoother transition in joint geometry reduced the stress concentration at the substrate-adhesive interface.
There is a long list of investigations that focused on the parametric optimization of adhesive lap joints by using appropriate change of spew and chamfer angles [16,[20][21][22][23][24][25]. Although the authors claimed that the combination of the latter two characteristics allowed them to achieve bonding in adhesive lap joints with the highest strength, the given values are not clearly defined.
Sancaktar and Simmons [26] carried out a parametric study using finite element analysis for three different model adhesives: rubber toughened film epoxy, a styrene-butadiene-styrene block copolymer, and a metal filled brittle epoxy adhesive. The experiments show that the increase in the joint strength with the introduction of notches which was in agreement with the finite element analysis.
Despite the parametric optimization carried out, non-parametric optimization, including topological optimization and optimization of the form, can be used to determine the strength of the adhesive joint.
Algorithms based on non-parametric optimization are aimed at structural optimization without taking into account the a priori chosen variables. Therefore, the basic advantage of the non-parametric optimization lies in defining the best form/topology without any a priori knowledge of the final construction design. However, the use of non-parametric approaches to the geometry of the adhesive or solder joint is not widely known and is treated rather marginally in the available literature.
The first attempt to optimize adhesive joints using the adhesive variable was described in [13]. The work was focused on achieving joints with maximum light in the static conditions by changing the profiles of adhesives. It was shown how the optimization of the form causes a significant reduction of stresses in the adhesive layer. The evolutionary method of structural optimization was applied to optimize the form of the adhesive layer in order to achieve minimum stresses [15].
Kaye and Heller [27] developed an automated sensitivity-based optimization procedure for the optimal design of free-form bonded repairs and lap joints to reduce adhesive stresses. Significant improvements were presented in relation to conventional designs.
Ejaz et al. [28] studied the applicability of non-parametric structural optimization algorithms for the optimization of adhesive joints of the types: single lap, double lap and double lap strap. Stress reduction in the adhesive layers was achieved.
In [29], ANSYS software was used to analyze the structure effect and shape of a soldered joint on fatigue life due to elastoplastic deformation of the electronic package; these factors influenced strongly the structure and shape of the solder.
Tian and Wang [30] proposed schemes for analyzing the shape and reliability of a soldered joint. By changing the pad's size, the shapes of soldered joints with different volumes of solder were predicted. Using the finite element method, the characteristics of the stresses and strains distribution in soldered joints under thermal load were analyzed.
In article [31], the shape's and height's influence of the soldered joint on the thermal load time was studied. Experimental data indicate that the solder shape stands for the dominant factor affecting the crack initiation time.
In [32], the shape and location optimization of soldered joints under vibration and shock conditions was studied in order to increase their reliability and durability. Based on the experiment, a finite element model was constructed that is close to the assembly model. Then, the modified genetic algorithm (global optimization) was used to determine the shape and location of soldered joints under shock loads. The optimum results showed that a solder joint with an optimal arrangement and shape had a lower maximum deformation, which increased the reliability of the solder joint.
In general, the optimization was focused on either the form of neighborhood elements or the form and position of the adhesive layers. The alternative strategy includes introduction of the variation/graduation of the properties along the thickness or length of the adhesive layers. It includes the modification of material properties or the geometry of adhesives which are called the functionally graded adhesive joints. In order to achieve the required properties of the adhesive, various approaches were introduced: (i) softening the brittle adhesive with rubber; (ii) supplementing additional stiffness of the flexible adhesive with glass microspheres or silica nanoparticles [24,25]; (iii) mixing of various bonded adhesives and control of the polarization level by induction heating [26][27][28]. Other variants of possible joint modifications were described in [30].
In general, the optimization was focused on either the form of neighborhood elements or the form and position of the adhesive or solder layers. One of the approaches indicates the use of a softer adhesive at the boundary ends and a more rigid one in the joint center. The second concept presents the contraction of the adhesive at the joint ends which allows for preserving neutral stresses and decreasing peak stresses in the adhesive [32,33].
In recent years, a series of works devoted to the study of adhesive properties and occurring stresses have been published [34][35][36]. The studies were carried out on the basis of analytical approaches [37,38], numerical modeling [39,40] as well as a combination of both techniques [41][42][43].
Studies based on continuous change of properties use a linearly elastic material model [35,36]. Independently of the approach for getting gradient properties of the adhesive layer, a control of the production process was employed to obtain the required optimized distribution/change of material properties. Owing to the complexity of the problem, many researchers began to reduce stress/deformation by using bi-adhesive with different stiffnesses of layers, allowing for easier control [36,41].
Brebia and Kassab [44] combined many papers with recent information regarding computer-aided structural design optimization software and integrated packages for structural optimization. The presented subjects included optimal control, shape optimization, linear and non-linear structural optimization, component reliability and optimal design of reinforced members.
París et al. [45] presented a strategy to deal with topology optimization based on minimum weight with stress constraints for topology optimization of continuum structures. The stiffness of the resulting structure was maximized for a given load case.
Qiu and Li [46] derived the KS (Keisselmeier-Steinhauser) global constraint function to reduce the members of stress constraints.
Stolpe and Svanberg [47] considered the discretized zero-one continuum topology optimization to find the optimal distribution of two linearly elastic materials in order to achieve the compliance minimum. They used a material interpolation model based on a certain rational function parameterized by a positive scalar such that the compliance was a convex function.
Le et al. [48] implemented an effective algorithm to resolve the stress-constrained topology optimization task. The applied technique combined the density filter for length scale control, solid isotropic material with penalization, and a new definition of stress to remove the singularity of stress and control local stress levels.
On the basis of the state of the art, it can be stated that even if non-parametric optimization was used, the existing optimization in the mentioned works and outside them concerned either the form of the bonded elements or the form and location of the adhesive layer.
The overall solder joint reliability was estimated by the combination of operating conditions and system design. The working environment determined the temperature extremes that the structure must withstand as well as the possibility of mechanical stress. The above studies show that under the influence of temperature and mechanical loads, structural failure occurs due to peak tangential stresses near the edges of the joint. Therefore, the present work focuses specifically on reducing shear stresses.
To the best of our knowledge, methods of topological optimization to define the optimal structure of the solder or adhesive layers in the presence of thermoelasticity were not used in the existing literature.
Despite the complexity of manufacturing topologically optimal connection layers, the structure type optimal for obtaining the maximum bond strength for a specific design is of great interest, since it determines the target solution that must be achieved.
In this paper, a mathematical model is proposed and a methodology was developed for solving a solder class of thermoelastic problems of topological optimization with the aim to obtain the optimal structure and required gradient properties of the solder layer to decrease the level of shear stresses generated in it.
Solder joint failures are a common failure mode observed in electronic packages [49]. A number of problems were considered to reduce the shear stresses of joints with silver solder in thermoelastic three-layer packets subjected to mechanical and thermal loads by topological optimization.

Problem Statement
Thermal stresses are one of the most important factors in the production of Micro Electro Mechanical Systems (MEMS), and have negative consequences due to the occurrence of harmful cracks. They also reduce the length of their exploitation with the required functional properties. The difference in the linear heat transfer coefficients associated with different materials of the package implies the occurrence of both thermal stresses and deformations. Particular attention should be paid to solder between package layers that guarantee electrical contact between the MEMS electronic components. The change in temperature within maximum and minimum values has many negative effects, including the occurrence of stoppages in MEMS.
In the following section, we present, due to its simplicity, a plane problem of thermoelasticity. We consider a package composed of an elastic non-homogenous body Ω exhibiting phase stress/strain state with regard to the x axis and consisting of three parts, as shown schematically in Figure 1 (parts Ω i are coupled and index i denotes the part number (i = 1,2,3)). The upper layer 1 stands for the foundation, the central layer 3 stands for the solder layer, while the bottom layer 2 plays the role of the base. The whole package is embedded into the temperature field T(x), x = {x 1 , x 2 }. We denote the temperature change θ(x) = T(x) -T 0 with respect to its initial value T 0 . The Young modulus and the temperature coefficient of the linear extension in Ω i are denoted by E i (x) and α i (x), respectively. The 2D elastic part of Ω is bounded by the closed surface Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 . It is assumed that we consider linear, elastic and isotropic material.

Problem Statement
Thermal stresses are one of the most important factors in the production of Micro Electro Mechanical Systems (MEMS), and have negative consequences due to the occurrence of harmful cracks. They also reduce the length of their exploitation with the required functional properties. The difference in the linear heat transfer coefficients associated with different materials of the package implies the occurrence of both thermal stresses and deformations. Particular attention should be paid to solder between package layers that guarantee electrical contact between the MEMS electronic components. The change in temperature within maximum and minimum values has many negative effects, including the occurrence of stoppages in MEMS.
In the following section, we present, due to its simplicity, a plane problem of thermoelasticity. We consider a package composed of an elastic non-homogenous body Ω exhibiting phase stress/strain state with regard to the x axis and consisting of three parts, as shown schematically in Figure 1 (parts Ωi are coupled and index i denotes the part number (i = 1,2,3)). The upper layer 1 stands for the foundation, the central layer 3 stands for the solder layer, while the bottom layer 2 plays the role of the base. The whole package is embedded into the temperature field T(x),x = {x1, x2}. We denote the temperature change θ(x) = T(x) -T0 with respect to its initial value T0. The Young modulus and the temperature coefficient of the linear extension in Ωi are denoted by Ei(x) and αi(x), respectively. The 2D elastic part of Ω is bounded by the closed surface 1 2 3 .      It is assumed that we consider linear, elastic and isotropic material. The boundary conditions are shown in Figure 1. Boundary Г1 refers to the rigid clamping, boundary Г2 is under the load of intensity t, whereas boundary Г3 is free. The field of displacement (u1, u2) implies the following equilibrium equation: where ij is the stress tensor. Linear deformations and displacements are coupled through the following relations The stress-strain relations, in the case of plane stress state, are estimated via the Duhamel-Neumann principle The boundary conditions are shown in Figure 1. Boundary Γ 1 refers to the rigid clamping, boundary Γ 2 is under the load of intensity t, whereas boundary Γ 3 is free. The field of displacement (u 1 , u 2 ) implies the following equilibrium equation: where σ ij is the stress tensor. Linear deformations and displacements are coupled through the following relations The stress-strain relations, in the case of plane stress state, are estimated via the Duhamel-Neumann principle where δ ij is the Kronecker's delta, E(x) is the Young modulus and α(x) is the linear coefficient of temperature expansion. In the case of a plane deformable state, the following substitution should be made: E(x) on the E(x) by relation (4) and α(x) on the α(x) by relation (5), where ν is the Poisson's ratio. Observe that the fields of displacement and temperature are coupled through Equation (3).

Formulation of Topological Optimization Problem
Analysis of vibrations in the solder joints shows that small solder thickness means the occurrence of shear stresses. They are concentrated at the ends of the solder and they achieve minimum at the solder center. This is why the goal of optimization is to decrease the peak values of stresses σ 12 in the solder layer. In fact, the optimization of the solder structure is focused on searching for the best solder distribution along Ω 3 in order to obtain the minimum peak values of stresses σ 12 .
Since a small number of constraints in the optimization problem plays a crucial role in computational processes, we start with only one constraint introduced on maximum stresses where r(x) is the construction variable and σ is the given value of the maximum shear stress in the solder. However, the maximum operators are not differentiable, which does not allow us to derive analytical sensitivity formulas under the occurrence of Equation (6). In order to smooth the mentioned operators, continuous functions matching local stresses to one global constraint were introduced [45,46,50].
In particular, París et al. [45] proposed the following global aggregation functions: where µ is the penalty parameter having usually large values. Then, inequality (6) can be recast to the following form: where g 1max = 1/µ ln(mes(Ω 3 )) which is obtained for σ 12 = σ in all points of space Ω 3 . The second most commonly used aggregation function for local constraints is the p norm approximation governed by the following formula: which tends to the norm (6) for p → ∞ . The numerical test using aggregation functions (7) and (9) showed that the application of function (9) ensures faster smoothing of shear stresses. That is why in our further study, we use function (9). In order to optimize the solder layer topology, we divide space Ω into finite elements. Projection variable r(x) is defined only in space Ω 3 and is coupled with Young modulus E(x) and with β(x) = E(x)α(x) as well as with the volume material density ρ(x) of each element from Ω 3 under the following relations (RAMP scheme [47]): where a, b are the penalty parameters used for the security of compact material distribution, which are selected as a result of a computational experiment, r(x) is the field of design variables 0 < r 0 ≤ r(x) ≤ 1, and r 0 is a small number which guarantees non-zero stiffness of the finite elements. For r(x) = 1 the whole volume is fully filled by the solder material.
The problem of topological optimization for the task of decreasing the peak values of stresses σ 12 in the solder layer can be formulated in the following way: under the constraint and subjected to two isoparametric constraints put on the physical density ρ(x): where γ denotes the proportion of the solder material used to build the structure and mes(Ω 3 ) is the area of the solder space.
In order to project the field of design variables into space 0/1, when the boundaries of optimal structure should be precisely defined, a Heaviside filter (HF) is used [48]. HF is a Heaviside step function which projects the field of design variables (now called the transitional field of design variables and denoted by ξ(x)) into the real field of design variables r(x). To simplify the computation of two gradients when conducting a sensitivity analysis to find a solution to the optimization problem, we use the smoothened Heaviside function of the following form: which is continuous. Parameter k ≥ 0 stands for the curvature of projection which is linear for k = 0, and tends to the unit Heaviside step function for k → ∞. Further, in order to remove the effect of "chess plane", we define a target function in the approximation process in the form of a linear combination of function (11) and additionally introduced penalty function of the following form: The second term stands for the penalty function, h 0 is the initial size of the mesh of finite elements and h max stands for the running mesh size. The quantity 0 ≤ q ≤ 1 (given coefficient) makes it possible to balance the target function and the penalty function. A solution to the so far defined optimization problem has been presented using the method of moving asymptotes [51]. We have used FEM with linear triangle FE and non-uniform mesh compressed in the vicinity of the adhesion layer. Finite element calculations have been performed using Comsol Multiphysics v. 5.3.0.316. An analysis of the problem solutions for a different number of finite elements has been carried out, the refinement of the partition grid does not lead to a significant improvement in the result, which indicates the results stability. The used number of FEs is 14000 with 1640 nodes. The number of degrees of freedom is 3280.
-duralumin, -bronze, -made from silver solder While solving the problem of optimization within the scheme of approximation of material based on (10), the following coefficients are fixed: a = 3 and b = 5. In the target function (15), we used p = 7 and q = 0.25. The values of these parameters were obtained as a result of numerical experiments for the problems under consideration. Figure 3 shows examples of the distribution of shear stresses along the upper and lower boundary of space Ω3 for the case of fully filled solder material (curves 1 and 4) and for the case of topologically optimal structure for γ = 0.6 (curves 2 and 3). In Figure 3, and later in the paper, the vertical axis shows the values of shear stresses 12Pa, whereas the horizontal axis shows the coordinate x1[mm] along the solder layer. From the distribution of curves, it can be concluded that in the case study under consideration, the optimization process implied a decrease of the peak values of shear stresses by more than 3 (4) times for the upper (lower) solder boundary. While solving the problem of optimization within the scheme of approximation of material based on (10), the following coefficients are fixed: a = 3 and b = 5. In the target function (15), we used p = 7 and q = 0.25. The values of these parameters were obtained as a result of numerical experiments for the problems under consideration. Figure 3 shows examples of the distribution of shear stresses along the upper and lower boundary of space Ω 3 for the case of fully filled solder material (curves 1 and 4) and for the case of topologically optimal structure for γ = 0.6 (curves 2 and 3). In Figure 3, and later in the paper, the vertical axis shows the values of shear stresses σ 12 Pa, whereas the horizontal axis shows the coordinate x 1 [mm] along the solder layer.
-duralumin, -bronze, -made from silver solder While solving the problem of optimization within the scheme of approximation of material based on (10), the following coefficients are fixed: a = 3 and b = 5. In the target function (15), we used p = 7 and q = 0.25. The values of these parameters were obtained as a result of numerical experiments for the problems under consideration. Figure 3 shows examples of the distribution of shear stresses along the upper and lower boundary of space Ω3 for the case of fully filled solder material (curves 1 and 4) and for the case of topologically optimal structure for γ = 0.6 (curves 2 and 3). In Figure 3, and later in the paper, the vertical axis shows the values of shear stresses 12Pa, whereas the horizontal axis shows the coordinate x1[mm] along the solder layer. From the distribution of curves, it can be concluded that in the case study under consideration, the optimization process implied a decrease of the peak values of shear stresses by more than 3 (4) times for the upper (lower) solder boundary. From the distribution of curves, it can be concluded that in the case study under consideration, the optimization process implied a decrease of the peak values of shear stresses by more than 3 (4) times for the upper (lower) solder boundary. Table 1 gives references to the figures on which the optimal distribution of the solder structure along the area Ω 3 for different values of γ, i.e., the amount of the solver solder is obtained. In Figure 4, and Figures 6, 8, 10, areas filled with solder with a Young's modulus value of E 3 are indicated in red, and voids are indicated in white. The transition colors correspond to the gradient properties of the solder with the Young's modulus varying from 0 to E 3 . The third column of the table gives the averaged value of shear stress along the solder joint area The percentage decrease in the peak values of shear stress in the optimal structure with respect to the peak values of shear stresses ∆ = σ Heonm 12 −σ onm 12 σ Heonm

12
· 100% of the input non-optimal construction is presented in the fourth column.
Based on the results shown in Table 1, the following conclusions were drawn. The smallest value of the averaged tangential stresses S in the solder layer is obtained for the minimum value of the solder amount γ. It is in agreement with physical expectations, i.e., the contact area of the first and third layer is reduced, and consequently the difference of their heating properties plays a less important role. In all cases the gradient of the optimal solder construction is exhibited, i.e., its non-homogeneity along both thickness and length.

Case Study 2
We consider a thermoelastic three-layer package free from external mechanical loads shown in Figure 5. In contrast to the previous example, areas Ω 1 and Ω 2 have different geometric dimensions as well as different sizes and forms of the solder area Ω 3 . Other parameters and material characteristics are selected in the same way as in case study 1.
Based on the results shown in Table 1, the following conclusions were drawn. The smallest value of the averaged tangential stresses S in the solder layer is obtained for the minimum value of the solder amount γ. It is in agreement with physical expectations, i.e., the contact area of the first and third layer is reduced, and consequently the difference of their heating properties plays a less important role. In all cases the gradient of the optimal solder construction is exhibited, i.e., its nonhomogeneity along both thickness and length.

Case Study 2
We consider a thermoelastic three-layer package free from external mechanical loads shown in Figure 5. In contrast to the previous example, areas Ω1 and Ω2 have different geometric dimensions as well as different sizes and forms of the solder area Ω3. Other parameters and material characteristics are selected in the same way as in case study 1. We investigated the topologically optimal microstructure of the solder with respect to length l of layer Ω3 (it is shown in the first column of Table 2). The second column presents the form of area Ω3 with the optimal structure of the solder material distribution (Figure 6).

Case Study 2
We consider a thermoelastic three-layer package free from external mechanical loads shown in Figure 5. In contrast to the previous example, areas Ω1 and Ω2 have different geometric dimensions as well as different sizes and forms of the solder area Ω3. Other parameters and material characteristics are selected in the same way as in case study 1. We investigated the topologically optimal microstructure of the solder with respect to length l of layer Ω3 (it is shown in the first column of Table 2). The second column presents the form of area Ω3 with the optimal structure of the solder material distribution (Figure 6).  We investigated the topologically optimal microstructure of the solder with respect to length l of layer Ω 3 (it is shown in the first column of Table 2). The second column presents the form of area Ω 3 with the optimal structure of the solder material distribution (Figure 6).   As can be seen from Table 2, the introduction of additional areas of the solder by extending the length of the layer does not give the expected additional decrease in the peak value of shear stresses. The fundamental solder structure in all considered cases can be found between the upper and bottom package layers while the remaining zones of the solder do not affect the level of tangential stresses and can be neglected. Therefore, the most optimal microstructure is the one which along the length of the solder layer is equal to length of the upper layer l = 2.5 mm. Figure 7 presents distributions of shear stresses of the non-optimal (curves 1, 4) and optimal (curves 2, 3) solder structure with respect to the upper and bottom layer boundary Ω3 for l = 2.5 mm. As can be seen in Figure 7, as a result of the optimization, a double decrease in the peak shear stresses The allowed values of the solder amount are shown in the third column of the Table 2: value γ = 1 corresponds to the initial non-optimal construction with fully filled solder area Ω 3 (γ = 0.75). In the fourth column, the average value of shear stress S follows from formula (14). In the fifth column, the percentage decrease of the peak shear stresses in the optimal structure with respect to the peak value of counterpart stresses in the input non-optimal construction is given (see Figure 7).
As can be seen from Table 2, the introduction of additional areas of the solder by extending the length of the layer does not give the expected additional decrease in the peak value of shear stresses. The fundamental solder structure in all considered cases can be found between the upper and bottom package layers while the remaining zones of the solder do not affect the level of tangential stresses and can be neglected. Therefore, the most optimal microstructure is the one which along the length of the solder layer is equal to length of the upper layer l = 2.5 mm. Figure 7 presents distributions of shear stresses of the non-optimal (curves 1, 4) and optimal (curves 2, 3) solder structure with respect to the upper and bottom layer boundary Ω 3 for l = 2.5 mm. As can be seen in Figure 7, as a result of the optimization, a double decrease in the peak shear stresses at the upper and bottom solder boundaries was achieved.
In addition to the solder in the form of a straight layer Ω 3 , we consider other possible geometric forms of the solder area. The first column of Table 3 gives data in which the input forms of area Ω 3 are compared with optimal gradients of the microstructures of the solder distribution ( Figure 8).
The second column of Table 3 shows the given area of space Ω 3 . The third column shows the values of the allowed amount of solder γ. The fourth column presents the values of shear stresses S computed using formula (16). The fifth column gives the percentage decrease in the peak values of shear stresses in the optimal microstructure with respect to the peak values of counterpart shear stresses in a non-optimal construction.
(c) -duralumin, -bronze, -made from silver solder, -lack of solder, -gradient solder  As can be seen from Table 2, the introduction of additional areas of the solder by extending the length of the layer does not give the expected additional decrease in the peak value of shear stresses. The fundamental solder structure in all considered cases can be found between the upper and bottom package layers while the remaining zones of the solder do not affect the level of tangential stresses and can be neglected. Therefore, the most optimal microstructure is the one which along the length of the solder layer is equal to length of the upper layer l = 2.5 mm. Figure 7 presents distributions of shear stresses of the non-optimal (curves 1, 4) and optimal (curves 2, 3) solder structure with respect to the upper and bottom layer boundary Ω3 for l = 2.5 mm. As can be seen in Figure 7, as a result of the optimization, a double decrease in the peak shear stresses at the upper and bottom solder boundaries was achieved.
In addition to the solder in the form of a straight layer Ω3, we consider other possible geometric forms of the solder area. The first column of Table 3 gives data in which the input forms of area Ω3 are compared with optimal gradients of the microstructures of the solder distribution ( Figure 8). The second column of Table 3 shows the given area of space Ω3. The third column shows the values of the allowed amount of solder γ. The fourth column presents the values of shear stresses S computed using formula (16). The fifth column gives the percentage decrease in the peak values of shear stresses      As can be seen from Table 3, the introduction of additional solder parts located outside the domain between the package layers did not result in an additional reduction of peak shear stresses. The fundamental solder structure is located between upper and bottom layers of the package, and the remaining zones of the solder do not play any important role and can be neglected. Therefore, the most optimal structure is the one which occurs in the third row of Table 3.

Topological Optimization of Non-Symmetric Packages under Thermal and Mechanical Loads: Numerical Results
The packages can be connected using soldering in different ways: by using a continuous solder, flange bevel solder, square groove solder, double butt solder, arc solder, fillet solder, etc. In the solder joints, the Stresses are distributed non-uniformly. Stress peaks occur in the points of geometric discontinuities. In general, the linkage by solder involves the occurrence of discontinuities at the ends of the solder. The occurring non-uniformities generate bending moments originating from eccentric loads and caused by non-uniform distribution of moments around the solder layer. These moments imply the occurrence of disruptive normal stresses in the solder layer. There are ways to reduce the harmful effects of the emerging eccentric loads in the solder layers. For instance, it is known that it is useful to apply the narrowing of edges of solder layers and introduce pre-bending. The decrease of the maximum shear stresses can be achieved by increasing the length and thickness of solder, the thickness of the combined layers and reducing the solder modulus. In the given example, all geometric parameters of the linkage remain constant, while a decrease in the maximum values of stresses at the ends of the solder is achieved by topological optimization of the solder structure.
Let us consider the construction shown in Figure 9. In contrast to the previous cases, despite thermal excitation, the mechanical load along the x 1 axis is introduced. The non-homogenous distribution of the moments around the welding layer results in the appearance of additional shear thermal stresses in the welding joint.
Materials 2020, 13, x FOR PEER REVIEW 13 of 18 the maximum shear stresses can be achieved by increasing the length and thickness of solder, the thickness of the combined layers and reducing the solder modulus. In the given example, all geometric parameters of the linkage remain constant, while a decrease in the maximum values of stresses at the ends of the solder is achieved by topological optimization of the solder structure. Let us consider the construction shown in Figure 9. In contrast to the previous cases, despite thermal excitation, the mechanical load along the x1 axis is introduced. The non-homogenous distribution of the moments around the welding layer results in the appearance of additional shear thermal stresses in the welding joint. The obtained results are included in Table 4. The first column presents input temperature θ [°C] The second column contains mechanical load F [N]. The third column refers to the figures, where the optimal weld microstructure distribution in space Ω3 is obtained. The average values of shear stresses S are shown in the fifth column, while the sixth column gives the decrease in peak values of shear stresses (in percent) in the optimal microstructure with respect to the peak values of the counterpart shear stresses in the input non-optimal construction. Table 4. Сomparing results of the topological optimization under thermal and mechanical loads with the original construction. The obtained results are included in Table 4. The first column presents input temperature θ[ • C] The second column contains mechanical load F[N]. The third column refers to the figures, where the optimal weld microstructure distribution in space Ω 3 is obtained. The average values of shear stresses S are shown in the fifth column, while the sixth column gives the decrease in peak values of shear stresses (in percent) in the optimal microstructure with respect to the peak values of the counterpart shear stresses in the input non-optimal construction.   Below, we present graphs of shear stresses in various cases of the mentioned loads: only the mechanical load is applied (Figure 11а) and only thermal load is applied (Figure 11b,c). Curves 2 and 3 of stress values 12 in all figures correspond to optimal constructions, while curves 1 and 4 correspond to the input (non-optimized) constructions. We begin with the optimization process of the weld microstructure subjected only to the mechanical load. In a given case, the averaged value S for both the optimized and input microstructure of the welding layer remains unchanged, which requires balancing of extension load F. However, the peak stress values are reduced by more than 4.5 Below, we present graphs of shear stresses in various cases of the mentioned loads: only the mechanical load is applied (Figure 11a) and only thermal load is applied (Figure 11b,c). Curves 2 and 3 of stress values σ 12 in all figures correspond to optimal constructions, while curves 1 and 4 correspond to the input (non-optimized) constructions. We begin with the optimization process of the weld microstructure subjected only to the mechanical load. In a given case, the averaged value S for both the optimized and input microstructure of the welding layer remains unchanged, which requires balancing of extension load F. However, the peak stress values are reduced by more than 4.5 (1.5) times in the lower (upper) boundary.
-duralumin, -bronze, -made from silver solder, -lack of solder, -gradient solder Below, we present graphs of shear stresses in various cases of the mentioned loads: only the mechanical load is applied (Figure 11а) and only thermal load is applied (Figure 11b,c). Curves 2 and 3 of stress values 12 in all figures correspond to optimal constructions, while curves 1 and 4 correspond to the input (non-optimized) constructions. We begin with the optimization process of the weld microstructure subjected only to the mechanical load. In a given case, the averaged value S for both the optimized and input microstructure of the welding layer remains unchanged, which requires balancing of extension load F. However, the peak stress values are reduced by more than 4.5 (1.5) times in the lower (upper) boundary.  The action of thermal load only (F = 0 N, θ = 100 • C, Figure 11b, and F = 0 N, θ = −100 • C, Figure 11c) implies a decrease of both S and peak values σ 12 by 4 (4.5) times on the upper (lower) boundary of the solder layer for θ = 100 • C and by 4.5 (4) times on the upper (lower) boundary of the solder layer for θ = −100 • C.
Simultaneous action of both mechanical and thermal load for F = 3000 N, θ = 100 • C (Figure 12a) on the upper (lower) boundary of the solder welding layer reduces the peak shear stresses σ 12 by 3.7 (3) times. In the case of F = 3000 N, θ = −100 • C ( Figure 12b) the counterpart values for the upper (lower) boundary correspond to the decrease of peak stresses by 4 (3) times.
This work continues a series of publications by authors dedicated to topological optimization [52][53][54]. Figure 11. Graphs of the shear stress distribution along upper (1, 2) and bottom (3,4) boundary of the solder layer under mechanical (a) and thermal (b,c) load.
The action of thermal load only (F = 0 N, θ = 100 °C, Figure 11b, and F = 0 N, θ = −100 °C, Figure   11c) implies a decrease of both S and peak values 12 by 4 (4.5) times on the upper (lower) boundary of the solder layer for θ = 100 °C and by 4.5 (4) times on the upper (lower) boundary of the solder layer for θ = −100 °C. Simultaneous action of both mechanical and thermal load for F = 3000 N, θ = 100 °C (Figure 12а) on the upper (lower) boundary of the solder welding layer reduces the peak shear stresses 12  by 3.7 (3) times. In the case of F = 3000 N, θ = −100 °C ( Figure 12b) the counterpart values for the upper (lower) boundary correspond to the decrease of peak stresses by 4 (3) times.
This work continues a series of publications by authors dedicated to topological optimization [52][53][54].

Concluding Remarks
Based on the investigations carried out, the following conclusions can be formulated. 1. A novel topological optimization methodology of structural soldered joints is developed. 2. The proposed methodology yielded stresses decrease for an optimal structure. 3. The main advantage of using topological optimization is that when designing structural elements, you do not need to know their geometry in advance. 4. The numerical implementation of the developed methodology for strain calculations is performed by the finite element method using the method of moving asymptotes. 5. In all cases, the constructed optimized structures demonstrate a significant decrease in the stress peaks in the solder layer at a smooth distribution of stresses along it, compared with a uniform solder layer. 6. The reliability is confirmed by the coincidence of numerical solutions obtained on the basis of the developed methodology for topological optimization with experimental data [6][7][8]. Peak stresses occur in small areas of the solder layer located near the connection ends; this coincides with the experimental results not reported here. 7. The developed methodology based on topological optimization can be used in the design of macro and micro structures.

Concluding Remarks
Based on the investigations carried out, the following conclusions can be formulated.

1.
A novel topological optimization methodology of structural soldered joints is developed.

2.
The proposed methodology yielded stresses decrease for an optimal structure. 3.
The main advantage of using topological optimization is that when designing structural elements, you do not need to know their geometry in advance.

4.
The numerical implementation of the developed methodology for strain calculations is performed by the finite element method using the method of moving asymptotes.

5.
In all cases, the constructed optimized structures demonstrate a significant decrease in the stress peaks in the solder layer at a smooth distribution of stresses along it, compared with a uniform solder layer.

6.
The reliability is confirmed by the coincidence of numerical solutions obtained on the basis of the developed methodology for topological optimization with experimental data [6][7][8]. Peak stresses occur in small areas of the solder layer located near the connection ends; this coincides with the experimental results not reported here. 7.
The developed methodology based on topological optimization can be used in the design of macro and micro structures.