Thermal deformation calculation method for multi-chip embedded board-level packaging

Abstract

The invention belongs to the technical field of chip packaging, and particularly relates to a thermal deformation calculation method for multi-chip embedded board-level packaging. On the basis of heat conduction and elastic mechanics theories, a three-dimensional finite difference thermal mechanical numerical model is constructed, and stress deformation generated by heating of a chip after the chip is embedded into a package and thermal expansion of a package material is described. The coupling problem of multi-chip heating is also considered. Compared with an analytical method, the finite difference method thermomechanical numerical model constructed by the method has higher universality during calculation, a complex analytical model does not need to be deduced in the model deduction process, and the time cost in the modeling process is lower. Compared with a finite element method, the finite difference method needs fewer computing resources, the configuration requirement for a computer is lower, and the time of the computing and solving process is shorter. And human resources and time cost can be saved under the condition that the calculation accuracy is ensured.

Description

Technical Field

[0001] The invention belongs to the technical field of chip packaging, and in particular relates to a thermal deformation calculation method for a multi-chip embedded board-level packaging. Background Art

[0002] As the demand for power density and miniaturization of electronic devices continues to increase, traditional single-chip packaging technology can no longer meet the requirements of high-performance electronic devices. Multi-chip FOPLP has become a development trend due to its excellent electro-thermal performance and cost-effectiveness, but it faces thermal management problems and mechanical stress problems in practical applications.

[0003] In the field of large-scale fan-out panel-level packaging (FOPLP) technology, thermal and mechanical stresses during switching can bring huge failure risks to integrated chips and packaging structures. Therefore, strict thermal and mechanical evaluation and design of these devices are crucial.

[0004] Numerical solutions for the temperature distribution can be obtained by a variety of methods. One approach is to directly apply the separation of variables technique to solve the Fourier equations, which is particularly effective for problems with well-defined boundary conditions. Another approach is to use the finite difference method (FDM), which provides a discrete solution that is well suited for complex geometries and non-uniform boundary conditions. The mechanical parameters can be obtained by both analytical and numerical methods. The transfer matrix method is used to analytically solve the thermodynamics of two-dimensional laminated beams [1], [2]. The coefficients in the general solution of the stress function are determined based on the Fourier sine series expansion of the boundary equations on the upper and lower surfaces of the beam and the linear relationship between the reaction force and displacement at the intermediate support of the beam [3]. This method has been used to analyze the stress and displacement distribution of multi-span beams with continuously varying thickness. The complexity of the model or boundary conditions greatly increases the difficulty of obtaining an analytical solution. The double Laplace transform technique (LTT) combined with the Adomia discretization method (ADM) has been shown to be an effective solution when thermal and mechanical loads are applied to the body. Specifically, the time reversal of LTT is achieved by applying the Durbin technique, which provides an accurate method for inverting the Laplace transform in the time domain [4]. Compared with thermal physics models, mechanical physics models are often more complex. On the one hand, the motion state and constraints of the mechanical model itself introduce a variety of modeling equations, such as the Euler-Bernoulli beam equation and Timoshenko beam theory that describe the vibration of the spacecraft. The Ritz method in the Newmark average acceleration method [5] can be used to perform dynamic analysis on an axial functionally graded material (FGM) cantilever beam subjected to sinusoidal dynamic loads. On the other hand, different mechanical properties usually require different theories. High-order shear deformation beam theory can be used for bending beams and free vibration beams. For the nonlinear free vibration of FGM beams, the local coordinated deformation assumption (LCDA) is proposed to solve the nonlinear frequency of FGM beams [6]. Finite element models are usually used for accurate solutions of complex problems [7]. It is mainly used to solve differential equations and has a variety of extended models. The scaled boundary finite element method (SBFEM) uses a semi-analytical method to solve differential equations and can discretize the model into a larger subdomain without losing accuracy [8]. The deflection analysis of composite beams can also be carried out using quasi-fitting finite element method with higher-order theory [9].

[0005] At present, the thermomechanical modeling methods are generally divided into analytical methods and numerical methods. The analytical method requires the derivation of the general solution of the Laplace equation and the constitutive equation. The solution process is complicated, and the general solution is different for different boundary conditions. After changing the boundary conditions, the equation needs to be re-derived and solved. Although the finite element method in the numerical method can solve the thermomechanical properties of complex models, it will consume a lot of computing time and resources. Finite difference does not require a complex differential function solution process, and the computing time and computing resources occupied are also lower than the finite element method. In order to ensure the accuracy of the calculation results and reduce the calculation cost, the finite difference method is very suitable for solving thermodynamic models. Summary of the invention

[0006] The object of the present invention is to provide a method for calculating thermal deformation of a multi-chip embedded board-level package that meets the market demand for miniaturized, high-power density electronic equipment.

[0007] The thermal deformation calculation method of multi-chip embedded panel-level package provided by the present invention first constructs a three-dimensional finite difference numerical calculation model that can accurately predict the thermal mechanical properties of multi-chip FOPLP, which is used to describe the stress deformation caused by the working heat generated by the chip after it is embedded in the package and the thermal expansion of the packaging material, including heat conduction analysis modeling and finite difference method thermoelastic numerical modeling; the specific steps are as follows:

[0008] (I) Heat conduction analysis modeling

[0009] For FOPLP structure, the packaging process usually involves stacking whole layers of materials. After the chip is fixed on the carrier, a whole piece of plastic encapsulation material is pressed on the chip, and then a whole copper layer is pressed on the plastic encapsulation material, and blind holes are drilled on the copper layer. The blind holes are used to achieve electrical connection between the chip and the copper layer. In order to distinguish the electrode area on the copper layer, a specific pattern can be etched on the copper layer. There is also a whole piece of plastic encapsulation material pressed on the copper layer and a whole piece of heat sink pressed on the copper layer. Therefore, it can be simplified to a rectangular structure composed of multiple stacked materials, such as Figure 1 As shown. The chip in the middle layer is the heating surface (the small square in the figure). This plane is selected as the three-axis origin plane in the heat conduction calculation, and the Fourier heat conduction equation of its three-dimensional structure is reduced to the Laplace equation, which is:

[0010]

[0011] where T(x,y,z) is the temperature at point (x,y,z). The general solution to the Laplace equation is usually expressed as θ(x,y,z):

[0012] θ(x,y,z)=T(x,y,z)-T air ;

[0013] T air is the ambient temperature; specifically, the general solution can be expressed as:

[0014]

[0015]

[0016] Where λ = mπ / L x ,m=1,2,3…,δ=nπ / L y ,n=1,2,3…, B i =φ(ξ)A i,i=1,2,3; The heat flux of the Q chip; for the convenience of calculation, φ(ξ) is used to represent B i and A i ξ can be replaced by λ, δ or β. x is the length of the package, L y is the width; considering the package thickness (e.g. 0.78 mm), which is much smaller than the cross-sectional length (e.g. 18 mm), the convective heat transfer from the sidewalls to the external environment can be ignored. Therefore, the heat dissipation to the environment is mainly carried out through the top and bottom surfaces. At the interface between different layers, the temperature and heat flux remain in a continuous state. The coefficient of each parameter is interrelated. It depends on the thermal power output by the chip, the thermal conductivity of the constituent materials, and the thickness of each structural component. Figure 1 The boundary conditions and their corresponding mathematical expressions are also given to complete the heat conduction model, specifically:

[0017] (1) The four walls are adiabatic boundaries:

[0018]

[0019] The upper and lower surfaces conduct convection heat exchange with the air:

[0020]

[0021] Among them, t p is the distance from the chip to the upper surface, t n is the distance from the chip to the bottom surface.

[0022] In the area where the chip generates heat:

[0023]

[0024] Among them, k p is the thermal conductivity of the upper layer, k n is the thermal conductivity of the lower layer, a, b are the length and width of the chip, and Q chip is the heat flux.

[0025] Outside the chip's heat generating area, the total amount of heat conducted in the up and down directions is the total amount of heat dissipated from the surface:

[0026]

[0027] Here, A c is the entire package cross-sectional area, A s is the cross-sectional area of ​​the chip, and h is the convection coefficient between the surface and the air.

[0028] Continuity conditions between different layers, temperature continuity, temperature gradient continuity, displacement continuity:

[0029] θi (x,y,z i )=θ i+1 (x,y,z i ),#(9)

[0030]

[0031] For multiple heat sources, the linear temperature rise superposition results in a composite heating effect for the entire system. This phenomenon enables the discrete heat dissipation of each heat source to be calculated individually and subsequently aggregated to determine the collective temperature distribution. Assuming there are n heat sources within the structure, the temperature distribution at different locations can be expressed as:

[0032] T(x,y,z)=θ 1 (x,y,z)+θ 2 (x,y,z)+…+θ n (x,y,z)+T air ,#(11)

[0033] (II) Finite difference method for thermoelastic numerical modeling

[0034] When the temperature changes, the thermal expansion of the material affects the deformation and stress distribution within the structure.

[0035] For thermoelasticity, the constitutive equation can be expressed in tensor form as:

[0036]

[0037] Among them, ε ij represents the sum of elastic strain and thermal strain, i,j = 1, 2, 3. G is the shear modulus, ν is the Poisson's ratio, T 0 is the initial temperature, T is the actual temperature, (TT 0 ) represents the temperature change, α is the linear thermal expansion coefficient, δ ij is the Kronecker symbol, σ kk and σ ij is the tensor form of stress.

[0038] The thermoelastic constitutive equation expressed in terms of stress tensor is:

[0039] σ ij =2με ij +[λε kk -α(3λ+2μ)(TT 0 )]δ ij #(13)

[0040] Where μ = G, λ is the Lame constant, defined as ε kk is the tensor representation of the total strain.

[0041] The geometric equation in tensor form is given by:

[0042]

[0043] Among them, u i,j is the tensor form of displacement. Without considering external forces, the tensor form of equilibrium equation is expressed as:

[0044] σ ij,j =0#(15)

[0045] Substituting the geometric equation (14) into the constitutive equation (12) and then into the equilibrium equation (15), the differential expressions of the three displacement components are obtained:

[0046]

[0047]

[0048] First, the finite difference method is used to solve the second-order partial differential equation of the three-dimensional displacement (u, v, w). Assuming that the bottom surface is fixed on the externally connected circuit board and there is no displacement change on the bottom surface, the expression of the second-order partial differential equation of the three-dimensional displacement (u, v, w) solved by the finite difference method is as follows, corresponding to equation (16):

[0049]

[0050] The corresponding two equations (17) and (18) are:

[0051]

[0052] and

[0053]

[0054] Δx, Δy, Δz are the step sizes of the grid in the X, Y, and Z directions respectively.

[0055] That is, the partial differential equations are converted into the corresponding difference equations to obtain a linear algebraic system, which can be expressed as:

[0056] Au=b,

[0057] Where A represents the coefficient matrix, whose elements are all the coefficients in the above formula, including b represents a constant term, the elements of which include

[0058] (III) Solution of deformation and stress

[0059] The coefficient matrix A is a large sparse matrix. The dimension of the matrix is ​​3 dots. all×3dots all .dots all represents the total number of all grid points in the discretized structure. Since each point has three degrees of freedom (u, v, w), the matrix dimension is multiplied by 3. In order to solve all unknown variables in the three equations in one step, a comprehensive parameter assignment of the sparse matrix is ​​required (i.e., assigning values ​​to the matrix A). Here, the linear indexing of the matrix is ​​converted to the subscript to assign values ​​to the sparse matrix.

[0060] The solution for the displacement vector (u,v,w) can be obtained by solving the system:

[0061] u=A\b.

[0062] Substituting this solution into the geometric equation (14) will determine the strain distribution at each point in the structure. At the same time, by applying the constitutive equation (13), the stress distribution of the entire structure can also be obtained.

[0063] The present invention can perform rapid and accurate simulation calculations on the thermomechanical characteristics of board-level packaging. The difference between the temperature predicted by the numerical model of the present invention and the temperature obtained by finite element analysis is within 1%, and the difference in mechanical response is within 5%. The thermomechanical calculation of board-level packaging by the method proposed by the present invention can save time on model building and computing resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] Figure 1 Schematic diagram of the FOPLP thermomechanical model.

[0065] Figure 2 Temperature distribution results of FOPLP numerical thermomechanical model. (a) upper surface, (b) chip layer, (c) lower surface; FOPLP finite element simulation temperature distribution results, (d) upper surface, (e) chip layer, (f) lower surface; comparison of results of two methods, (g) upper surface, (h) chip layer, (i) lower surface.

[0066] Figure 3 The deformation distribution of the upper surface in the Z direction predicted by the thermoelastic numerical model. Among them, (a) the finite difference numerical method proposed by the present invention; (b) the finite element simulation; (c) the comparison of these two thermoelastic prediction methods. DETAILED DESCRIPTION

[0067] The method proposed in the present invention is compared with the simulation results of finite element software.

[0068] The case uses a silicon carbide board-level packaging device with a package size of 20mm*20mm. The structure and thickness of the device from top to bottom are: solder mask layer 20um, RDL layer 70um, blind hole layer 80um, chip layer 200um, solder layer 20um, carrier layer 230um, plastic layer 130um, heat sink 30um. Assume that during operation, the upper and lower bridge arms work alternately with a duty cycle of 50%, the switching frequency is 100kHz, and the on-current is constant at 10A. According to the chip data sheet, the on-resistance of a single chip is 40mΩ. Without considering the switching loss, it is calculated that the thermal power of each chip is 5W. Under these conditions, the temperature distribution of the upper and lower surfaces and the chip layer in the structure is as follows Figure 2 shown. Figure 2 (a)(b)(c) are the temperature distribution of the upper surface, chip layer and lower surface obtained by the calculation method proposed in the present invention; (d)(e)(f) are the temperature distribution of the upper surface, chip layer and lower surface obtained by the calculation method proposed by the finite element simulation software; (g)(h)(i) are the data comparison and error of the temperature distribution of different layers of the two methods. It can be seen that the heat is concentrated at the location of each layer of the chip, and the temperature of the chip layer is the highest. In addition, the temperature variation within the same layer is relatively small. The numerical thermal model effectively captures the temperature distribution characteristics of FOPLP. We conducted a comparative analysis along the longitudinal section of the center of the chip, and the results showed that the difference between the temperature predicted by the numerical model of the present invention and the temperature obtained by the finite element analysis is within 1%. This confirms the accuracy of the thermal model proposed in the present invention.

[0069] Figure 3 The deformation distribution of the upper surface in the Z direction calculated by the thermoelastic model is shown in FIG. Figure 3 (a) is the deformation distribution in the Z direction obtained by the calculation method proposed in the present invention; (b) is the deformation distribution in the Z direction of the lower surface obtained by the calculation method proposed by the finite element simulation software; (c) is the data comparison and error of the deformation distribution in the Z direction of the lower surface of the two methods. The thermomechanical model based on the finite difference method proposed in the present invention accurately reflects the mechanical response caused by the internal heating of the FOPLP chip. Figure 3 (c) It can be seen that the error between the proposed method and the finite element simulation is within 5%. The main source of error is the mesh discretization at the edges. Increasing the mesh resolution can significantly improve the accuracy of the numerical results and provide a more accurate approximation of the true value.

[0070] References:

[0071] [1]Z.Zhang,D.Zhou,H.Fang,J.Zhang,and X.Li,“Analysis of laminatedbeams with temperature-dependent material properties subjected to thermal andmechanical loads,”Composite Structures,vol.227,p.111304,Nov.2019,doi:10.1016 / j.compstruct.2019.111304.

[0072] [2]Y.Xu and D.Zhou,“Two-dimensional thermoelastic analysis of beamswith variable thickness subjected to thermo-mechanical loads,”AppliedMathematical Modelling,vol.36,no.12,pp.5818–5829,Dec.2012,doi:10.1016 / j.apm.2012.01.048.

[0073] [3]Y.Xu and D.Zhou,“Elasticity solution of multi-span beams withvariable thickness under static loads,”Applied Mathematical Modelling,vol.33,no.7,pp.2951–2966,Jul.2009,doi:10.1016 / j.apm.2008.10.027.

[0074] [4]J.N.Sharma and R.Kaur,“Flexural response of thermoelastic thinbeam resonators due to thermal and mechanical loads,”International Journal ofMechanical Sciences,vol.101–102,pp.170–179,Oct.2015,doi:10.1016 / j.ijmecsci.2015.07.014.

[0075] [5] D. “Forced Vibration Responses of Axially FunctionallyGraded Beams by using Ritz Method,”Journal of Applied and ComputationalMechanics,vol.7,no.1,pp.109–115,Jan.2021,doi:10.22055 / jacm.2020.34865.2491.

[0076] [6]K.Xie,Y.Wang,X.Fan,and T.Fu,“Nonlinear free vibration analysis offunctionally graded beams by using different shear deformation theories,”Applied Mathematical Modelling,vol.77,pp.1860–1880,Jan.2020,doi:10.1016 / j.apm.2019.09.024.

[0077] [7]G.Lin,P.Zhang,J.Liu,and J.Li,“Analysis of laminated composite andsandwich plates based on the scaled boundary finite element method,”CompositeStructures,vol.187,pp.579–592,Mar.2018,doi:10.1016 / j.compstruct.2017.11.001.

[0078] [8]N.Garg,A.Saputra,M.Donough,C.Song,A.W.Phillips,and B.G.Prusty,“Application of scaled boundary finite element method for three-dimensionalmodeling of bi-axial bending in thick laminated composite plates,”Mechanicsof Advanced Materials and Structures,Dec.2022,Accessed:Aug.28,2024.[Online].Available:https: / / www.tandfonline.com / doi / abs / 10.1080 / 15376494.2021.1989526

[0079] G.Shi,K.Y.Lam,and T.E.Tay,“On efficient finite element modeling ofcomposite beams and plates using higher-order theories and an accuratecomposite beam element,”Composite Structures,vol.41,no.2,pp.159–165,Feb.1998,doi:10.1016 / S0263-8223(98)00050-6。

Claims

1. A method for calculating thermal deformation of a multi-chip embedded board-level package, characterized in that: Firstly, a three-dimensional finite difference numerical calculation model is constructed to predict the thermomechanical properties of multi-chip FOPLP. It is used to describe the stress and deformation caused by the heat generated by the chip after it is buried in the package and the thermal expansion of the packaging material. It includes heat conduction analysis modeling and finite difference thermoelastic numerical modeling. The specific steps are as follows: (I) Heat conduction analysis modeling For the FOPLP structure, the packaging process involves stacking a whole layer of materials, which is simplified into a rectangular structure consisting of multiple stacked materials. The heating surface of the chip is selected as the axis origin plane in the heat conduction calculation. The Fourier heat conduction equation of its three-dimensional structure is reduced to the Laplace equation, which is: Where T(x,y,z) is the temperature at point (x,y,z); the general solution of the Laplace equation θ(x,y,z) is: Where λ = mπ / L x ,m=1,2,3…,δ=nπ / L y ,n=1,2,3…, B i =φ(ξ)A i ,i=1,2,3; φ(ξ) represents B i and A i ξ can be replaced by λ, δ or β; L x is the length of the package, L y is the width of the package; considering that the package thickness is much smaller than the cross-sectional length, the convective heat transfer from the sidewalls to the external environment is negligible; therefore, the heat dissipation to the environment is mainly carried out through the top and bottom surfaces; at the interface between different layers, the temperature and heat flux remain in a continuous state; the coefficient of each parameter is interrelated, and it depends on the thermal power output by the chip, the thermal conductivity of the constituent materials, and the thickness of each structural component; For multiple heat sources, the linear temperature rise superposition leads to a composite heating effect of the entire system; this phenomenon enables the discrete heat dissipation of each heat source to be calculated individually and subsequently aggregated to determine the collective temperature distribution; assuming there are n heat sources within the structure, the temperature distribution at different locations is expressed as: T(x,y,z)=θ1(x,y,z)+θ2(x,y,z)+…+θ n (x,y,z)+T air #(3) Among them, T air is the ambient temperature; (II) Finite difference method for thermoelastic numerical modeling When the temperature changes, the thermal expansion of the material affects the deformation and stress distribution within the structure; For thermoelasticity, the constitutive equation is expressed in tensor form as: Among them, ε ij represents the sum of elastic strain and thermal strain, i, j = 1, 2, 3; G is the shear modulus, v is Poisson's ratio, T0 is the initial temperature, T is the actual temperature, (T-T0) represents the temperature change, α is the linear thermal expansion coefficient, δ ij is the Kronecker symbol; σ kk and σ ij is the tensor form of stress; The thermoelastic constitutive equation expressed in terms of stress tensor is: s ij =2with ij +[le kk -α(3λ+2μ)(T-T0)]δ ij #(5) Where μ = G, λ is the Lame constant, defined as ε kk is the total strain tensor; The geometric equation in tensor form is given by: Among them, u i,j is the tensor form of displacement; without considering external forces, the tensor form of the equilibrium equation is expressed as: s ij,j =0#(7) Substituting the geometric equation (6) into the constitutive equation (4) and then into the equilibrium equation (7), we obtain the differential expressions of the three displacement components: First, the finite difference method is used to solve the second-order partial differential equation of the three-dimensional displacement (u, v, w); assuming that the bottom surface is fixed on the externally connected circuit board and there is no displacement change on the bottom surface, the expression of the second-order partial differential equation of the three-dimensional displacement (u, v, w) solved by the finite difference method is as follows, corresponding to equation (8): The corresponding two equations (9) and (10) are: and Δx, Δy, Δz are the step sizes of the grid in the X, Y, and Z directions respectively; That is, the partial differential equations are converted into the corresponding difference equations to obtain a linear algebraic system, which can be expressed as: Au=b, Where A represents the coefficient matrix, whose elements are all the coefficients in the above formula, including b represents a constant term, the elements of which include (III) Solution of deformation and stress The coefficient matrix A is a large sparse matrix with a dimension of 3 dots. all ×3dots all .dots all represents the total number of all grid points in the discretized structure; each point has three degrees of freedom (u, v, w), and the matrix dimension is multiplied by 3; in order to solve all unknown variables in the three equations in one step, a comprehensive parameter assignment A of the sparse matrix is ​​required, and here the linear index of the matrix is ​​converted to assign values ​​to the sparse matrix; The solution for the displacement vector (u,v,w) is obtained by solving the system: u=A\b; Substituting this solution into the geometric equation (6) will determine the strain distribution at each point in the structure; by applying the constitutive equation (5), the stress distribution of the entire structure is obtained.

Images