CN111723486B

CN111723486B - Double-sided PCB structure steady-state thermal analysis method

Info

Publication number: CN111723486B
Application number: CN202010572678.2A
Authority: CN
Inventors: 张亚斌; 董秀成; 吴昌东; 郭嘉; 张彼德; 陈琦
Original assignee: Xihua University
Current assignee: Xihua University
Priority date: 2020-06-22
Filing date: 2020-06-22
Publication date: 2021-05-04
Anticipated expiration: 2040-06-22
Also published as: CN111723486A

Abstract

The invention discloses a steady-state thermal analysis method for a double-sided PCB structure, which comprises the following steps: s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and obtaining an insulating layer temperature analytic solution matrix equation; s2, correcting the boundary condition of the substrate insulating layer, and taking the heating rate of the component as the known heat source condition to be recorded into a matrix equation expression of an analytic solution; s3, approximating a heat balance equation of the discrete metal layer based on a finite volume method, and summarizing to obtain a metal layer heat diffusion lumped matrix equation; s4, constructing a thermal analysis coupling equation set of the double-sided PCB structure; and S5, calculating and inducing to obtain a lumped matrix equation describing the linear relation between the unit potentials of the metal layer, calculating the Joule heating distribution into the metal layer thermal diffusion matrix equation, and performing iterative calculation with the temperature distribution to obtain the temperature distribution of the surface of each layer of the PCB under the premise of calculating the Joule heating of the circuit. And S6, verifying the calculation accuracy of the thermal analysis method of the double-sided PCB structure by adopting a COMSOL modeling operation result.

Description

Double-sided PCB structure steady-state thermal analysis method

Technical Field

The invention belongs to the technical field of PCB thermal analysis, and particularly relates to a steady-state thermal analysis method for a double-sided PCB structure.

Background

The PCB is a main carrier of electronic components, has become an important component of electronic equipment, and the reliability of the PCB is related to the reliable operation of the whole electronic equipment. One of the main factors affecting the reliable operation of the PCB and the components thereon is the temperature. The design of the PCB structure ensures that the heating of the components on the PCB structure can be effectively diffused and can exchange heat with the external environment, namely, the temperature rise of the components is ensured to be within a safe range, and the PCB also becomes a main medium for heat dissipation of the components. Therefore, under the condition of considering the actual PCB structure, the circuit layout and the thermal boundary condition, the thermal analysis of the PCB can realize the effect of assisting in optimizing the PCB design and analyzing the reliability of the electronic system.

The PCB thermal simulation software widely applied in the market at present is FlothERM of Mentor company and Icepak of ANSYS company. The former is based mainly on approximate equivalent thermal conductivity related to the local copper coverage of the PCB in the algorithm principle, and the latter is based on the finite element method. The COMSOL and Sigrity PowerDC based on the finite element method can also perform thermal simulation analysis on the basis of completing PCB and device modeling.

Byron Blackmore of Mentor corporation in 2009 disclosed a method adopted by FloTHERM to derive local equivalent thermal conductivity from the local copper coverage of PCBs. According to the data announced by the company website, the occupation rate of Flotherm in the PCB thermal analysis market in the year is up to 88%. In the past decade, Flotherm has also been used by many national scholars to develop relevant application studies. However, a white paper disclosed in 2018 by Mentor corporation reveals that the method based on local copper coating rate adopted previously may have the problems of overestimating local thermal conductivity and underestimating actual temperature rise. The white paper indicates that another empirical algorithm is added in the Flotherm XT v3.2 version of software to analyze the thermal distribution of the PCB, the new algorithm is a solution method based on the local transverse thermal conductivity summarized after the analysis of 140 PCB products, and the thermal conductivity is still related to the copper coating rate; however, the white paper also indicates that empirical algorithms may underestimate actual thermal conductivity and thus overestimate actual temperature rise. Furthermore, Flotherm is expensive to use, and has a key authorization (Lisence) cost of 60 to 200 million RMB in China.

Finite Element Method (FEM) is a numerical discrete method widely used, and is characterized in that a continuous unknown function is approximated to a finite sub-function (generally in a polynomial form containing unknown variables which can be approximated to linear change in the discrete domain) effective in the discrete control domain in a discrete finite space or time, the finite sub-functions are substituted into an original continuous equation to establish the relation between corresponding variables of each discrete sub-domain, and then the approximate distribution of the unknown variables can be solved by using a given boundary condition. The software such as ANSYS Icepak and COMSOL Multiphysics based algorithms are used by many scholars to perform thermal simulation analysis of power devices, PCBs, and electronic systems.

The results of questionnaires on the use experiences of multi-physical field analysis software such as ANSYS, COMSOL and the like in the industry and the academia of the university of Necklands, Lester university, UK, Song Xueguan and the like of the university of the big connective workers in combination with the university of Necklands, show that more than 60% of users are satisfied with the operation accuracy of the software, but the same more than 60% of users consider that the operation speed of the software is slow. The CFD software based on the finite element method generally needs to disperse the whole structure to be analyzed in the simulation operation process, which may be one of the main reasons for the unsatisfactory operation speed. Meanwhile, the price of the software is expensive, taking COMSOL as an example, the price of the software is between 10 and 40 thousands according to the selected functional module.

From the development trend of the thermal analysis method of the PCB, the method for establishing the overall transverse and longitudinal equivalent thermal conductivity of the PCB in the early stage cannot be used for thermal analysis of the PCB as a general method due to poor structural applicability, otherwise, larger calculation errors can be brought. Therefore, methods for deriving equivalent thermal conductivity and finding local discretization orthotropic thermal conductivity based on the local copper cladding rate of each layer of the PCB are developed subsequently, and common characteristics of the methods are that errors caused by using the local equivalent thermal conductivity are expected to be reduced by discretizing the PCB structure and enhancing analysis of the local line distribution and the laminated board structure of the PCB.

Disclosure of Invention

The invention aims to provide a steady-state thermal analysis method of a double-sided PCB structure aiming at the defects in the prior art so as to solve the problem of large thermal analysis error of the prior PCB structure.

In order to achieve the purpose, the invention adopts the technical scheme that:

a steady-state thermal analysis method for a double-sided PCB structure comprises the following steps:

s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation;

s2, calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the boundary condition of the substrate insulating layer, and calculating the heating rate of the component as the known heat source condition into a matrix equation expression of an analytic solution;

s3, approximating a heat balance equation of a discrete metal layer based on a finite volume method, summarizing to obtain a metal layer heat diffusion lumped matrix equation, and counting the heat conduction between the metal layer and a surface heating device thereof and the heat conduction of the metal via hole;

s4, constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the corrected insulating layer temperature analytic solution matrix equation and the metal layer thermal diffusion matrix equation based on numerical value dispersion;

s5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and calculating potential distribution in the metal layer, so as to further calculate current density distribution and Joule heating distribution of the metal layer, and calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with temperature distribution to calculate temperature distribution of the surface where each layer of the PCB is positioned under the condition of Joule heating of a circuit;

and S6, verifying the calculation accuracy of the thermal analysis method of the double-sided PCB structure by adopting a COMSOL modeling operation result.

Preferably, step S1 is to construct a steady-state thermal balance laplace equation of the double-sided PCB insulating layer, and solve a fourier series analytic solution coefficient based on the boundary condition of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation, including:

constructing a stable heat balance Laplace equation of the double-sided PCB insulating layer, and obtaining boundary conditions of the periphery, the upper surface and the lower surface of the PCB insulating layer:

wherein T is the amount of temperature change from ambient temperature or the temperature at ambient temperature of 0 ℃, q_iu(x, y) and q_id(x, y) is a distribution function of the heat flux density introduced by the upper and lower surfaces, h_uAnd h_dHeat transfer coefficients of the upper and lower surfaces of the insulating layer, k_iIs the thermal conductivity of the insulating layer, L_x、L_y、L_inThe lengths of the insulating layer in the x direction and the y direction and the thickness of the insulating layer in the z direction are respectively measured under a Cartesian coordinate system;

the expression form of Fourier series analytic solution corresponding to the steady-state heat balance equation is as follows:

wherein, C₁，C₂，C_n,m，C_γn,mIs a Fourier series coefficient, beta_n、μ_m、γ_n,mCharacteristic values in a Fourier series;

decomposing a heat balance equation and a boundary condition by adopting a heat effect superposition principle:

wherein theta and eta are superposition component variables of the insulating layer heat balance equation solution;

substituting Fourier series solution into lower surface z corresponding to theta as L_inThe relation between partial coefficient and characteristic value in solution is obtained under the condition of thermal boundary, and when h is reached_dWhen not zero, C₂And C_γn,mThe corresponding expression is:

wherein H_uAnd H_dDeriving intermediate parameters of the process for the analytic solution; when h is generated_dIs zero, C₁,C_γn,mComprises the following steps:

substituting the series solution into the thermal boundary condition of the upper surface to obtain:

multiplying both sides of the equation by cos (. beta.)_nx)cos(μ_my) and the integral is carried out in the surface heating area, according to the characteristics of the trigonometric function, the integral after multiplying different series numbers is zero, so that the left side of the equation only contains cos²(β_nx)cos²(μ_my) is not zero after integration, and the right side of the equation is used for each heating unit q in the heating area_iu,iIntegrating and then summing, and calculating to obtain a coefficient C₁And C_n,m：

Wherein d is_qIntegration Domain Sq for grid cell Length after discretization of surface Heat Source grid_iu,iI.e. representing an incoming heat flow density of q_iu,iSurface discrete unit area of δ_nAnd delta_mTo resolve dimensionless intermediate parameters of the derivation process, q_iu,1～q_iu,N1The unit heat flux density is obtained after discretization of the heat source grid on the surface of the insulating layer; if the heat flux density transmitted by a certain heat source is more uniform, grid dispersion is not carried out on the longitudinal heat flux density transmitted by the area where the heat source is located; further, the summation operation of Fourier series can be expressed by means of vector operation, and for the above analytic solution, q_iu,1～q_iu,N1I.e. the column vector q of the heat flow density is constructed_iuThen, the analytic solution of any point in the single-sided insulating layer can be expressed as:

solving the temperature distribution of the designated area is the thermal resistance row vector R_u(x,y,z)Expanding into a thermal resistance matrix R according to the coordinates of the corresponding distribution points of the region_uWherein each row is corresponding to a thermal resistance row vector for solving the temperature of a certain point, and then R is added_uMultiplying the temperature distribution vector by the corresponding heat flow density vector q to obtain a temperature distribution vector of the specified area, namely the heat effect generated by the corresponding q in the specified area; finally, superposing Fourier series analytic solutions in a matrix form corresponding to theta and eta of the specified region to obtain the temperature distribution of the region;

thus, the matrix equation corresponding to the analytic solution temperature distribution vector of the upper and lower surfaces of the insulating layer in the double-sided PCB will have the following form:

wherein, the distribution vectors of the heat flux density transmitted into the insulating layer by the upper and lower surface metal layers are q_M,u、q_M,dWherein q is_uAnd q is_dThe distribution vector q of heat flux density transferred between the upper and lower surface heating components and the covering region thereof_uComprising q_iu(x, y) heat flux density distribution after dispersion, q_dComprising q_id(x, y) a heat flow density distribution after discretization; r_Mu,uIs q_M,uAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Mu,dIs q_M,dAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Md,uIs q_M,uAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Md,dq_M,dAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_iu,uIs q_uAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_iu,dIs q_dAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_id,uIs q_uAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_id,dIs q_dAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; due to q_uAnd q is_dThe heat flow transmitted through the surface of the metal layer is counted in the heat diffusion matrix equation corresponding to the metal layer, and only q is counted in the analytic solution matrix equation_iu(x, y) and q_id(x, y) heat flow transferred between the associated top and bottom surface heat generating components and the surface area of the insulating layer directly covered thereby, R_iu,u、R_iu,d、R_id,uAnd R_id,dWill include a partial all zero column vector.

Preferably, step S2 includes calculating the thermal effect of the device based on the thermal resistance parameter of the device, and modifying the thermal boundary condition of the substrate insulating layer to obtain a modified insulating layer temperature analytic solution matrix equation expression, which includes:

setting the heat flux density vector of the covering area of the upper surface heating device as q_h,u：

q_h,u＝q_u+h_uT_J,u

Wherein, T_J,uA temperature distribution vector for an area occupied by the upper surface device; because the covered area of the component does not exist direct connection with the external environmentIs in h_uThe calculated heat exchange under the condition, but the part of the heat exchange is already taken into account in the analytic solution derivation process of the step S1, so the occupied area of the component is required to be determined according to h_uCompensating for the multiple calculated heat exchanges; in addition, the thermal resistance parameters provided by the components can be used for deducing the heat transfer between the related components and the PCB surface, and the specific deduction process is as follows:

the two formulas are combined to obtain:

q_h,u＝C_P,uq_J,u-C_T,uT_J,u

wherein the temperature distribution of the upper surface element device surface is T_topThe heat transfer coefficient of the surface to the outside is h_topThe equivalent heat conduction corresponding to the heat sink existing on the surface thereof is also taken into account_top(ii) a Mean temperature of the core is T_j(ii) a The average thermal resistance between the surface and the core is psi_jt(ii) a The average thermal resistance between the junction center and the bottom of the device is psi_jc(ii) a Equivalent thermal resistance of solder mask coating is psi_c,u(ii) a The heat flux density transferred from the core to the surface of the component is q_jt(ii) a The heating volume density vector corresponding to the upper surface element device is P_J,u(ii) a The diagonal matrix composed of height parameters of the elements is C_J,u(ii) a The density vector of the heating surface corresponding to the upper surface element device is q_J,u；C_P,uFor taking into account the thermal resistance parameter of the device_J,uThe coefficient of (d); c_T,uFor taking into account the corresponding T after the thermal resistance parameter of the device_J,uThe coefficient of (a);

with q_h,uAlternative upper surface temperature distribution analytical solution T_uQ in (1)_u：

Wherein, T_J,dThe temperature distribution vector of the coverage area of the lower surface heating device is obtained; r_iu,uC_P,uq_J,uThe term is a vector of known constants, and T_uIncluding T_J,uAll elements of the vector, so the matrix coefficient can be transformed by adding zero to the ascending dimension to realize the T in the matrix equation_uInstead of T_J,uTo achieve the goal of unifying vectors, the equivalent transformation equation is as follows:

C_Ju,uT_u＝R_iu,uC_T,uT_J,u

wherein, C_Ju,uTo realize T_J,uAnd T_uAnd (3) equivalently transforming the coefficients of the related terms to further obtain:

T_u＝R_Mu,uq_M,u+R_Md,uq_M,d+R_id,uq_d+R_iu,uC_P,uq_J,u-C_Ju,uT_u

similarly, when there is a device on the lower surface, the heat flow transferred between the device on the lower surface and the surface unit, which takes into account the analytical solution, is also corrected, i.e. q in the above formula_dItem replacement:

wherein the heat flux density vector of the compensated lower surface covered by the heating device is q_h,d(ii) a The density vector of the heating surface corresponding to the lower surface component is q_J,d；C_P,dFor taking into account the thermal resistance parameter of the device_J,dThe coefficient of (a); and C_Ju,uFunction of (C) is similar to the definition_Jd,uFor realizing the temperature distribution T of the coverage area of the lower surface heating device_J,dAnd T_dCoefficients of the correlation term equivalent transform;

in the same way, the T corrected by the thermal boundary condition can be obtained_dCorresponding analytic solution matrix equation, and after correction, T_uAnd T_dCorresponding coupletThe analytical solution matrix equation becomes:

wherein, with C_Ju,uAnd C_Jd,uFunction of (C) is similar to the definition_Ju,d、C_Jd,dRespectively for realizing T in the lower surface analytical solution matrix equation_J,uAnd T_uRelated item and T_J,dAnd T_dThe correlation term is equivalent to the transformed coefficients.

Preferably, the step S3 approximates the discrete metal layer heat balance equation based on the finite volume method, generalizes the metal layer heat diffusion matrix equation, and includes the heat conduction between the metal layer and the surface heat generating device thereof, and the heat conduction of the metal via, and includes:

dispersing the metal layer by adopting a structured square grid, and dispersing a heat balance equation of heat conduction in the metal layer based on second-order Taylor series approximation to obtain an upper surface discrete metal layer heat balance equation:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_c,uThe joule heating rate distribution per unit volume of the upper surface metal layer; t is_c,uIs the upper surface metal cell temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_u,uIs T_c,uThe longitudinal heat flow density of the upper surface of the corresponding unit can be the density of heat flow transferred between the corresponding unit and the heating device; q. q.s_d,uIs T_c,uThe longitudinal heat flux density transmitted between the lower surface of the corresponding unit and the insulating layer; t is_E,uIs T_c,uThe temperature of the metal unit adjacent to the east side of the corresponding unit; t is_W,uIs T_c,uThe temperature of the corresponding unit west side adjacent metal unit; t is_N,uIs T_c,uThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_S,uIs T_c,uCorrespondence sheetThe temperature of the adjacent metal unit on the side of the Yunan; psi_Nv,zIs the unit thermal resistance of the metal via; t is_c,dvThe temperature of the lower surface metal layer unit connected through the via hole is shown;

similarly, a discrete heat balance equation corresponding to the lower surface metal layer unit can be obtained:

wherein q is_c,dThe joule heating rate distribution per unit volume of the lower surface metal layer; t is_c,dIs the lower surface metal unit temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_d,dIs T_c,dThe longitudinal heat flux density of the lower surface of the corresponding unit can be the density of heat flux transmitted between the corresponding unit and the heating device; q. q.s_u,dIs T_c,uThe longitudinal heat flux density transmitted between the upper surface of the corresponding unit and the insulating layer; t is_E,dIs T_c,dThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_W,dIs T_c,dThe temperature of the corresponding unit west side adjacent metal unit; t is_N,dIs T_c,dThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_S,dIs T_c,dThe temperature of the adjacent metal unit on the south side of the corresponding unit; t is_c,uvThe temperature of the lower surface metal layer unit connected through the via hole is shown;

and (3) inducing a thermal diffusion matrix equation of the upper and lower surface metal layers based on numerical value dispersion:

wherein q is_MJ,uAnd q is_MJ,dRespectively corresponding unit joule heating vectors of the upper surface metal layer and the lower surface metal layer; h is_topThe heat transfer coefficient between the surface of the component and the outside is shown; m_MV,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; m_MV,dBy detachment of the metal layer unit from the lower surfaceThe coefficient corresponding to the temperature of the metal layer unit in the heat dissipation balance equation is formed; m_V,dThe coefficient of the temperature of the lower surface metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; m_V,uThe coefficient corresponding to the temperature of the upper surface metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; due to q_h,u、q_h,dThe compensated heat flow directly transmitted through the surface of the insulating layer is counted in an analytic solution matrix equation, and only the heat flow transmitted between the heating device and the part of the covering metal layer is counted in a lumped matrix equation of heat diffusion, so that C_M,u、C_M,dWill also include partial all-zero column vectors;

q is to be_h,uAnd q is_h,dSubstituting the expression of (a) into the thermal diffusion matrix equation to obtain:

wherein, C_Mu,uAnd C_Md,dRespectively for realizing T in the above equation_J,uAnd T_uRelated item and T_J,dAnd T_dThe correlation term is equivalent to the transformed coefficients.

Preferably, in step S4, the method includes constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the modified insulating layer temperature analytic solution matrix equation and the metal layer numerical value discrete-based thermal diffusion matrix equation, where the method includes:

the thermal analysis coupling equation set of the double-sided PCB structure is as follows:

wherein, T_u,T_d,q_M,u,q_M,dIf the vector is unknown, and other vector and matrix coefficients are known, the unknown temperature vector and the unknown heat flux density vector can be obtained by combining the four matrix equations.

Preferably, in step S5, the method approximates a discrete metal layer current continuity equation based on a finite volume method, calculates and summarizes a lumped matrix equation describing a linear relationship between unit potentials of the metal layer, and calculates a potential distribution in the metal layer, so as to further calculate a current density distribution and a joule heating distribution of the metal layer, and calculate the joule heating distribution into a metal layer thermal diffusion matrix equation, and perform iterative calculation with a temperature distribution to calculate a temperature distribution of a surface where each layer of the PCB is located under the condition of joule heating of the circuit, including:

based on second-order Taylor series approximation, obtaining an approximate discrete equation corresponding to a current continuity equation of a square area which is surrounded by the sections of Pn, Pw, Ps and Pe and takes the vertex Pc of the metal layer unit as the center:

wherein, V_Pc、V_PE、V_PW、V_PNAnd V_PSPc and P are respectively_E、P_W、P_NAnd the potential at PS; wherein σ_Pe、σ_Pw、σ_PnAnd σ_PsRespectively, respectively

Electrical conductivities at Pe, Pw, Pn and Ps;

the above is to apply the finite volume method to the control volume centered at the vertex;

the electrical conductivity of a metal varies with the temperature of the material, and its relationship to temperature can be represented by the following equation:

where ρ is the resistivity, ρ₀For resistivity at initial ambient temperature, Δ T_MTemperature rise of discrete metal units compared to ambient temperature, α_TAlpha for copper, for a material dependent temperature coefficient_TIs positive, i.e. its resistivity increases with increasing temperature;

the conductivity of the midpoint of the cross section is considered to be right and left or upper and lower sides thereofThe conductivity corresponding to the equivalent resistance after the metal unit synthesis is the reciprocal of the sum of the two unit resistivities, sigma_PsConductivity σ Using C1 cell and C2 cell_C1And σ_C2Represents:

therefore, the electrical conductivity of each unit of the metal layer is calculated according to the temperature distribution obtained in step S4; the conductivity corresponding to Pc can be approximately found by applying the conductivities corresponding to the four midpoints surrounding Pc:

for metal cells at line boundary locations, it may happen that the analyzed control volume cells are non-square, and the corresponding approximately discrete current continuity equation is:

wherein σ_Pm1、σ_Pm2、σ_Pm3And σ_Pm4Are each P_m1、P_m2、P_m3And P_m4The electrical conductivity of (d); v_Pvt1、V_Pvt2、V_Pvt3And V_Pvt4Potentials at Pvt1, Pvt2, Pvt2 and Pvt2, respectively;

in the PCB design process, 45 degree diagonal lines are designed, but they can be reconstructed by a square grid, and the discrete approximate current continuity equation of a small square control volume cell of quarter cell area with Ptr2 as the vertex at the line boundary is:

wherein σ_PmeAnd σ_PmnAre each P_meAnd P_mnThe electrical conductivity of (d); v_Ptr2、V_Pc2And V_Pc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation shows that the current continuity equation is the same as that corresponding to the triangular control unit surrounding Ptr2 after filling the sawtooth gaps of the boundary, i.e. it proves that the partial sawtooth boundaries formed discretely by the square grid and the boundary of the 45-degree inclined line can be approximately equivalent in current continuity;

while the control volume units at the line corners are not directly equivalent to the actual line corners, the control volume units with the cross sections of Pmn3, Pmw3 and Pms3 surrounding the boundary vertex Pb3 have the length d of the cross section of the lower side Pms3 in the sawtooth boundary_cA/2, while the length in the actual line is d_cTherefore, such corners need to be individually processed in the algorithm to account for the approximate current continuity equation corresponding to their actual boundaries:

wherein σ_Pmw3、σ_Pmn3And σ_Pms3Are each P_mw3、P_mn3And P_ms3The electrical conductivity of (d); v_Pb3、V_Pb4、V_Pc4And V_Pc5Potentials at Pb3, Pb4, Pc4 and Pc5, respectively;

for the metal layer units connected by the via holes, the corresponding approximately discrete current continuity equation can also be obtained according to the above vertex-centered finite volume method, and assuming that vertices of adjacent metal layer units, of which the metal layer unit vertex Pc is connected by the via holes, are Pu and Pd, the corresponding approximately discrete current continuity equation becomes:

wherein, V_PuAnd V_PdAre each P_uAnd P_dThe potential of (d); sigma_PuAnd σ_PdElectrical conductivity at Pu and Pd, respectively;

finally, a summary matrix equation of approximate discrete current continuity equations corresponding to the upper and lower surface metal layers can be obtained:

wherein, V₁And V₂Respectively corresponding to the summary vectors of all the discrete unit top potential with unknown potential of the upper surface metal layer and the lower surface metal layer; m_σ,1And M_σ,2Respectively as V in the upper surface matrix equation₁With V in the lower surface matrix equation₂The coefficient matrix of (a), is formed by conductivity parameters; k_p,1And K_p,2Known terms in two matrix equations, such as corresponding terms describing known endpoint potentials or cross-sectional currents, are included respectively; m_{σV,1_2}Is V in the upper surface matrix equation₂The coefficient matrix of (A) is composed of conductivity parameters corresponding to the lower surface metal layer unit connected by the via hole and the upper surface metal unit discrete heat balance equation, for V₂In cells not connected by vias, M_{σV,1_2}The corresponding column in (1) is an all-zero column; m_{σV,2_1}For V in the lower surface matrix equation₁The coefficient matrix of (a) is formed by conductivity parameters corresponding to the upper surface metal layer unit connected by the via hole in the discrete thermal equilibrium equation of the lower surface metal unit, for V₁In cells not connected by vias, M_{σV,2_1}The corresponding column in (1) is an all-zero column; the above matrix equation can therefore also incorporate the approximate discrete current continuity equation for metal cells not connected to a via; v can be obtained by simultaneously solving the two matrix equations₁And V₂；

After the potential distribution of the top point in the line is obtained, the potentials corresponding to the center of each square cell and the center of the cross section around the cell are obtained, and the potential V at the center of the cell C2 is obtained by applying the finite volume method to the small diamond region around the center of the cell C2_C2：

Wherein, V_Pvt4Potential at cell C2 apex Pvt 4; similarly, the potential V of the point Ps at the center of the cross section can be obtained_Ps：

Wherein, V_C1Potential at the center of cell C1;

thus, the lateral current density and the vertical current density corresponding to each cell center can be approximated, and the lateral current density J corresponding to the cell center can be approximated for the cell C1_x,C1And longitudinal current density J_y,C1Can be expressed as:

then C1 corresponds to a total current density of:

for the triangular unit at the oblique line boundary, assuming that the potential is linearly distributed in the region, the current density of the corresponding triangular region is obtained, and the linear potential distribution corresponding to the triangle formed by three points Ptr2, Ptr1 and Pc1 can be expressed as follows:

wherein, the third-order square matrix in the equation is composed of operation formulas corresponding to the coordinates of three vertexes under a Cartesian coordinate system, S_ΔIs a triangular area, so that the coefficient a corresponding to the potential linear distribution expression can be obtained₀、a₁And a₂Then, the transverse current density J within the triangle is known_x,ΔAnd longitudinal current density J_y,ΔComprises the following steps:

wherein σ_x,ΔConductivity corresponding to the triangular region;

on the premise of obtaining the current density corresponding to the metal cell, the joule heating rate per unit volume can be obtained, for example, for cell C1, the joule heating rate per unit volume q can be obtained_C1Can be expressed as:

q_C1＝E_C1·J_C1＝(J_x,C1 ²+J_y,C1 ²)/σ_C1＝J_C1 ²/σ_C1

wherein E is_C1The electric field strength at the center of cell C1;

the joule heating rate is related to the electrical conductivity, and the electrical conductivity is affected by the temperature, so when the joule heating of the line is counted, iterative operation between the temperature distribution and the joule heating rate distribution is required to be carried out, when the error of the iterative operation is smaller than a set error, the temperature distribution and the joule heating rate distribution can be approximately calculated, and the error can be set as the maximum value of the allowable unit temperature change in two iterative calculations. q. q.s_C1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the upper and lower surface metal layers_MJ,uAnd q is_MJ,d。

The steady-state thermal analysis method for the double-sided PCB structure provided by the invention has the following beneficial effects:

the method for coupling numerical value dispersion and temperature analysis solution can avoid integral dispersion and solution of the integral PCB structure by solving the temperature distribution temperature analysis solution of the designated areas on the upper surface and the lower surface of the substrate insulating layer, can realize approximate representation of the thermal diffusion between the metal circuit layers through the metal through holes and the metal layers on the basis of a finite volume method, can approximately correct the homogeneous structure hypothesis of the temperature analysis solution on the substrate insulating layer containing the through holes, can simultaneously count the boundary condition that the surface heat flow density distribution of the substrate insulating layer is subjected to the thermal diffusion effect of the metal layers, and further realizes the mutual coupling solution of the temperature analysis solution matrix equation and the metal layer thermal diffusion matrix equation to obtain the surface temperature distribution of the PCB.

Drawings

Fig. 1 is a schematic diagram of structured discrete metal layer unit heat conduction.

FIG. 2 is a cross-sectional view of a via.

Fig. 3 is a schematic diagram of a discrete unit of the circuit 1.

Fig. 4 is a schematic diagram of a discrete element of the circuit 2.

FIG. 5 is a schematic diagram of a zigzag boundary equivalent to a 45 degree diagonal boundary.

Fig. 6 is a schematic diagram of a corner case where a 45 degree boundary may occur.

FIG. 7 shows a general structure and flow of a corresponding procedure of a dual-sided PCB structure coupled thermal analysis modeling method.

Fig. 8 shows the layout and temperature distribution of the power MOSFET pad leads.

Fig. 9 is a circuit layout of upper and lower surface pad pins and upper surface temperature distribution after the addition of thermally conductive vias.

Fig. 10 is a diagram of a double-sided wiring structure.

Fig. 11 is a top surface line joule heating rate distribution.

Fig. 12 is a lower surface line joule heating rate distribution.

Fig. 13 is a top surface line temperature distribution.

Fig. 14 is a bottom surface line temperature distribution.

Detailed Description

The following description of the embodiments of the present invention is provided to facilitate the understanding of the present invention by those skilled in the art, but it should be understood that the present invention is not limited to the scope of the embodiments, and it will be apparent to those skilled in the art that various changes may be made without departing from the spirit and scope of the invention as defined and defined in the appended claims, and all matters produced by the invention using the inventive concept are protected.

According to one embodiment of the application, the steady-state thermal analysis method for the double-sided PCB structure comprises the following steps:

The above steps will be described in detail below according to one embodiment of the present application.

Firstly, thermal analysis premise hypothesis is required to be carried out on the structural characteristics of the double-sided PCB:

the double-sided PCB structure is mainly composed of a core substrate insulating layer and metal circuit layers covered on the upper and lower surfaces of the core substrate insulating layer. A common substrate layer material is FR-4 glass epoxy, which has a thermal conductivity of about 0.3W/mK and a crystallization temperature of about 120 ℃, and when this temperature is exceeded, both the mechanical and insulating properties of FR-4 will be significantly reduced. The metal layer is mainly made of copper, and the thermal conductivity is about 390W/mK. While the thickness of the substrate insulating layer (>0.5mm) is much larger than the thickness of the metal layer (generally less than 0.1 mm). The PCB surface is covered with a soldemask solder resist coating and a silk screen printing layer for marking devices and printing trademarks and other information. Because the area of the silk-screen layer on the surface of the PCB is low, the influence of the silk-screen layer on the heat transfer is ignored. In addition, because the solder resist coating is thin (typically 1mil) and has low thermal conductivity (0.25W/mK), its equivalent thermal resistance is taken into account in the heat transfer coefficient between the PCB surface and the external environment. The Heat Transfer Coefficient (HTC) is also one of the important thermal boundary conditions for deriving a temperature-resolved solution.

However, the existence of the components will destroy the assumed condition that the upper and lower surfaces of the PCB and the external environment exchange heat under a certain HTC, and the PCB and the external environment also exchange heat through the surfaces of the components, so that heat flow transmission also exists between the components and the PCB, and the package of the components will also affect the assumed thermal boundary condition.

Although the insulating material occupies most of the PCB structure, the large difference in thermal conductivity between the metal layer and the insulating layer results in the metal layer wiring having much higher thermal diffusivity than the insulating material even in the micrometer-sized thickness, which cannot be ignored in the thermal analysis. However, since the metal layer is thin and has high thermal conductivity, the temperature gradient between the upper and lower surfaces thereof can be ignored, and the insulating layer is closely attached to the metal layer, so that it can be assumed that the temperature distribution of the metal layer is the same as that of the surface region of the insulating layer in contact with the metal layer. The substrate insulating layer is a main constituent material of the PCB, but since it is not a heat source, the temperature distribution inside the substrate insulating layer is not a main solution target.

Therefore, in the modeling method, the temperature distribution analytic solution of the surface of the insulating layer is used for providing the temperature distribution expression of the metal layer and the heating device which are in contact with the insulating layer, namely the temperature distribution expression is used as one of the coupling variables of the numerical discrete analysis and the analytic solution; meanwhile, a numerical value dispersion method based on finite volume is used for counting the heat diffusion inside the metal layer and the heat conduction process between the metal layers through the metal through holes, and an approximate steady-state heat diffusion matrix equation of the metal layer is given, wherein the heat diffusion matrix equation contains the heat flow density distribution actually transmitted by the metal layer to the insulating layer, and the distribution is also used as one of coupling variables; and finally solving the coupling variable through a simultaneous analytic solution matrix equation and a thermal diffusion matrix equation. And finally, obtaining the core temperature of the heating element by using the thermal resistance parameters of the heating element.

In addition, because the heat dissipation areas of the four sides of the PCB are generally much smaller than the heat dissipation areas of the upper and lower surfaces of the PCB, and the main heat sources are both on the upper and lower surfaces, the heat exchange between the four sides of the PCB and the outside is neglected. The method is therefore not suitable for thermal analysis of PCBs of the type in which heat transfer is achieved primarily through the side and where the side has a plug section, etc. in an application.

The following detailed description of the algorithm is made under the above-described assumption:

step S1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer, and obtaining an insulating layer temperature analytic solution matrix equation, wherein the equation specifically comprises the following steps:

the corresponding steady state heat balance Laplace equation of the insulating layer in the double-sided PCB and the boundary conditions of the periphery, the upper surface and the lower surface are as follows:

wherein T is the temperature variation compared with the ambient temperature or the temperature when the ambient temperature is 0 ℃, and the density distribution function of the heat transmitted from the upper surface and the lower surface is q_iu(x, y) and q_id(x,y)，h_u(W/m²K) And h_d(W/m²K)Heat Transfer Coefficients (HTC), k, of the upper and lower surfaces of the insulating layer, respectively_i(W/mK) is insulation layer thermal conductivity, L_x、L_y、 L_inThe length of the insulating layer in the x and y directions and the thickness of the insulating layer in the z direction in a cartesian coordinate system, respectively.

according to the principle of thermal effect superposition, the thermal effect generated by the upper surface device and the lower surface device or other heat sources in the insulating layer can be regarded as the result of superposition of the thermal effects after the upper surface heat source and the lower surface heat source act independently, and at the moment, the equation (1) is changed into:

therefore, only the analytic solution corresponding to the surface temperature distribution of the insulating layer with the single-sided structure needs to be deduced, and the analytic solution corresponding to the double-sided structure can be obtained through superposition. Only the solving process of the thermal effect corresponding to the analytic solution theta generated by the upper surface heat source will be described, and the solving process of eta is similar to the solving process. Substituting the order solution 2 into the lower surface (z ═ L)_in) Corresponding thermal boundary conditions can be obtained, when h is reached, the corresponding relation between partial coefficient and characteristic value in solution can be obtained_dWhen not zero, C₂And C_γn,mThe corresponding expression is:

when h is_dIs zero, C₁,C_γn,mCan be expressed as:

discussion of h only_dGeneral derivation of other coefficients in an example analytic solution corresponding to a single-sided structure when non-zero, h_dThe derivation process for zero time is similar and will not be described here. The number of stages is substituted into the thermal boundary condition of the upper surface (z is 0), and the following results are obtained:

multiplying both sides of the above equation by cos (. beta.)_nx)cos(μ_my) and integration is performed in the surface heating area, the integration after multiplication of different numbers of stages is zero according to the characteristics of trigonometric function, so that the left side of the equation only contains cos²(β_nx)cos²(μ_my) is not zero after integration, and the right side of the equation is used for each heating unit q in the heating area_iu,iIntegrating and then summing, and calculating to obtain coefficient C₁And C_n,m：

Wherein d is_qIntegration Domain Sq for grid cell Length after discretization of surface Heat Source grid_iu,iI.e. representing an incoming heat flow density of q_iu,i(W/m²) Surface discrete unit area of (d)_nAnd delta_mTo resolve dimensionless intermediate parameters of the derivation process, q_iu,1～q_iu,N1In order to discretize the unit heat flux density of the heat source grids on the surface of the insulating layer, if the heat flux density transmitted by a certain heat source is uniform, grid discretization is not carried out on the longitudinal heat flux density transmitted by the area where the heat source is located; c₁The molecule in the expression is the heat quantity transmitted in the whole unit time, C_n,mThe double integral summation term in the expression can be directlyIntegration is performed in the area where the heat source is located.

Further, the summation operation of Fourier series can be expressed by means of vector operation, and for the above analytic solution, q_iu,1～q_iu,N1Namely, a column vector q of heat flux density is constructed_iuThen, the analytic solution of any point in the single-sided insulating layer can be expressed as:

solving the temperature distribution of the designated area is the thermal resistance row vector R_u(x,y,z)Expanding into a thermal resistance matrix R according to the coordinates of the corresponding distribution points of the region_uWherein each row is corresponding to a thermal resistance row vector for solving the temperature at a certain point, and then R is added_uMultiplying the temperature distribution vector by the corresponding heat flow density vector q to obtain a temperature distribution vector of the specified area, namely the heat effect generated by the corresponding q in the specified area; and finally, superposing Fourier series analytic solutions in a matrix form corresponding to theta and eta of the specified region to obtain the temperature distribution of the region.

wherein, the distribution vectors of the heat flux density transmitted into the insulating layer by the upper and lower surface metal layers are q_M,u、q_M,dWherein q is_uAnd q is_dThe distribution vector q of heat flux density transferred between the upper and lower surface heating components and the covering region thereof_uComprising q_iu(x, y) heat flux density distribution after dispersion, q_dComprising q_id(x, y) a heat flow density distribution after discretization; r_Mu,uIs q_M,uAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Mu,dIs q_M,dAt T_uThe defined area generating a thermal effectA corresponding thermal resistivity matrix; r_Md,uIs q_M,uAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Md,dq_M,dAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_iu,uIs q_uAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_iu,dIs q_dAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_id,uIs q_uAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_id,dIs q_dAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; due to q_uAnd q is_dThe heat flow transmitted through the surface of the metal layer is counted in the heat diffusion matrix equation corresponding to the metal layer, and only q is counted in the analytic solution matrix equation_iu(x, y) and q_id(x, y) heat flow transferred between the associated top and bottom surface heat generating components and the surface area of the insulating layer directly covered thereby, R_iu,u、R_iu,d、R_id,uAnd R_id,dWill include a partial all-zero column vector;

step S2, calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the thermal boundary condition of the substrate insulating layer, and obtaining a corrected insulating layer temperature analytic solution matrix equation expression, which specifically comprises the following steps:

the thermal boundary conditions according to which the analytical solution of the surface temperature distribution of the insulating layer of the substrate is solved are mainly the heat flux density distribution transmitted by the upper surface and the lower surface of the insulating layer of the substrate and the heat transfer coefficient for representing the heat exchange between the surface and the external environment. However, the thermal boundary conditions ignore the fact that the partial region covered by the heat generating component is only in heat exchange with the component, rather than directly with the external environment, and do not account for heat exchange between the surface of the component and the external environment. To correct the assumed thermal boundary conditions, the area occupied by the components must be based on h_uThe multi-calculated heat exchange heat flow density is compensated, and on the other hand, the thermal resistance parameters provided by the components can be used for deducing the position between the related components and the PCB surfaceTo heat transfer.

Setting the heat flux density vector of the covering area of the upper surface heating device as q_h,uIt can be expressed as follows:

q_h,u＝q_u+h_uT_J,u (10)

wherein, T_J,uA temperature distribution vector for an area occupied by the upper surface device; because the covered area of the component does not exist in the range h between the component and the external environment_uThe calculated heat exchange under the condition, but the part of the heat exchange is taken into account in the analytical solution derivation process of the aforementioned step S1, so the occupied area of the component is determined according to h_uCompensating for the multi-calculated heat exchange; in addition, the thermal resistance parameters provided by the components can be used for deducing the heat transfer between the related components and the PCB surface, and the specific deduction process is as follows:

the two formulas are combined to obtain:

q_h,u＝C_P,uq_J,u-C_T,uT_J,u

wherein the temperature distribution of the upper surface element device surface is T_topThe heat transfer coefficient of the surface to the outside is h_topThe equivalent heat conduction corresponding to the heat sink existing on the surface thereof is also taken into account_top(ii) a Mean temperature of the core is T_j(ii) a The average thermal resistance between the surface and the core is psi_jt(ii) a The average thermal resistance between the junction center and the bottom of the device is psi_jc(ii) a Equivalent thermal resistance of solder mask coating is psi_c,u(ii) a The heat flux density transferred from the core to the surface of the component is q_jt(ii) a The heating volume density vector corresponding to the upper surface element device is P_J,u(ii) a The diagonal matrix composed of height parameters of the components is C_J,u(ii) a The density vector of the heating surface corresponding to the upper surface element device is q_J,u；C_P,uFor taking into account the thermal resistance parameter of the device_J,uThe coefficient of (a); c_T,uFor taking into account the corresponding T after the thermal resistance parameter of the device_J,uThe coefficient of (a);

Wherein, T_J,dThe temperature distribution vector of the coverage area of the lower surface heating device is obtained; r_iu,uC_P,uq_J,uThe term is a vector of known constants, and T_uIncluding T_J,uAll elements of the vector can be transformed by adding zero to the matrix coefficients in a rising-dimension manner to realize the use of T in the matrix equation_uInstead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_Ju,uT_u＝R_iu,uC_T,uT_J,u (14)

T_u＝R_Mu,uq_M,u+R_Md,uq_M,d+R_id,uq_d+R_iu,uC_P,uq_J,u-C_Ju,uT_u (15)

in the same way, the T corrected by the thermal boundary condition can be obtained_dCorresponding analytic solution matrix equation, and after correction, T_uAnd T_dThe corresponding simultaneous analytical solution matrix equation becomes:

wherein, with C_Ju,uAnd C_Jd,uFunction of (C) is similar to the definition_Ju,d、C_Jd,dRespectively for realizing T in the lower surface analytical solution matrix equation_J,uAnd T_uRelated item and T_J,dAnd T_dCoefficients of the correlation term equivalent transform;

s3, approximating a discrete metal layer heat balance equation based on a finite volume method, summarizing to obtain a metal layer heat diffusion matrix equation, and calculating the heat conduction between the metal layer and a surface heating device thereof and the heat conduction of the metal via hole therein, wherein the method specifically comprises the following steps:

s3.1, approximating a heat balance equation of a discrete metal layer based on a finite volume method;

the steady state thermal equilibrium equation for the metal layer has the form of poisson's equation:

wherein k is_m(W/mK) represents the metal layer thermal conductivity; q. q.s_c(x,y)(W/m³) Is a joule heating rate distribution function per unit volume of the metal layer.

To calculate the thermal diffusion in the metal layer, the above thermal equilibrium equation will be discretized mainly using a finite volume method, and the metal layer will be discretized using a structured square grid, as shown in fig. 1. Wherein the contribution of the metal layer to the heat diffusion is realized by the heat conduction among the metal layer units, and the heat flow density conducted between the central unit and the four surrounding units is q as shown in fig. 1_n,q_s,q_wAnd q is_e(W/m²) And the heat exchange between the top of the device and the component is expressed as heat flow density q_u(W/m²) The heat conduction with the bottom insulating layer is formed by q_d(i.e., constituting q in the aforementioned analytical solution equation)_M,uElements of a vector) (W/m)²) Represents, but onlyJoule heat in the interior of the element is generated by the volume heat generation rate q_c(W/m³) And (4) representing. When Joule heating of certain lines or elements is not considered in the thermal analysis, the corresponding q_cI.e. zero.

Performing triple volume integration on the above equation in a unit, and applying the divergence theorem, we can obtain:

wherein q is_cIs the unit volumetric heat rate, d_cIs the unit side length of the structured and discrete metal unit, d_mMultiplying both ends of the equation by k for the thickness of the metal layer_mAnd decomposing the unit closed surface integral into surface integrals on six surfaces thereof, so as to obtain:

and then based on the approximation of second-order Taylor series, the partial derivative term in the above formula is linearized, such as the temperature T of the central unit_cAnd its east cell temperature T_EThe corresponding Taylor-series expansion is:

wherein, the position (x)_e,y_e) The midpoint of the connecting line of the corresponding central unit and the right unit; o (d)_c ³) Is the sum of all terms containing partial derivatives of the third order and above, and d contained therein_cThe power exponent of the power term is greater than or equal to 3, since d is generally_cSmaller value of (a) can be represented by O (d)_c ³) Are omitted. And then subtracting the two expressions to approximately obtain the partial derivative of the temperature on the corresponding section:

similarly, other partial derivatives in equation 20 can be processed approximately linearly, and finally, the heat dissipation balance equation can be realized, and an approximate linear relationship between the temperatures of adjacent units is obtained:

therefore, the temperature relationship between adjacent units is established by the equation after dispersion, so that the subsequent dispersion heat balance equation corresponding to all the metal layer units can be conveniently induced into a lumped matrix equation representing the metal layer heat diffusion, and when the units are connected with other adjacent metal layer units through metal via holes, the q-axis heat resistance can be based on the metal via hole heat resistance_uOr q_dAnd (5) approximate linearization processing. The discrete heat balance equations corresponding to all the metal layer units are collected to form a lumped matrix equation representing the heat diffusion of the metal layers and the heat conduction between the layers through the via holes, so that the coupling solution of the equation with the analytic solution matrix is realized conveniently.

S3.2, calculating the heat conduction of the metal via hole in a metal layer separation heat dissipation balance equation;

there are typically metal vias between metal layers in a PCB structure, and fig. 2 shows a cross-sectional view of the vias, where the longitudinal thermal resistance can be expressed as follows:

wherein the thickness of the insulating layer is L_in(m) via radius r_v(m) via thermal conductivity k_v(W/mK)，d_pw(m) is the thickness of the inner wall of the via hole, the unit of thermal resistance is K/W, as the via hole material is generally copper, the thermal conductivity is higher, the via hole is smaller, the thermal conduction effect is still not negligible, if the via hole is not designed, the occupied space is completely filled into the FR-4 resin material of the PCB insulating layer, and the longitudinal thermal resistance corresponding to the part is calculated as:

wherein k is_in(W/mK) is the thermal conductivity of the FR-4 material, which is about 0.3. Due to the large difference of the thermal conductivity (k) between copper and FR-4 material in the normal working temperature range_cu: k_FR-4≈390：0.3>1000) And the ratio of the thermal resistance of the metal via to the thermal resistance of FR-4 is as follows:

it can be seen that as long as r_v<2000d_pw(d_pwTypically 1-4 mil), the thermal resistance of a thin-walled metal via will be lower than the thermal resistance of the same via volume filled with FR-4, and this condition is met in most PCB via designs; on the other hand, when the via hole is a solid via hole filled with copper completely, the via hole is generally used as a heat conducting via hole due to low thermal resistance; therefore, the thermal conduction contribution of the vias must be accounted for in the thermal analysis of the PCB structure.

Therefore, on the premise of taking the heat conduction of the metal via hole into account, the heat balance equation of the discrete metal layer on the upper surface can be obtained:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_c,uThe joule heating rate distribution per unit volume of the upper surface metal layer; t is_c,uIs the upper surface metal cell temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_u,uIs T_c,uThe longitudinal heat flow density of the upper surface of the corresponding unit can be the density of heat flow transferred between the corresponding unit and the heating device; q. q.s_d,uIs T_c,uThe longitudinal heat flux density transmitted between the lower surface of the corresponding unit and the insulating layer; t is_E,uIs T_c,uThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_W,uIs T_c,uThe temperature of the corresponding unit west side adjacent metal unit; t is_N,uIs T_c,uThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_S,uIs T_c,uThe temperature of adjacent metal units on the south side of the corresponding unit; psi_Nv,zIs the unit thermal resistance of the metal via; t is_c,dvThe temperature of the lower surface metal layer unit connected through the via hole is shown;

similarly, a discrete heat balance equation corresponding to the lower surface metal layer unit can also be obtained:

wherein q is_c,dThe joule heating rate distribution per unit volume of the lower surface metal layer; t is_c,dIs the lower surface metal unit temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_d,dIs T_c,dThe longitudinal heat flux density of the lower surface of the corresponding unit can be the density of heat flux transferred between the corresponding unit and the heating device; q. q.s_u,dIs T_c,uThe longitudinal heat flux density transmitted between the upper surface of the corresponding unit and the insulating layer; t is_E,dIs T_c,dThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_W,dIs T_c,dCorresponding to the temperature of the adjacent metal unit on the west side of the unit; t is_N,dIs T_c,dThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_S,dIs T_c,dThe temperature of the adjacent metal unit on the south side of the corresponding unit; t is_c,uvThe temperature of the lower surface metal layer unit connected through the via hole is shown;

step S3.3, the upper and lower surface metal layer thermal diffusion matrix equation based on numerical value dispersion is summarized:

wherein q is_MJ,uAnd q is_MJ,dRespectively corresponding unit joule heating vectors of the upper surface metal layer and the lower surface metal layer; h is_topHeat of the surface of the component and the outsideA transfer coefficient; m_MV,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; m_MV,dThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; m_V,dThe coefficient of the discrete heat balance equation of the upper surface metal layer unit corresponding to the temperature of the lower surface metal layer unit is formed; m_V,uThe coefficient corresponding to the temperature of the upper surface metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; due to q_h,u、 q_h,dThe compensated heat flow directly transmitted through the surface of the insulating layer is counted in an analytic solution matrix equation, and only the heat flow transmitted by the heating device and the part of the covering metal layer of the heating device is counted in a lumped matrix equation of heat diffusion, so that C_M,u、C_M,dWill also include partial all-zero column vectors;

wherein, C_Mu,uAnd C_Md,dRespectively for realizing T in the above equation_J,uAnd T_uRelated item and T_J,dAnd T_dCoefficients of the correlation term equivalent transform;

step S4, according to the corrected insulation layer temperature analytic solution matrix equation and the metal layer thermal diffusion matrix equation based on numerical value dispersion, constructing a thermal analysis coupling equation set of the double-sided PCB structure, which specifically comprises the following steps:

from the temperature-resolved solution matrix equation 17, and the thermal-analysis coupling equation set corresponding to the numerical-dispersion-based thermal-diffusion matrix equation 30, the double-sided PCB structure can be expressed as:

wherein, T_u,T_d,q_M,u,q_M,dIf the temperature vector is unknown, and other vector and matrix coefficients are known, the unknown temperature vector and the unknown heat flow density vector can be obtained by combining the four matrix equations.

S5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and obtaining potential distribution in the metal layer, so that current density distribution and Joule heating distribution of the metal layer can be further obtained, the Joule heating distribution can be included in a metal layer thermal diffusion matrix equation, iterative calculation is carried out between the Joule heating distribution and temperature distribution, and temperature distribution of the surface where each layer of the PCB is included under the premise of Joule heating of a circuit is obtained, and the method specifically comprises the following steps:

s5.1, approximating a discrete metal layer current continuity equation based on a vertex-centered finite volume method, which specifically comprises:

the current continuity equation in the metal layer can be expressed as:

wherein, J (A/m)²) For current density, σ (S/m) is the temperature-dependent conductivity, E (V/m) is the electric field strength, and V (V) is the electric potential; further applying the Gaussian divergence theorem, the volume and closed surface integral form of the equation can be obtained:

for a square region surrounded by the cross sections of Pn, Pw, Ps and Pe with the metal layer unit vertex Pc connected by the non-metal via hole as the center in fig. 3, the integral corresponding to the closed surface can be represented by the sum of the integrals corresponding to the cross sections around the region:

wherein σ_Pe、σ_Pw、σ_PnAnd σ_PsElectrical conductivities at Pe, Pw, Pn and Ps, respectively; j. the design is a square_Pe、J_PwThe transverse current densities at Pe and Pw, respectively; j. the design is a square_PnAnd J_PsThe longitudinal current densities at Pn and Ps, respectively;

the above equation also describes that the total current flowing into the region is equal to the total current flowing out;

to vertex Pc and P_EThe potential of (c) is developed using a taylor series:

wherein, V_PeIs P_eThe potential of (d); (x)_Pe，y_Pe) Is P_eCoordinates of (c); d_cIs the unit length of the metal layer unit; further subtracting the two to obtain a partial differential expression based on the second-order taylor series approximation:

therefore, an approximate discrete equation corresponding to a current continuity equation of a square region surrounded by the cross section where Pn, Pw, Ps and Pe are located can be obtained by using the metal layer unit vertex Pc shown in fig. 3 as the center (not connected with the metal via):

wherein, V_Pc、V_PE、V_PW、V_PNAnd V_PSPc and P are respectively_E、P_W、P_NAnd the potential at PS;

the conductivity of the midpoint of the cross section is considered to be the conductivity corresponding to the equivalent resistance of the synthesized metal units on the left and right sides or the upper and lower sides, i.e. the reciprocal of the sum of the resistivities of the two units, e.g. sigma_PsThe conductivities of the C1 cell and the C2 cell can be used_C1And σ_C2Represents:

the electrical conductivity of each cell of the metal layer is therefore calculated from the temperature distribution obtained in claim 5; the conductivity corresponding to Pc can be approximately found by applying the conductivities corresponding to the four midpoints surrounding Pc:

for metal cells at line boundary positions, it may happen that the control volume cell to be analyzed is non-square, such as the control cell No. 2 (which is the smaller square control cell No. 1 below) surrounded by dark black multi-point lines shown in fig. 4, and its corresponding approximately discrete current continuity equation is:

in the PCB design process, 45 degree inclined lines are often designed, but the square grid can be better reconstructed, for example, for a small square control volume unit with a quarter cell area with Ptr2 as the vertex at the boundary of the upper left corner line of fig. 5, the discrete approximate current continuity equation is:

wherein σ_PmeAnd σ_PmnAre each P_meAnd P_mnThe electrical conductivity of (d); v_Ptr2、V_Pc2And V_Pc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation indicates that the current continuity equation is the same as the current continuity equation corresponding to the triangle control cell surrounding Ptr2 after filling the boundary sawtooth gap. That is, it is proved that the partial zigzag boundary formed discretely in the square grid and the boundary of the 45-degree inclined line are approximately equivalent in current continuity.

The control volume units at the line corners are not directly equivalent to the actual line corners, for example, the control volume units with the cross sections of Pmn3, Pmw3 and Pms3 around the boundary vertex Pb3 have the length d of the cross section of the lower side Pms3 in the sawtooth boundary_cA/2, while the length in the actual line is d_cTherefore, such corners need to be individually processed in the thermal analysis method to account for the approximate current continuity equation corresponding to their actual boundaries:

wherein σ_Pmw3、σ_Pmn3And σ_Pms3Are each P_mw3、P_mn3And P_ms3The electrical conductivity of (d); v_Pb3、V_Pb4、V_Pc4And V_Pc5Potentials at Pb3, Pb4, Pc4 and Pc5, respectively; fig. 6 shows a possible corner situation for a 45 degree diagonal line and its corresponding control volume range.

For the metal layer units connected by the vias, the corresponding approximately discrete current continuity equation can also be obtained according to the above vertex-centered finite volume method, and assuming that the vertices of the adjacent metal layer units connected by the metal layer unit vertex Pc through the vias shown in fig. 3 are Pu and Pd, the corresponding approximately discrete current continuity equation becomes:

s5.2, summarizing to obtain a lumped matrix equation for describing the linear relation among the electric potentials of the metal layer units, and solving the electric potential distribution of the top points of the metal layer units, wherein the lumped matrix equation specifically comprises the following steps:

summarizing and summarizing the discrete current continuity equations corresponding to the metal layer units of various types to obtain a matrix equation of an approximate discrete current continuity equation corresponding to each metal layer, and finally obtaining a summarized matrix equation of the approximate discrete current continuity equations corresponding to the upper and lower surface metal layers:

wherein, V₁And V₂Respectively corresponding to the summary vectors of all the discrete unit top potential with unknown potential of the upper surface metal layer and the lower surface metal layer; m_σ,1And M_σ,2Respectively V in the upper surface matrix equation₁With the lower surface matrix equationV₂The coefficient matrix of (a), is formed by conductivity parameters; k_p,1And K_p,2Known terms in two matrix equations, such as corresponding terms describing known endpoint potentials or cross-sectional currents, are included respectively; m_{σV,1_2}Is V in the upper surface matrix equation₂The coefficient matrix of (a) is formed by conductivity parameters corresponding to the lower surface metal layer unit connected by the via hole and the upper surface metal unit discrete heat balance equation, for V₂In cells not connected by vias, M_{σV,1_2}The corresponding column in (1) is an all-zero column; m_{σV,2_1}For V in the lower surface matrix equation₁The coefficient matrix of (a) is formed by the conductivity parameters corresponding to the upper surface metal layer unit connected by the via hole in the discrete thermal equilibrium equation of the lower surface metal unit, for V₁In cells not connected by vias, M_{σV,2_1}The corresponding column in (1) is an all-zero column; the above matrix equation can therefore also incorporate the approximate discrete current continuity equation for metal cells not connected to a via; v can be obtained by simultaneously solving the two matrix equations₁And V₂；

S5.3, solving the potential distribution of the center of the unit in the metal layer and the midpoint of the cross section around the unit, wherein the method specifically comprises the following steps:

after the electrical potential distribution at the vertices of the line is determined, the corresponding electrical potentials at the center of each square cell and at the center of the cross-section around the cell can be determined, for example, by applying the finite volume method to the small diamond region around the center of cell C2 in FIG. 3 to determine the electrical potential V at the center of cell C2_C2：

Wherein, V_C1Potential at the center of cell C1;

step S5.4, further solving the current density distribution and the Joule heating distribution of the metal layer, and calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with the temperature distribution to obtain the temperature distribution of the surface where each layer of the PCB is located under the premise of calculating the Joule heating of the circuit, which specifically comprises the following steps:

on the basis of the potential distribution of the metal layer, the transverse current density and the longitudinal current density corresponding to the center of each cell can be approximately obtained, for example, for the cell C1, the transverse current density J corresponding to the center of the cell_x,C1And longitudinal current density J_y,C1Can be expressed as:

then C1 corresponds to a total current density of:

for the triangular cells at the oblique line boundary, the electric potential is assumed to be linearly distributed in the region, so as to obtain the current density of the corresponding triangular region, for example, the linear electric potential distribution corresponding to the triangle enclosed by three points Ptr2, Ptr1 and Pc1 can be expressed as follows:

wherein σ_x,ΔConductivity corresponding to the triangular region;

q_C1＝E_C1·J_C1＝(J_x,C1 ²+J_y,C1 ²)/σ_C1＝J_C1 ²/σ_C1 (52)

wherein E is_C1The electric field strength at the center of cell C1; q. q.s_C1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the upper and lower surface metal layers_MJ,uAnd q is_MJ,d。

The joule heating rate is related to the electrical conductivity, and the electrical conductivity is affected by the temperature, so when the joule heating of the line is counted, the iterative operation between the temperature distribution and the joule heating rate distribution needs to be performed, when the error of the iterative operation is smaller than the set error, the temperature distribution and the joule heating rate distribution can be approximately calculated, and the error can be set as the maximum value of the unit temperature change allowed in the two iterative calculations, as shown in fig. 7.

The general program structure for realizing the thermal analysis modeling of the double-sided PCB structure is as follows:

the general structure and flow of the PCB structure thermal analysis program which can be written by using the method are shown in FIG. 7, and the program can be divided into an initialization program and a main program. If an analysis model is newly established in an initialization program, a circuit diagram is firstly imported, all components such as a circuit, a device and a bonding pad need to be identified and position confirmed in the process of analyzing the circuit by the program, and then a coefficient matrix required by summarizing a coupling equation set can be calculated, wherein the coefficient matrix comprises coefficient matrixes required by analyzing an antipyretic resistance coefficient matrix and a numerical value discrete equation and taken into account of heat conduction of a metal layer and thermal resistance parameters of the device and a radiator. If the existing model is modified, the step of setting/modifying system parameters can be directly carried out, and finally, the data is stored.

The main program comprises two parts of direct coupling calculation temperature distribution and iterative calculation temperature distribution, and the direct coupling calculation calculates a thermal analysis result according to a coupling equation set by calling each coefficient matrix and a heating density vector of a known device. If the joule heating of part of the metal layer circuit needs to be taken into account, the influence of the temperature on the metal conductivity needs to be taken into account, so that iterative operation between the joule heating and the temperature distribution is needed to obtain the steady-state circuit temperature distribution. In the iterative calculation part, the Joule heating density distribution of the circuit is calculated according to the assumed temperature distribution of the circuit, then the coupling solution of the corresponding temperature distribution is calculated, the new Joule heating density distribution of the circuit is calculated according to the temperature distribution, then the new coupling solution is calculated, and the iterative calculation is carried out in sequence until the maximum difference value of the unit temperatures in the two temperature distributions is lower than the given error limit value.

According to one embodiment of the present application, a specific PCB structure is analyzed to verify the computational accuracy of the present algorithm:

and the MATLAB is utilized to complete the double-sided PCB structure thermal analysis modeling program capable of analyzing the two-dimensional line structure diagram based on the method. In the program, the thermal parameters such as the size of the PCB, the thickness of the metal layer and the insulating layer, the ambient temperature, the thermal conductivity and the surface heat transfer coefficient of the material, the thermal resistance of surface coatings such as a solder mask, the thermal resistance of the device and the like can be set, and the functions of setting the electrical parameters such as the power consumption of the device, the current and the voltage of the double-end line, the branch current of the multi-branch line, the terminal potential of. The previous studies of the inventor have concluded that a sufficient requirement for the known electrical parameters of the line, which are required for the potential distribution and current density analysis of the multi-branch line, is that the sum of the known branch currents or the number of branch terminal potentials needs to be equal to the number of line terminals, and at least one terminal potential is known.

a. Power MOSFET thermal analysis modeling example

An example of thermal analysis of a local MOS pin pad structure on a PCB and further optimization of the structure is first given. A power MOS with the model number of STB20N95K5 is used as a main circuit switch on a double-sided PCB of a certain vehicle-mounted DC-DC converter, and specific parameters of a pin pad structure of the MOS are shown in a table below.

TABLE 1 MOSFET pad Pin line parameters

Size of substrate

Thickness of copper foil

Via hole wall thickness

Pore diameter

Heat generation rate

Bottom cold plate temperature

Type of package

22x25x2(mm)

35μm

25.4μm(1mil)

1mm

3.4W

60℃

D²PAK

The original structure and the thermal analysis temperature distribution result of the MOS pin pad are shown in fig. 8, which is the most important heating device in the converter PCB, and here, the contribution of other parts of the PCB to the MOS heat conduction is ignored in the analysis process, and only the heat conduction from the PCB to the bottom cold plate is taken into account, so as to simulate the theoretical temperature rise extreme value of the MOS under the working condition. And simultaneously modeling the structure in COMSOL and carrying out thermal analysis so as to check the reasonability and the accuracy of the modeling method through comparison. As can be seen from the analysis results, the thermal distribution of the analysis results obtained based on the thermal analysis method is basically consistent with the COMSOL analysis results, and the temperature difference of the two hot spots accounts for less than 1% of the total hot spot temperature rise. However, the hot spot temperature of the lead bonding pad region exceeds 195 ℃, and the MOS cannot work normally, so that the lead bonding pad structure needs to be improved.

Because the bottom is a constant temperature cold plate, in order to enhance the heat conduction from the MOS to the bottom, some metal through holes are added in the pin pad area, and a layer of copper-clad metal layer is also added at the bottom, the structure and the thermal analysis temperature distribution result are shown in fig. 9, and the related analysis parameters and the structure discrete parameters corresponding to the coupled thermal analysis program and the COMSOL model are given at the same time.

TABLE 2 MOSFET pad pin line simulation calculation parameters

Therefore, under the condition of adding the metal heat conduction through hole, the temperature rise of the hot spot of the welding pad area is reduced by over 90 ℃, the core temperature of the welding pad area does not exceed 106 ℃ according to the calculation of the related thermal resistance parameters of the MOS package, and the welding pad area is in a safe working range. Compared with simulation results of COMSOL, the simulation results show that on one hand, the operation principle of the coupling thermal analysis is reasonable and accurate, and on the other hand, on the premise of having similar precision, the total number of numerical discrete units of the PCB in the coupling thermal analysis related by the invention is only 7495 of the pad area, which is far lower than the number of units of the PCB which is wholly discrete by COMSOL based on a finite element analysis method, so that the method has obvious advantages in numerical analysis efficiency.

b. Modeling example for joule heating thermal analysis of double-sided PCB (printed Circuit Board) circuit

In the following, the joule heating and the temperature distribution of a double-sided line having a via hole are simulated, calculated and compared, and the analysis of joule heating will involve the iterative operation mentioned in the foregoing program structure. The basic initial parameters of the circuit are listed in tables 3 and 4, wherein the negative current represents the current flowing into the corresponding end point, the circuit structure diagrams and the corresponding simulation calculation results are shown in fig. 10 to 14, fig. 11 and 12 are the joule heating rates per unit volume of the upper and lower surfaces, respectively, and fig. 13 and 14 are the temperature distributions of the upper and lower surfaces, respectively.

TABLE 3 two-sided line basic parameters

Size of substrate

Thickness of copper foil

Via hole wall thickness

Pore diameter

Pad1

Pad2

Pad3

Via1(2nd layer)

Via2

70x34x1.6(mm)

18μm

25.4μm(1mil)

2mm

-3A

-5A

-3A

12V

TABLE 4 double-sided line simulation calculation parameters

Therefore, the coupling thermal analysis operation result based on the method keeps higher consistency with the COMSOL operation result, and the effectiveness of an iterative operation mode related to joule heating of a line is verified.

The thermal analysis method for the double-sided PCB structure is mainly based on a Fourier series temperature analytical solution and a finite volume numerical value dispersion method, and can realize that only the metal layer and the corresponding surface area of the insulating layer are dispersed according to a two-dimensional circuit diagram, but not the whole PCB structure is dispersed; by carrying out numerical value discrete approximation on the heat balance equation of the metal layer, the linear relation between the temperatures of the metal layer units is utilized to approximately represent the heat diffusion inside the metal layer units and the heat conduction process between the metal layer layers through the metal via holes; and the heating rate of the device is calculated in the coupling solution by counting the thermal resistance parameters of the device.

The method can compensate and correct the thermal boundary condition related to the assumption of the contact area of the device, correct the homogeneity assumption of the insulating layer of the PCB substrate and the thermal boundary condition including the surface heat flux density distribution, and finally obtain the average temperature of the core of the device and the temperature distribution of other hot spot areas. Programming realization and simulation calculation of the method are carried out in an MATLAB environment, and the theoretical accuracy of the modeling method is verified by comparing with the COMSOL modeling operation result.

While the embodiments of the invention have been described in detail in connection with the drawings, it should not be construed as limiting the scope of the invention. Various modifications and changes may be made by those skilled in the art without inventive step within the scope of the appended claims.

Claims

1. A steady-state thermal analysis method for a double-sided PCB structure is characterized by comprising the following steps:

s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation, wherein the equation comprises the following steps:

wherein T is the amount of temperature change from ambient temperature or the temperature at ambient temperature of 0 ℃, q_iu(x, y) and q_id(x, y) is a heat flow density distribution function introduced by the upper and lower surfaces, h_uAnd h_dHeat transfer coefficients of the upper and lower surfaces of the insulating layer, k_iIs the thermal conductivity of the insulating layer, L_x、L_y、L_inThe lengths of the insulating layer in the x direction and the y direction and the thickness of the insulating layer in the z direction are respectively measured under a Cartesian coordinate system;

wherein H_uAnd H_dDeriving intermediate parameters of the process for the analytic solution; when h is generated_dIs zero, C₁,C_γn，mComprises the following steps:

multiplying both sides of the equation by cos (. beta.)_nx)cos(μ_my) and the integral is carried out in the surface heating area, according to the characteristics of the trigonometric function, the integral after multiplying different series numbers is zero, therefore, the left side of the equation only contains cos²(β_nx)cos²(μ_my) is not zero after integration, and the right side of the equation is used for each heating unit q in the heating area_iu,iIntegrating and then summing, and obtaining after operationCoefficient C₁And C_n,m：

solving the temperature distribution of the designated area is the thermal resistance row vector R_u(x,y,z)Expanding into a thermal resistance matrix R according to the coordinates of the corresponding distribution points of the region_uWherein each row is corresponding to a thermal resistance row vector for solving the temperature at a certain point, and then R is added_uMultiplying the temperature distribution vector by the corresponding heat flow density vector q to obtain a temperature distribution vector of the specified area, namely the heat effect generated by the corresponding q in the specified area; finally, superposing Fourier series analytic solutions in a matrix form corresponding to theta and eta of the specified region to obtain the temperature distribution of the region;

wherein, the distribution vectors of the heat flux density transmitted into the insulating layer by the upper and lower surface metal layers are q_M,u、q_M,dWherein q is_uAnd q is_dRespectively the heat flux density distribution vector q transmitted between the upper and lower surface heating components and the covering region thereof_uComprising q_iu(x, y) heat flux density distribution after dispersion, q_dComprising q_id(x, y) a heat flow density distribution after discretization; r_Mu,uIs q_M,uAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Mu,dIs q_M,dAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Md,uIs q_M,uAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_Md,dq_M,dAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_iu,uIs q_uAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_iu,dIs q_dAt T_uA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_id,uIs q_uAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_id,dIs q_dAt T_dA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; due to q_uAnd q is_dThe heat flow transmitted through the surface of the metal layer is counted in the heat diffusion matrix equation corresponding to the metal layer, and only q is counted in the analytic solution matrix equation_iu(x, y) and q_id(x, y) heat flow transferred between the upper and lower surface heat-generating components associated therewith and the surface areas of the insulating layer directly covered thereby, R_iu,u、R_iu,d、R_id,uAnd R_id,dWill include a partial all-zero column vector;

s5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and calculating potential distribution in the metal layer, so as to further calculate current density distribution and Joule heating distribution of the metal layer, and calculate the Joule heating distribution into a metal layer thermal diffusion matrix equation, and perform iterative calculation with temperature distribution to calculate the temperature distribution of the surface of each layer of the PCB under the condition of calculating Joule heating of the circuit;

2. The steady-state thermal analysis method for the double-sided PCB structure as recited in claim 1, wherein the step S2 includes the steps of calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the thermal boundary condition of the substrate insulating layer, and obtaining the corrected insulating layer temperature analytic solution matrix equation expression, which includes:

q_h,u＝q_u+h_uT_J,u

Wherein, T_J,uA temperature distribution vector for an area occupied by the upper surface device; because the covered area of the component does not exist in the range h between the component and the external environment_uThe calculated heat exchange under the condition, but the part of the heat exchange is taken into account in the analytic solution derivation process of the step S1, so the occupied area of the component is determined according to h_uCompensating for the multi-calculated heat exchange; in addition, can benefitThe thermal resistance parameter provided by the component is used for deducing the heat transfer between the relevant component and the PCB surface, and the specific deduction process is as follows:

the two formulas are combined to obtain:

q_h,u＝C_P,uq_J,u-C_T,uT_J,u

Wherein, T_J,dThe temperature distribution vector of the coverage area of the lower surface heating device is obtained; r_iu,uC_P,uq_J,uThe term is a vector of known constants, and T_uIncluding T_J,uAll elements of the vector can be transformed by adding zero to the matrix coefficients in a rising dimension way to realize the use of T in the matrix equation_uInstead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_Ju,uT_u＝R_iu,uC_T,uT_J,u

T_u＝R_Mu,uq_M,u+R_Md,uq_M,d+R_id,uq_d+R_iu,uC_P,uq_J,u-C_Ju,uT_u

wherein the heat flux density vector of the compensated lower surface covered by the heating device is q_h,d(ii) a The density vector of the heating surface corresponding to the lower surface component is q_J,d；C_P,dFor taking into account the thermal resistance parameter of the device_J,dThe coefficient of (a); c_Jd,uFor realizing the temperature distribution T of the covering area of the lower surface heating device_J,dAnd T_dCoefficients of the correlation term equivalent transform;

wherein, C_Ju,d、C_Jd,dRespectively for realizing T in the lower surface analytical solution matrix equation_J,uAnd T_uRelated item and T_J,dAnd T_dThe correlation term is equivalent to the transformed coefficients.

3. The steady-state thermal analysis method for the double-sided PCB structure of claim 1, wherein the step S3 approximates the discrete metal layer thermal equilibrium equation based on the finite volume method, generalizes the metal layer thermal diffusion matrix equation, and includes the following steps:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_c,uThe joule heating rate distribution per unit volume of the upper surface metal layer; t is_c,uIs the upper surface metal cell temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_u,uIs T_c,uThe longitudinal heat flow density of the upper surface of the corresponding unit can be the density of heat flow transferred between the corresponding unit and the heating device; q. q.s_d,uIs T_c,uThe longitudinal heat flux density transmitted between the lower surface of the corresponding unit and the insulating layer; t is_E,uIs T_c,uThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_W,uIs T_c,uThe temperature of the corresponding unit west side adjacent metal unit; t is_N,uIs T_c,uThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_S,uIs T_c,uThe temperature of the adjacent metal unit on the south side of the corresponding unit; psi_Nv,zIs the unit thermal resistance of the metal via; t is_c,dvIs the temperature of the lower surface metal layer unit connected by the via holeDegree;

wherein q is_c,dThe joule heating rate distribution per unit volume of the lower surface metal layer; t is_c,dIs the lower surface metal unit temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_d,dIs T_c,dThe longitudinal heat flux density of the lower surface of the corresponding unit can be the density of heat flux transferred between the corresponding unit and the heating device; q. q.s_u,dIs T_c,uThe longitudinal heat flux density transmitted between the upper surface of the corresponding unit and the insulating layer; t is_E,dIs T_c,dThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_W,dIs T_c,dThe temperature of the corresponding unit west side adjacent metal unit; t is_N,dIs T_c,dThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_S,dIs T_c,dThe temperature of the adjacent metal unit on the south side of the corresponding unit; t is_c,uvThe temperature of the lower surface metal layer unit connected through the via hole is shown;

wherein q is_MJ,uAnd q is_MJ,dRespectively corresponding unit joule heating vectors of the upper surface metal layer and the lower surface metal layer; h is_topThe heat transfer coefficient between the surface of the component and the outside is taken as the coefficient; m_MV,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; m_MV,dThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; m_V,dDiscrete thermal equilibrium equation from top surface metal layer unitThe temperature coefficient of the middle corresponding lower surface metal layer unit is formed; m_V,uThe coefficient corresponding to the temperature of the upper surface metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; due to q_h,u、q_h,dThe compensated heat flow directly transmitted through the surface of the insulating layer is already counted in the analytic solution matrix equation, and only the heat flow transmitted between the heating device and the covering metal layer part of the heating device is counted in the lumped matrix equation of heat diffusion, so C_M,u、C_M,dWill also include partial all-zero column vectors;

4. The steady-state thermal analysis method for the double-sided PCB structure of claim 1, wherein the step S4 is implemented for constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the modified insulation layer temperature analytic solution matrix equation and the metal layer numerical value discrete-based thermal diffusion matrix equation, and comprises:

5. A method for steady-state thermal analysis of a dual-sided PCB structure according to claim 1, wherein in step S5, a discrete metal layer current continuity equation is approximated based on a finite volume method, a lumped matrix equation describing a linear relationship between unit potentials of the metal layer is obtained by calculation and induction, and a potential distribution in the metal layer is obtained, so that a current density distribution and a joule heating distribution of the metal layer can be further obtained, and the joule heating distribution can be included in a metal layer thermal diffusion matrix equation, and an iterative calculation can be performed between the joule heating distribution and the temperature distribution, so as to obtain a temperature distribution of a surface where each layer of the PCB is located under the condition of the joule heating of the circuit, including:

wherein, V_Pc、V_PE、V_PW、V_PNAnd V_PSPc and P are respectively_E、P_W、P_NAnd the potential at PS; wherein σ_Pe、σ_Pw、σ_PnAnd σ_PsElectrical conductivities at Pe, Pw, Pn and Ps, respectively;

where ρ is the resistivity, ρ₀For resistivity at initial ambient temperature, Δ T_MTemperature rise of discrete metal units compared to ambient temperature, α_TAlpha for copper, for a material dependent temperature coefficient_TIs positive, i.e. its resistanceThe rate increases with increasing temperature;

the conductivity of the midpoint of the cross section is considered to be the conductivity corresponding to the equivalent resistance of the synthesized metal units at the left side, the right side, the upper side and the lower side, namely the reciprocal of the sum of the resistivities of the two units, sigma_PsConductivity σ Using C1 cell and C2 cell_C1And σ_C2Represents:

for the metal units at the line boundary positions, the analyzed control volume units are non-square, and the corresponding approximate discrete current continuity equation is as follows:

during the PCB design process, 45 degree inclined lines are designed, but a square grid can be reconstructed, and the discrete approximate current continuity equation of a square control volume unit with a quarter cell area and Ptr2 as a vertex at the line boundary is as follows:

while the control volume units at the line corners are not directly equivalent to the actual line corners, the control volume units with the cross sections of Pmn3, Pmw3 and Pms3 surrounding the boundary vertex Pb3 have the length d of the cross section of the lower Pms3 in the sawtooth boundary_cA/2, while the length in the actual line is d_cTherefore, such corners need to be individually processed in the thermal analysis method to account for the approximate current continuity equation corresponding to their actual boundaries:

wherein, V₁And V₂Respectively corresponding to the summary vectors of all the discrete unit top potential with unknown potential of the upper surface metal layer and the lower surface metal layer; m_σ,1And M_σ,2Respectively V in the upper surface matrix equation₁With V in the lower surface matrix equation₂The coefficient matrix of (a), is formed by conductivity parameters; k_p,1And K_p,2Known terms in two matrix equations are respectively included, and corresponding terms of known endpoint potential or section current are included; m_{σV,1_2}Is V in the upper surface matrix equation₂The coefficient matrix of (a) is formed by conductivity parameters corresponding to the lower surface metal layer unit connected by the via hole and the upper surface metal unit discrete heat balance equation, for V₂In cells not connected by vias, M_{σV,1_2}The corresponding column in (1) is an all-zero column; m_{σV,2_1}For V in the lower surface matrix equation₁The coefficient matrix of (a) is formed by conductivity parameters corresponding to the upper surface metal layer unit connected by the via in the discrete thermal equilibrium equation of the lower surface metal unit, for V₁In cells not connected by vias, M_{σV,2_1}The corresponding column in (1) is an all-zero column; the above matrix equation can therefore also incorporate the approximate discrete current continuity equation for metal cells not connected to a via; v can be obtained by simultaneously solving the two matrix equations₁And V₂；

After the peak potential distribution in the line is found, each square can be foundPotentials at the center of the cell and at the center of the cross-section around the cell were determined by finite volume method using a small diamond region around the center of the cell C2 to determine the potential V at the center of the cell C2_C2：

Wherein, V_C1Potential at the center of cell C1;

then C1 corresponds to a total current density of:

wherein, the third-order square matrix in the equation is composed of operation formulas corresponding to the coordinates of three vertexes under a Cartesian coordinate system, S_ΔIs a triangular area, and a coefficient a corresponding to the potential linear distribution expression can be obtained₀、a₁And a₂Then, the transverse current density J within the triangle is known_x,ΔAnd longitudinal current density J_y,ΔComprises the following steps:

wherein σ_ΔConductivity corresponding to the triangular region; on the premise of obtaining the current density corresponding to the metal cell, the joule heating rate per unit volume can be obtained, for example, for cell C1, the joule heating rate per unit volume q can be obtained_C1Can be expressed as:

q_C1＝E_C1·J_C1＝(J_x,C1 ²+J_y,C1 ²)/σ_C1＝J_C1 ²/σ_C1

wherein E is_C1The electric field strength at the center of cell C1;

the Joule heating rate is related to the conductivity, and the conductivity is influenced by the temperature, so when the Joule heating of the circuit is counted, iterative operation between the temperature distribution and the Joule heating rate distribution is required, when the error of the iterative operation is smaller than a set error, the temperature distribution and the Joule heating rate distribution can be approximately calculated, and the error can be set as the maximum value of the allowable unit temperature change in two iterative calculations; q. q.s_C1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the upper and lower surface metal layers_MJ,uAnd q is_MJ,d。