CN111723507B

CN111723507B - Multilayer PCB structure steady-state thermal analysis method

Info

Publication number: CN111723507B
Application number: CN202010573239.3A
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: CN111723507A

Abstract

The invention discloses a steady-state thermal analysis method for a multilayer PCB structure, which comprises the following steps: s1, constructing a steady-state heat balance equation of the PCB substrate insulating layer, and obtaining a matrix expression of the analytic solution of the substrate insulating layer; s2, correcting the boundary condition of the substrate insulating layer, and taking the heating rate of the component as a known heat source condition to be recorded into a matrix expression of an analytic solution; s3, calculating and summarizing to obtain a metal layer thermal diffusion lumped matrix equation, and calculating the thermal conductivity of the metal via hole in the thermal diffusion lumped matrix equation; s4, constructing a thermal analysis coupling equation set of the multilayer PCB structure; and S5, calculating and summarizing to obtain a lumped matrix equation describing the linear relation between the unit potentials of the metal layer, calculating the potential distribution in the metal layer, further calculating the current density distribution and the Joule heating distribution 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 condition of the Joule heating of the circuit.

Description

Multilayer 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 multilayer PCB structure.

Background

The PCB is not only a carrier of an electronic device, but also a main medium for heat dissipation of components, the PCB-level thermal analysis can realize the functions of assisting in optimizing the PCB design and analyzing the reliability of an electronic system, and if the PCB-level thermal analysis is combined with chip-level and packaging-level thermal analysis programs for simulation, the objective prediction of the hot spots in the integrated circuit can be realized under the condition of considering the actual PCB structure, circuit layout and thermal boundary conditions. The heat conduction and heat dissipation design of the PCB circuit can be assisted by establishing a PCB local circuit or a whole equivalent thermal analysis model in thermal analysis software, on one hand, a more practical heat dissipation circuit layout or heat dissipation structure can be designed, on the other hand, the heating condition of the device can be mastered in time, and the purpose of finally assisting and optimizing the PCB design is achieved.

At present, Flotherm software under the Mentor flag and Sigrity powerDC software under the Cadence flag are mainly used for directly realizing the PCB thermal analysis function, belong to Electronic Design Automation (EDA) design software, and can be used for assisting the design work of the PCB. Meanwhile, Computer Aided Engineering (CAE) software ANSYS Icepak, COMSOL Multiphysics and other CFD software for thermal analysis can also realize modeling and thermal simulation of PCBs and devices.

However, related EDA or CAE software aiming at PCB thermal analysis is not available at home at present. The occurrence of a complete export ban (including components, software, technical services and the like) on Chinese communication in the united states in 2018 reminds Chinese people of the need to develop electronic auxiliary design software with independent intellectual property rights.

The Yotanovich m.m. of the university of luo, canada, and the like, have proposed a method for calculating the lumped equivalent thermal conductivity of PCBs in the early stage, and the method is generally applicable to thermal analysis of PCBs with more uniform copper-clad distribution of metal layers, and can achieve the purpose of simplifying the operation amount, but the method may generate larger errors.

Byron Blackmore of Mentor company proposes to deduce the distribution of local equivalent thermal conductivity according to the local copper coating rate of each layer of a PCB, and then carries out thermal analysis on the whole PCB, which is also the basis of the algorithm of the currently widely applied PCB thermal analysis software Mentor Flotherm, and through the comparison of the simulation calculation and the temperature measurement experiment result of a PCB product, although the error of some hot spots exceeds 5%, the precision can be improved by improving the discrete unit density of each layer of the PCB. Furthermore, Flotherm is expensive to use, and has a key authorization (Lisence) cost of 60 to 200 million RMB in China.

M.b.dourouz of Cisco, inc, proposed a method of deriving the thermal conductivity of its locally discretized orthogonal anisotropy (x, y, z directions in cartesian coordinates) from the wiring diagram of the PCB and modeling the PCB using ANSYS Icepak, which by comparison indicated that the conventional method of analyzing the equivalent longitudinal and lateral uniform thermal conductivity of the PCB as a whole is less versatile. Meanwhile, the thermal analysis precision can be improved by improving the discrete resolution of the heat flow density distribution; and it was found that ignoring the line joule heating may underestimate the error of the device temperature rise by approximately 10% at a ratio of the device heating rate to the line joule heating rate on the PCB of approximately 12: 1.

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 multilayer PCB structure aiming at the defects in the prior art so as to solve the problem of large thermal analysis error of the conventional PCB structure.

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

a multi-layer PCB structure steady state thermal analysis method, comprising:

s1, constructing a steady-state thermal balance equation of the PCB substrate insulating layer, and solving Fourier series analytic solution coefficients based on boundary conditions of the steady-state thermal balance equation to obtain a matrix expression of the analytic solution of the substrate insulating layer;

s2, calculating the thermal effect of the component in the analytic solution based on the thermal resistance parameters, 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 expression of the analytic solution;

s3, approximating a heat balance equation of the discrete metal layer based on a finite volume method, calculating and inducing to obtain a metal layer heat diffusion lumped matrix equation, and adding the heat conduction of the metal via hole in the heat diffusion lumped matrix equation;

s4, constructing a thermal analysis coupling equation set of the multilayer PCB structure based on the assumption of longitudinal discrete distributed thermal resistance of the prepreg layer and according to the analytic solution matrix of the insulating layer of the substrate, the thermal diffusion lumped matrix equation of the metal layer, the current density distribution and the Joule heat distribution;

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 that current density distribution and Joule heating distribution of the metal layer can be further calculated, the Joule heating distribution can be calculated into a metal layer thermal diffusion matrix equation, iterative calculation is conducted between the Joule heating distribution and temperature distribution, and temperature distribution of the surface where each layer of the PCB is calculated on the premise of Joule heating of the circuit is calculated.

Preferably, step S1 is to construct a steady-state thermal equilibrium equation of the PCB substrate insulating layer, and solve the fourier series analytic solution coefficients based on the boundary conditions of the steady-state thermal equilibrium equation to obtain a matrix expression of the analytic solution of the substrate insulating layer, including:

the multilayer PCB structure is formed by mutually laminating a PCB substrate with metal circuits covered on the upper and lower surfaces and a prepreg for bonding, namely the number of metal layers in the PCB and the number of substrate layers are in a relation of 2: 1; and assuming that the temperature distribution of the metal layer is the same as that of the corresponding area of the surface of the insulating layer in contact with the metal layer;

constructing a steady-state heat balance equation corresponding to a substrate insulating layer in the PCB and calculating to obtain a boundary condition of an analytic solution:

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 of longitudinal transfer of the upper and lower surfaces, h_uAnd h_dRespectively, the heat transfer coefficients of the upper and lower surfaces of the insulating layer of the substrate, k_iIs the thermal conductivity of the insulating layer of the substrate, L_x、L_y、L_inThe length of the insulating layer of the substrate in the x and y directions and in the z direction respectivelyA thickness in a direction; for PCB top substrate, h_dIs zero; for PCB bottom substrate, h_uIs zero; and for other substrates inside the PCB, h_uAnd h_dAre all zero;

decomposing the thermal equilibrium equation and the boundary condition by adopting a superposition principle:

wherein theta is a superposition component variable of a substrate insulating layer heat balance equation solution; eta is a superposition component variable of a substrate insulating layer heat balance equation solution;

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;

the fourier series corresponding to η is obtained by the same method as the fourier series corresponding to θ, and the fourier series solution is substituted into the lower surface corresponding to θ (z ═ L)_in) The 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_γnmThe 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, 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 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 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;

because the thickness of the insulating layer of the substrate is far greater than the thickness of the metal layer and the prepreg layer, and the thermal conductivity of the prepreg layer is less than one thousandth of that of copper, the thermal diffusion effect of the prepreg layer is far lower than that of the copper-clad circuit metal layer bonded on the upper surface and the lower surface of the prepreg layer, and the transverse thermal diffusion effect can be approximately ignored compared with the metal layer; therefore, for the thermal analysis method of the multilayer structure, assuming that the transverse heat conduction in the prepreg layer is neglected, the prepreg layer is approximately dispersed into longitudinal distributed thermal resistance, wherein the longitudinal thermal resistance corresponding to the metal via hole in the prepreg layer is also included;

the matrix equation form of the analytic solution vector of the temperature distribution of the upper surface and the lower surface of the nth layer of substrate insulating layer in the PCB is as follows:

t is adopted because the temperature distribution analytical solution of the 2n-1 layer and the 2n layer metal layer is also included at the same time_2n-1And T_2nAs the temperature distribution vector sign of the corresponding surface; wherein q is_M,u,nAnd q is_M,d,nThe distribution vectors of the longitudinal heat flux density transmitted by the upper surface metal layer and the lower surface metal layer of the insulating layer of the nth substrate are obtained; r_M,2n-1、R_M,2nAre each q_M,u,n、q_M,d,nAt T_2n-1、T_2nA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_{M,2n_2n-1}Is q_M,d,nAt T_2n-1Thermal resistance system corresponding to thermal effect generated in limited areaThe number matrix relates to the thermal effect generated by the longitudinal transfer heat flow between the 2 n-th metal layer and the n-th substrate insulating layer on the 2 n-1-th metal layer, so that the suffix is 2n _2 n-1; r_{M,2n-1_2n}Is q_M,u,nAt T_2nA thermal coefficient of resistance matrix corresponding to the thermal effect is generated in the limited area, namely the thermal effect generated by the longitudinal transmission heat flow between the 2n-1 th metal layer and the n-th substrate insulating layer on the 2 n-th metal layer is related, so that the lower patch is 2n-1_2 n; q. q.s_pre,n-1The heat flux density distribution vector is longitudinally transmitted between the nth layer substrate and the nth-1 layer substrate through the nth-1 layer prepreg layer; q. q.s_pre,nThe heat flux density distribution vector is longitudinally transmitted between the nth layer of substrate and the (n + 1) th layer of substrate through the nth layer of prepreg layer; r_I,2n-1And R_{I,2n_2n-1}Are each q_pre,n-1And q is_pre,nAt T_2n-1A thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_{I,2n-1_2n}And R_I,2nAre each q_pre,n-1And q is_pre,nAt T_2nA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area;

the matrix equation of the analytic solution vector of the temperature distribution of the upper surface and the lower surface of the insulating layer of the first layer of the substrate on the top of the PCB is as follows:

T₁＝R_M,1q_M,u,1+R_{M,2_1}q_M,d,1-R_{I,2_1}q_pre，1+R_iu,1q_u

T₂＝R_{M,1_2}q_M,u,1+R_M,2q_M,d,1-R_I,2q_pre,1+R_iu,2q_u

wherein q is_uIs the heat flux density distribution vector, q, transferred between the upper surface heating component and the covering region_uIncludes q_iu(x, y) a heat flow density distribution after discretization; r_iu,1And R_iu,2Are each q_uAt T₁And T₂The defined area produces a matrix of thermal resistivity corresponding to thermal effects due to q_uThe 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_iuCorrelationThe upper surface heating element and the surface area of the insulating layer directly covered by the upper surface heating element are in heat transfer, so that R_iu,1And R_iu,2Will include a partial all-zero column vector; the definitions of other temperature vectors T, heat flux density distribution vectors q and matrix coefficients R can refer to the definitions of relevant vectors and matrix coefficients in the analytic solution matrix equation of the temperature distribution of the upper surface and the lower surface of the insulating layer of the nth layer of substrate;

the matrix equation of the analytic solution vector of the temperature distribution of the upper surface and the lower surface of the Nth layer of the substrate insulating layer at the bottom of the PCB is as follows:

T_2N-1＝R_M,2N-1q_M,u,N+R_{M,2N_2N-1}q_M,d,N+R_I,2N-1q_pre,N-1+R_id,2N-1q_d

T_2N＝R_{M,2N-1_2N}q_M,u,N+R_M,2Nq_M,d,N+R_{I,2N-1_2N}q_pre,N-1+R_id,2Nq_d

wherein q is_dThe heat flux density vector transferred between the heating component on the lower surface and the covering area thereof comprises q_id(x, y) a heat flow density distribution after discretization; r_id,2N-1And R_id,2NAre each q_dAt T_2N-1And T_2NA thermal resistivity matrix corresponding to the thermal effect produced in the defined area, and R_iu,1And R_iu,2Similarly, R_id,2N-1And R_id,2NWill also include partial all-zero column vectors; the definitions of other temperature vectors T, heat flux density distribution vectors q and matrix coefficients R can refer to the definitions of relevant vectors and matrix coefficients in the analytic solution matrix equation of the temperature distribution of the upper surface and the lower surface of the insulating layer of the nth layer of substrate; p_J,dIs the vector of the heating rate of the lower surface heating component in unit volume.

Preferably, the step S2 includes calculating the thermal effect of the component in the analytic solution based on the thermal resistance parameter, modifying the boundary condition of the substrate insulating layer, and calculating the heat generation rate of the component as the known heat source condition into a matrix expression of the analytic solution, including:

the vector of the heat flux density of the area covered by the heating device on the upper surface, which is counted for compensation, is 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; to the area occupied by the components according to h_uCompensating for the multi-calculated heat exchange; the thermal resistance parameters provided by the components are used for deducing the heat transfer between the related components and the PCB surface, and an analytic solution is taken into account, specifically comprising the following steps:

the two formulas are combined to obtain:

wherein P is the time when the device heating is negligible_J,uAnd q is_J,uIs zero; the temperature distribution of the surface unit of the upper surface element device 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；

With q_h,uReplacement of PCB top substrate insulation layer surface temperature distribution analytic solution T₁Q in (1)_u：

Wherein, C_P,uTo take into account device heatPost-parametric mapping q_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); r_iu,1C_P,uq_J,uThe term is a vector of known constants, and T₁Including 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₁Instead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_Ju,1T₁＝R_iu,1C_T,uT_J,u

wherein, C_Ju,1To realize T_J,uAnd T₁And (3) equivalently transforming the coefficients of the related terms to further obtain:

T₁＝R_M,1q_M,u,1+R_{M,2_1}q_M,d,1-R_{I,2_1}q_pre，1+R_iu,1C_P,uq_J,u-C_Ju,1T₁

similarly, the compensated T for thermal boundary conditions is taken into account₂The corresponding analytical solution matrix equation expression may be transformed into:

T₂＝R_{M,1_2}q_M,u,1+R_M,2q_M,d,1-R_I,2q_pre,1+R_iu,2C_P,uq_J,u-C_Ju,2T₁

when there are devices on the lower surface, the above derivation and equivalent transformation method can be applied to T_2N-1And T_2NAnd correspondingly correcting heat flow transmitted between the lower surface device and the surface unit in an analytic solution matrix equation:

in the same way, T_2NThe corresponding analytical solution matrix equation will be modified to:

T_2N＝R_{M,2N-1_2N}q_M,u,N+R_M,2Nq_M,d,N+R_{I,2N-1_2N}q_pre,N-1+R_id,2NC_P,dq_J,d-C_Jd,2NT_2N

wherein the density vector of the heating surface corresponding to the lower surface component is q_J,d(ii) a And C_P,uAnd C_T,uSimilarly, C_P,dAnd C_T,dFor counting the corresponding q after the thermal resistance parameter of the lower surface device is counted_J,dAnd T_J,dThe coefficient of (a); and C_Ju,1Function of (C) is similar to the definition_Jd,2N-1、C_Jd,2NRespectively for realizing the temperature distribution T of the covering area of the lower surface heating device_J,dAnd T_2N-1And T_2NThe correlation term is equivalent to the transformed coefficients.

Preferably, the step S3 approximates the equation of the thermal equilibrium of the discrete metal layer based on the finite volume method and takes into account the thermal conductivity of the metal via, including:

the discrete heat balance equation corresponding to the 2n-1 layer metal layer unit on the upper surface of the nth layer PCB substrate is as follows:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_c,2n-1Is the volume Joule heating rate, T, of a certain unit in the 2n-1 th metal layer_c,2n-1Is the temperature of the metal unit, q_u,2n-1The density of heat flow longitudinally transmitted on the upper surface of the unit can be the density of heat flow transmitted between the heating device or the density of heat flow transmitted through a metal through hole in a prepreg layer contacted with the heating device; including its passage of q_d,2n-1Is the heat flux density, T, of the unit lower surface in longitudinal transfer_E,2n-1Is the temperature of the metal cell adjacent to the east side of the cell, T_W,2n-1Is the temperature, T, of the west adjacent metal unit of the unit_N,2n-1Is the temperature of the adjacent metal cell on the north side of the cell, T_S,2n-1The temperature of the adjacent metal unit on the south side of the unit; d_cIs the unit side length of the structured discrete metal unit; d_mIs the thickness of the metal layer; psi_Nv,zIs the unit thermal resistance of the metal via; t is_c,2nThe temperature of the lower surface 2 nth layer metal layer unit connected through the via hole is shown;

similarly, a discrete thermal balance equation corresponding to the 2 nth layer metal layer unit on the lower surface of the nth layer PCB substrate can be obtained:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_c,2nIs the volume Joule heating rate, T, of a certain unit in the 2 n-th metal layer_c,2nIs the temperature of the metal unit, q_u,2nDensity of heat flow transferred longitudinally to the upper surface of the unit, q_d,2nThe density of heat flow longitudinally transmitted by the lower surface of the unit can be the density of heat flow transmitted between the heating device or the density of heat flow transmitted through a metal through hole in a prepreg layer contacted with the heating device; t is_E,2nIs the temperature of the metal cell adjacent to the east side of the cell, T_W,2nIs the temperature, T, of the west adjacent metal unit of the unit_N,2nIs the temperature of the adjacent metal cell on the north side of the cell, T_S,2nThe temperature of the adjacent metal cell on the south side of the cell.

Preferably, the step S3 generalizes the discrete thermal balance equation corresponding to the metal layer unit to a metal layer thermal diffusion lumped matrix equation, which includes:

summarizing the discrete heat balance equations corresponding to all units of the same metal layer to obtain a lumped matrix equation representing the heat diffusion of the metal layer, wherein the equation form of the metal layer on the upper surface and the lower surface of the nth layer substrate in the PCB is as follows:

M_MV,2n-1T_2n-1-M_{V,2n-1_2n}T_2n＝q_MJ，2n-1+C_M,2n-1q_pre,n-1-q_M,u,n

M_MV,2nT_2n-M_{V,2n_2n-1}T_2n-1＝q_MJ，2n-C_M,2nq_pre,n-q_M,d,n

wherein q is_MJ,2n-1And q is_MJ,2nRespectively corresponding unit joule heating vectors of the 2n-1 th metal layer and the 2n th metal layer; m_MV,2n-1From discrete heat balance equation of layer 2n-1 metal layer unit corresponding to temperature of metal layer unitForming a coefficient; m_{V,2n-1_2n}The coefficient corresponding to the temperature of the 2n layer metal layer unit in the discrete heat balance equation of the 2n-1 layer metal layer unit is formed; m_MV,2nThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the 2 n-th metal layer unit is formed; m_{V,2n_2n-1}The coefficient corresponding to the temperature of the 2n-1 layer metal layer unit in the discrete heat balance equation of the 2n layer metal layer unit is formed; q. q.s_pre,n-1And q is_pre,nThe heat flux density of the substrate surface metal layer unit and the heat flux density of the substrate insulating layer surface unit are respectively included, and the matrix coefficient C corresponding to the heat flux density is obtained_M,2n-1And C_M,2nWill include a partial all-zero column vector; t is_2n-1And T_2nIncluding the temperature of all units of the 2n-1 th and 2n th metal layers and the temperature of the surface unit of the substrate insulating layer covered by partial non-metal layers, so that the temperature at M is higher than that of the surface unit of the substrate insulating layer_MV,2n-1、M_{V,2n-1_2n}、M_MV,2nAnd M_{V,2n_2n-1}Will include a partial all-zero column vector;

for the layer 1 substrate on the top of the PCB, the lumped matrix equation of the thermal diffusion corresponding to the metal layers on the upper and lower surfaces is:

M_MV,1T₁-M_{V,1_2}T₂＝q_MJ，1+C_M,uq_h,u-q_M,u,1

M_MV,2T₂-M_{V,2_2}T₁＝q_MJ，2-C_M,2nq_pre,1-q_M,d,1

due to q_h,uThe 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,uWill also include partial all-zero column vectors; q is to be_h,uSubstituting the expression into a thermal diffusion matrix equation corresponding to the layer 1 metal layer to obtain:

C_Mu,1to achieve T in the above equation_J,uAnd T₁Coefficients of the correlation term equivalent transform;

for the Nth layer of the substrate at the bottom of the PCB, the lumped matrix equation of the thermal diffusion corresponding to the metal layers on the upper and lower surfaces is as follows:

M_MV,2N-1T_2N-1-M_{V,2N-1_2N}T_2N＝q_MJ，2N-1+C_M,2N-1q_pre,N-1-q_M,u,N

M_MV,2NT_2N-M_{V,2N_2N-1}T_2N-1＝q_MJ，2N+C_M,dq_h,d-q_M,d,N

due to 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,dWill also include partial all-zero column vectors; q is to be_h,dSubstituting the expression into a thermal diffusion matrix equation corresponding to the 2N-th metal layer to obtain:

C_Md,2Nto achieve T in the above equation_J,dAnd T_2NThe correlation term is equivalent to the transformed coefficients.

Preferably, in step S4, based on the assumption that the thermal resistance is distributed in a longitudinal discrete manner in the prepreg layer, and according to the analytic solution matrix of the insulating layer of the substrate and the lumped matrix equation of thermal diffusion of the metal layer, a thermal analysis coupling equation set of the multilayer PCB structure is constructed, including:

temperature distribution T of upper and lower surfaces of n-th prepreg layer_2nAnd T_2n+1That is, the relationship between the temperature distributions of the two substrate surfaces corresponding to their bonding can be approximately expressed as:

T_2n-T_2n+1＝R_pre,nq_pre,n

wherein q is_pre,nThe thermal resistance vector of the longitudinal discrete distribution of the nth prepreg layer is shown;

thus, the system of thermal analysis coupling equations for an N-layer PCB substrate bonded by N-1 prepreg layers, i.e., a multilayer PCB structure with 2N (N >1) metal layers, will have the following form:

wherein, T₁～T_2N，q_M,u,1～q_M,u,N，q_M,d,1～q_M,d,N，q_pre,1～q_pre,N-1The vector is an unknown vector, and other vectors and matrix coefficients are known quantities;

the temperature corresponding to the metal unit in the temperature distribution coupling solution vector is the temperature rise of the metal unit compared with the ambient temperature, and the method can be used for calculating the conductivity influenced by the temperature, further calculating the Joule heating distribution of the metal layer under the temperature distribution, and calculating the iterative operation of the temperature distribution under the condition of Joule heating of the metal layer.

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 further calculate the joule heating distribution into a metal layer thermal diffusion matrix equation, and iteratively calculate the joule heating distribution 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 surrounded by the sections of the metal layer unit with the vertex Pc as the center and Pn, Pw, Ps and Pe:

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;

the finite volume method is applied to the control volume taking the vertex as the center, the conductivity of the midpoint of the section is the conductivity corresponding to the equivalent resistance synthesized by the metal units at the left side, the right side or the upper side and the lower side, namely the reciprocal of the sum of the resistivity of the two units, sigma_PsConductivity σ Using C1 cell and C2 cell_C1And σ_C2Represents:

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 electrical conductivity of each element of the metal layer is therefore calculated on the basis of the temperature distribution obtained in claim 6;

similarly, the conductivity corresponding to Pc can be approximately found by applying the conductivities corresponding to the four midpoints surrounding Pc:

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:

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 process of designing the PCB, a 45-degree inclined line is often designed, a square grid is used to reconstruct the 45-degree inclined line, and a discrete approximate current continuity equation corresponding to a small square control volume unit with a quarter cell area and a Ptr2 as a vertex at the boundary of the line is as follows:

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 the current continuity equation corresponding to the triangle control cell surrounding Ptr2 after filling the sawtooth gap of the boundary;

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 σ_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 case for a 45 degree diagonal line and its corresponding control volume range;

for the metal layer units connected by the via holes, the corresponding approximately discrete current continuity equation is obtained according to the finite volume method taking the vertex as the center, and if the vertices of the 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 an approximate discrete current continuity equation corresponding to each metal layer is obtained, and the current continuity matrix equation corresponding to the 2n-1 th metal layer has the following form:

M_σ，2n-1V_2n-1-M_{σV，2n-1_2n-2}V_2n-2-M_{σV，2n-1_2n}V_2n＝K_p,2n-1

wherein, V_2n-2、V_2n-1And V_2nRespectively summarizing vectors of all unknown top potential of the 2n-2 nd layer, the 2n-1 th layer and the 2n-2 nd layer metal layer; m_σ,2n-1For V in the matrix equation_2n-1The coefficient matrix is formed by conductivity parameters; m_{σV,2n-1_2n-2}And M_{σV,2n-1_2n}Is V in the matrix equation_2n-2And V_2nBy a 2 n-th coefficient matrix connected to the 2n-1 th metal layer unit through a via-the conductivity parameters of the 2 layers corresponding to the 2 n-th metal layer unit are formed for V not connected by the via_2n-2And V_2nCell top potential of (1), M_{σV,2n-1_2n-2}And M_{σV,2n-1_2n}The corresponding column in the matrix is an all-zero column, so the matrix equation can also be included in an approximate discrete current continuity equation of the metal unit which is not connected with the through hole; k_p,2n-1The parameter vector comprises known items in all matrix equations;

from this, a set of matrix equations can be listed that approximately describe the continuity of the current for all metal layers:

wherein, V₁～V_2NThe summary vectors of the top potential of the discrete units with unknown potentials respectively corresponding to the first to 2N metal layers can be obtained by simultaneously solving the 2N matrix equations; m_σ1～M_σ2NV in the current continuity matrix equation corresponding to the first to 2N metal layers respectively₁～V_2NCoefficient matrix of, M_{σV,1_2}～M_{σV,2N-1_2N}Respectively correspond to V in the first 2N-1 matrix equations₂～V_2NA coefficient matrix of (a); m_{σV,2_1}～M_{σV,2N_2N-1}Respectively correspond to V in the subsequent 2N-1 matrix equations₁～V_2N-1A coefficient matrix of (a); k_p,1～K_p,2NThe parameter vector comprises known items in all matrix equations;

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 can be obtained, and the potential V at the center of the cell C2 can be obtained by using a finite volume method in a small diamond region around the center of the cell C2_C2：

Wherein, V_Pvt4Potential of cell C2 at vertex Pvt4(ii) a 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 longitudinal current density corresponding to the center of each cell, for example, the lateral current density J corresponding to the center of the cell C1, are approximated_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_ΔFor the area of the triangle, the coefficient a corresponding to the linear distribution expression of the potential 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 corresponding to the metal cell can be obtained, and for the 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_C1The thickness of the metal layer can be multiplied to obtain the heating surface density of the corresponding metal layer surface unit, and the heating surface density of all units of the n-th metal layer forms the joule heating surface density vector q of the modified layer_MJ,n。

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

the temperature distribution of the surface of the insulating layer is obtained by calculating temperature analytic solution, so that the problem that the whole structure needs to be discretized by numerical analysis methods such as a finite element method is avoided.

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 thermal analysis process for PCB structures.

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 an embodiment of the present application, referring to fig. 1, the method for steady-state thermal analysis of a multilayer PCB structure of the present solution includes:

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 that current density distribution and Joule heating distribution of the metal layer can be further calculated, the Joule heating distribution can be calculated into a metal layer thermal diffusion matrix equation, iterative calculation is performed between the Joule heating distribution and temperature distribution, and temperature distribution of the surface where each layer of the PCB is calculated on the premise of Joule heating of a circuit is calculated;

s6, realizing the general program structure and flow of multilayer PCB structure thermal analysis.

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

Firstly, thermal analysis assumption needs to be carried out on the structural characteristics of the PCB:

the PCB structure is formed by mutually overlapping and pressing an insulating layer and a metal layer, the insulating layer is divided into two types, one type is a core substrate layer (the thickness is 1-2 mm), and the insulating layer plays a role of mechanical support; one type is a prepreg pre-layer (thickness less than 0.1mm) which serves as an adhesive between the substrate layers, the pre-layer material being an unhardened substrate material, and thus the thermal conductivity of the substrate layer and the pre-layer can be considered to be the same. The thermal conductivity of the common FR-4 glass epoxy resin insulating layer material in the PCB is 0.3W/mK, the crystallization temperature is about 120 ℃, and when the temperature is exceeded, the mechanical property and the insulating property of FR-4 are obviously reduced. The PCB applied in a high temperature environment may use ceramic or polyimide as an insulating layer material.

While the material of the metal layer is predominantly copper (thermal conductivity of about 390W/mK), typically 0.7mil (1mil is 25.4 μm, i.e., one thousandth of an inch, 0.7mil also corresponds to 0.5oz, in the electronics industry copper clad metal layers are 1oz, i.e., the thickness corresponding to a 1 ounce copper weight laid uniformly over a 1 square foot area), 1.37mil (35 μm, corresponding to 1oz), or other custom made thickness. In the multi-layer PCB structure, the upper surface and the lower surface of each layer of substrate are covered with metal layers, and then a plurality of substrates are bonded and laminated by using a prefabricated layer material to form the multi-layer PCB structure. The upper surface and the lower surface of the PCB are also covered with a soldemask solder resist coating and a silk screen printing layer for marking devices, printing trademarks and other information. The method mainly realizes the thermal analysis of the multilayer structure consisting of the insulating layer and the metal layer, and 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 the method; since solder resist coatings are thin (typically 1mil) and have low thermal conductivity (0.25W/mK), their equivalent thermal resistance can be taken into account in the heat transfer coefficient between the PCB surface and the outside world.

The PCB has the main heating source which is a component with relatively large power loss, and in order to realize large-scale production, the PCB product mainly adopts automatic welding, namely the packaging of the component is also mainly surface-mounted packaging, because the bottom of a part of the components with larger heat productivity is provided with a pad area or the base of the part of the components is provided with the pad area, the heat diffusion of the part of the components is directly realized through a metal layer, and other surface-mounted components with smaller heat productivity transmit heat through pins on one hand and through the occupied PCB surface on the other hand. The existence of the components can enable the heat exchange between the PCB and the outside to not only pass through the surface of the multilayer board but also pass through the surface of the device, the heat conduction relation related to the heat resistance parameter of the device is utilized in the method of the invention to account for the heat effect related to the device, and the assumed condition that the whole surfaces of the upper insulating layer and the lower insulating layer of the multilayer board and the outside exchange heat under a certain heat transfer coefficient is corrected, for example, the temperature analytic solution or the discrete heat balance equation corresponding to the PCB surface unit with the heating device can be supplemented and corrected.

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 with 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 insulating layer is a main constituent material of the PCB, but since it is not a heat source, the temperature distribution inside the insulating layer is not a main solution target.

Therefore, in the method, the temperature distribution analytical solution of the surface of the insulating layer is used for giving a temperature distribution expression of the insulating layer contacting with the metal layer and the heating device, namely the temperature distribution expression is used as one of the coupling variables of the numerical discrete analysis and the analytical solution; meanwhile, the heat conduction inside the metal layer and between the metal layers through the via holes is calculated by using a numerical value dispersion method based on a finite volume method, so as to provide an approximate steady-state thermal diffusion matrix equation containing the temperature distribution of the metal layer, wherein the heat flow density distribution actually transmitted to the insulating layer by the metal layer is contained, and the distribution is also used as one of coupling variables; and finally solving the solution of the coupling variable through a simultaneous analytical 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 area of the four sides of the PCB is generally much smaller than that 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 in the method. The method is therefore not suitable for thermal analysis of PCBs of the type where heat transfer is achieved primarily through the side and where the side has a plug-in portion, etc. in the application.

The following detailed description of the thermal analysis method is made under the above-described assumption:

step S1, constructing a steady-state thermal balance equation of the PCB insulating layer, and solving Fourier series analytic solution coefficients based on boundary conditions of the steady-state thermal balance equation to obtain a matrix expression of the analytic solution of the insulating layer, wherein the matrix expression specifically 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 of longitudinal transfer of the upper and lower surfaces, h_uAnd h_dRespectively, the heat transfer coefficients of the upper and lower surfaces of the insulating layer of the substrate, k_iIs the thermal conductivity of the insulating layer of the substrate, L_x、L_y、L_inThe length of the substrate insulating layer in the x direction and the y direction and the thickness of the substrate insulating layer in the z direction are respectively measured under a Cartesian coordinate system; for PCB top substrate, h_dIs zero; for PCB bottom substrate, h_uIs zero; and for other substrates inside the PCB, h_uAnd h_dAre all zero;

theta is a superposition component variable of a substrate insulating layer heat balance equation solution; eta, the superposition component variable of the substrate insulating layer heat balance equation solution;

only the solving process of theta corresponding to the Fourier series is described, and the solving method of eta corresponding to the Fourier series is the same as that of the solving method of the Fourier series, and is not described in detail. Substituting fourier series solution into the lower surface corresponding to θ (z ═ L)_in) The 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:

only description of h is provided herein_dSolving process of Fourier series when not zero, h_dThe solving method of zero-time Fourier series is the same as that of zero-time Fourier series, and is not described in detail. Substituting the series solution into the thermal boundary condition of the upper surface to obtain:

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 columns of heat flux density are constructedVector q_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; 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;

because the thickness (>0.5mm) of the insulating layer of the substrate is far larger than that of the metal layer (generally smaller than 0.1mm) and the prepreg layer (about 0.1mm), and the thermal conductivity of the prepreg layer is less than one thousandth of that of copper, the thermal diffusion effect of the prepreg layer is generally far lower than that of the copper-clad circuit metal layer bonded to the upper surface and the lower surface of the prepreg layer, and the transverse thermal diffusion effect can be approximately ignored compared with the metal layer. Therefore, for the thermal analysis method of the multilayer structure, one basic assumption to be adopted is to ignore the transverse heat conduction in the prepreg layer and approximately disperse the prepreg layer into longitudinal distributed thermal resistance, wherein the longitudinal thermal resistance corresponding to the metal via hole in the prepreg layer is also included;

the matrix equation form of the analytic solution vector of the temperature distribution of the upper surface and the lower surface of the nth layer of substrate insulating layer in the PCB is shown as follows, and T is adopted because the analytic solution of the temperature distribution of the 2n-1 layer and the 2n layer of metal layer is also included at the same time_2n-1And T_2nTemperature distribution vector sign as the corresponding surface:

wherein q is_M,u,nAnd q is_M,d,nIs as followsThe upper surface metal layer and the lower surface metal layer of the insulating layer of the n layers of substrates and the longitudinal heat flux density distribution vector transmitted by the upper surface metal layer and the lower surface metal layer are distributed; r_M,2n-1、R_M,2nAre each q_M,u,n、q_M,d,nAt T_2n-1、T_2nA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_{M,2n_2n-1}Is q_M,d,nAt T_2n-1A thermal coefficient of resistance matrix corresponding to the thermal effect is generated in the limited area, namely the thermal effect generated by the longitudinal transfer heat flow between the 2 n-th metal layer and the n-th substrate insulating layer on the 2 n-1-th metal layer is related, so that the suffix is 2n _2 n-1; r_{M,2n-1_2n}Is q_M,u,nAt T_2nA thermal coefficient of resistance matrix corresponding to the thermal effect is generated in the limited area, namely the thermal effect generated by the longitudinal transmission heat flow between the 2n-1 th metal layer and the n-th substrate insulating layer on the 2 n-th metal layer is related, so that the lower patch is 2n-1_2 n; q. q.s_pre,n-1The heat flux density distribution vector is longitudinally transmitted between the nth layer substrate and the nth-1 layer substrate through the nth-1 layer prepreg layer; q. q.s_pre,nThe heat flux density distribution vector is longitudinally transmitted between the nth layer of substrate and the (n + 1) th layer of substrate through the nth layer of prepreg layer; r_I,2n-1And R_{I,2n_2n-1}Are each q_pre,n-1And q is_pre,nAt T_2n-1A thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_{I,2n-1_2n}And R_I,2nAre each q_pre,n-1And q is_pre,nAt T_2nA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area;

the matrix equation form of the analytic solution vector of the temperature distribution of the upper surface and the lower surface of the insulating layer of the first layer substrate on the top of the PCB is shown as follows:

wherein q is_u(W/m²) Is the heat flux density distribution vector, q, transferred between the upper surface heating component and the covering region_uIncludes q_iu(x, y) a heat flow density distribution after discretization; r_iu,1(℃m²/W) and R_iu,2(℃m²W) are each q_uAt T₁And T₂The defined area produces a matrix of thermal resistivity corresponding to thermal effects due to q_uThe 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_iuHeat flow transferred between the associated top surface heat generating component and the surface area of the insulating layer directly covered thereby, R_iu,1And R_iu,2Will include a partial all-zero column vector; the definitions of other temperature vectors T, heat flux density distribution vectors q and matrix coefficients R can refer to the definitions of relevant vectors and matrix coefficients in the analytic solution matrix equation of the temperature distribution of the upper surface and the lower surface of the insulating layer of the nth layer of substrate;

the matrix equation form of the analytic solution vector of the temperature distribution of the upper surface and the lower surface of the Nth layer of the substrate insulating layer at the bottom of the PCB is as follows:

wherein q is_d(W/m²) The heat flux density vector transferred between the heating component on the lower surface and the covering area thereof comprises q_id(x, y) a heat flow density distribution after discretization; r_id,2N-1(℃m²/W) and R_id,2N(℃m²W) are each q_dAt T_2N-1And T_2NA thermal resistivity matrix corresponding to the thermal effect produced in the defined area, and R_iu,1And R_iu,2Similarly, R_id,2N-1And R_id,2NWill also include partial all-zero column vectors; the definitions of other temperature vectors T, heat flux density distribution vectors q and matrix coefficients R can refer to the definitions of relevant vectors and matrix coefficients in the analytic solution matrix equation of the temperature distribution of the upper surface and the lower surface of the insulating layer of the nth layer of substrate; p_J,d(W/m³) Is a unit volume heating rate vector of the lower surface heating component;

step S2, calculating the thermal effect of the component in the analytic solution based on the thermal resistance parameters, correcting the thermal boundary condition of the top or bottom substrate insulating layer, and obtaining a corrected insulating layer temperature analytic solution matrix equation, which specifically comprises:

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

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 an analytic solution is taken into account, and the specific deduction process is as follows:

the two formulas are combined to obtain:

Wherein, 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); r_iu,1C_P,uq_J,uThe term is a vector of known constants, and T₁Including 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₁Instead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_Ju,1T₁＝R_iu,1C_T,uT_J,u (16)

T₁＝R_M,1q_M,u,1+R_{M,2_1}q_M,d,1-R_{I,2_1}q_pre，1+R_iu,1C_P,uq_J,u-C_Ju,1T₁ (17)

T₂＝R_{M,1_2}q_M,u,1+R_M,2q_M,d,1-R_I,2q_pre,1+R_iu,2C_P,uq_J,u-C_Ju,2T₁ (18)

wherein the heat flux density vector of the compensated lower surface covered by the heating device is q_h,d；

T_2N＝R_{M,2N-1_2N}q_M,u,N+R_M,2Nq_M,d,N+R_{I,2N-1_2N}q_pre,N-1+R_id,2NC_P,dq_J,d-C_Jd,2NT_2N (20)

Step S3, approximating a discrete metal layer and its thermal equilibrium equation based on a finite volume method of second order Taylor series approximation and a structured square grid, and taking into account the heat transfer between the metal layer and the heat generating device and the heat conduction of the metal via hole, which specifically includes:

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.

For counting in the metal layerThermal diffusion will primarily use a finite volume method to discretize the thermal equilibrium equations described above, and a structured square grid to discretize the metal layer, 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 of the conduction between the central unit and the four surrounding units is q 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)²) The joule heat inside the cell is represented 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 the linear processing equation 20 can be approximated, and then the heat dissipation balance equation can be finally realized, and an approximate linear relationship between the temperatures of the 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 heat diffusion of the metal layer, and when the units are connected with other adjacent metal layer units through metal via holes, the q is also induced based on the heat resistance of the metal via holes_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.

S3.3, inducing a discrete heat balance equation corresponding to the metal layer unit into a metal layer heat diffusion lumped matrix equation;

wherein q is_MJ,2n-1And q is_MJ,2nRespectively corresponding unit joule heating vectors of the 2n-1 th metal layer and the 2n th metal layer; m_MV,2n-1The coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the 2n-1 layer metal layer unit is formed; m_{V,2n-1_2n}The coefficient corresponding to the temperature of the 2n layer metal layer unit in the discrete heat balance equation of the 2n-1 layer metal layer unit is formed; m_MV,2nThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the 2 n-th metal layer unit is formed; m_{V,2n_2n-1}The coefficient corresponding to the temperature of the 2n-1 layer metal layer unit in the discrete heat balance equation of the 2n layer metal layer unit is formed; q. q.s_pre,n-1And q is_pre,nThe heat flux density of the substrate surface metal layer unit and the heat flux density of the substrate insulating layer surface unit are respectively included, and the matrix coefficient C corresponding to the heat flux density is obtained_M,2n-1And C_M,2nWill include a partial all-zero column vector; t is_2n-1And T_2nIncluding the temperature of all units of the 2n-1 th and 2n th metal layers and the temperature of the surface unit of the substrate insulating layer covered by partial non-metal layers, so that the temperature at M is higher than that of the surface unit of the substrate insulating layer_MV,2n-1、M_{V,2n-1_2n}、M_MV,2nAnd M_{V,2n_2n-1}Will include a partial all-zero column vector;

C_Md,2Nto achieve T in the above equation_J,dAnd T_2NCoefficients of related term equivalent transform

T_2n-T_2n+1＝R_pre,nq_pre,n (37)

wherein, T₁～T_2N，q_M,u,1～q_M,u,N，q_M,d,1～q_M,d,N，q_pre,1～q_pre,N-1For unknown vectors, other vectors and matrix coefficients are known.

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

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:

similarly, an approximate discrete equation corresponding to the current continuity equation in the form of the surface integral can be obtained:

wherein, V_PW、V_PNAnd V_PSAre each P_W、P_NThe potential at PS;

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

the relationship of conductivity to temperature can be represented by the following equation:

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:

for metal cells at line boundary positions, it may happen that the analyzed control volume cell is non-square, such as the control cell No. 2 (which is followed by the smaller square control cell No. 1) surrounded by dark black multi-dotted 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 upper left line boundary in 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 shows the current continuity equation and boundaryThe current continuity equations for the triangular control cells surrounding Ptr2 after the saw tooth gap is filled are the same. 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 Pms3 below the control volume units 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 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, for example, the current continuity matrix equation corresponding to the 2n-1 th metal layer has the following form:

M_σ，2n-1V_2n-1-M_{σV，2n-1_2n-2}V_2n-2-M_{σV，2n-1_2n}V_2n＝K_p,2n-1 (52)

wherein, V_2n-2、V_2n-1And V_2nRespectively summarizing vectors of all unknown top potential of the 2n-2 nd layer, the 2n-1 th layer and the 2n-2 nd layer metal layer; m_σ,2n-1For V in the matrix equation_2n-1The coefficient matrix is formed by conductivity parameters; m_{σV,2n-1_2n-2}And M_{σV,2n-1_2n}Is V in the matrix equation_2n-2And V_2nThe coefficient matrix is composed of conductivity parameters corresponding to the 2n-2 layer and 2n layer metal layer units which are connected with the 2n-1 layer metal layer unit through the via hole, and for V which is not connected through the via hole_2n-2And V_2nCell top potential of (1), M_{σV,2n-1_2n-2}And M_{σV,2n-1_2n}The corresponding column in the matrix is an all-zero column, so the matrix equation can also be included in an approximate discrete current continuity equation of the metal unit which is not connected with the through hole; k_p,2n-1All known terms in the matrix equation, such as corresponding terms describing known endpoint potentials or section currents, are included as parameter vectors;

wherein, V₁～V_2NRespectively corresponding to the first layerThe summary vector of all the discrete unit top potential with unknown potential from the 2N layer metal layer can be obtained by simultaneously solving the 2N matrix equations; m_σ,1～M_σ,2NV in the current continuity matrix equation corresponding to the first to 2N metal layers respectively₁～V_2NCoefficient matrix of, M_{σV,1_2}～M_{σV,2N-1_2N}Respectively correspond to V in the first 2N-1 matrix equations₂～V_2NA coefficient matrix of (a); m_{σV,2_1}～M_{σV,2N_2N-1}Respectively correspond to V in the subsequent 2N-1 matrix equations₁～V_2N-1A coefficient matrix of (a); k_p,1～K_p,2NAll known terms in the matrix equation are included for the parameter vector.

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_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;

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, wherein the method 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 unit at the oblique line boundary, the electric potential can be assumed to be linearly distributed in the region, so as to find 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, 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 σ_x,ΔConductivity corresponding to the triangular region;

on the premise of obtaining the current density corresponding to the metal unit, the corresponding joule heating per unit volume can be obtainedRate, e.g. for cell C1, joule heating per unit volume rate q_C1Can be expressed as:

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

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, 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 joule heating surface density vector q of the layer-changed layer by the heating surface density of all the units of the n-th metal layer_MJ,n。

Step S6, a general program structure and a flow for implementing thermal analysis of a multilayer PCB structure, which specifically include:

the general structure and flow of programming a PCB structure thermal analysis program by using the method of the invention are shown in FIG. 7, the program can be divided into an initialization program and a main program, if an analysis model is newly established in the initialization program, the number of metal layer layers can be firstly set, then a circuit diagram is introduced, in the process of analyzing the circuit by the program, the identification and the position confirmation of each component part such as a circuit, a device, a bonding pad and the like are required, and then a coefficient matrix required by summarizing a coupling equation set can be generated, wherein the coefficient matrix comprises a coefficient matrix required by calculating, analyzing and solving a thermal resistance coefficient matrix and a numerical value discrete equation and is counted in the thermal conductivity of the metal layer and a through hole, the thermal. If the existing model is modified, the step of setting/modifying the system electrical parameters and thermal parameters can be directly carried out, and finally, the data are stored.

The main procedure for implementing the thermal analysis method can include two parts of coupling calculation and iterative calculation, wherein the coupling calculation calculates the thermal analysis result according to the coupling equation system by calling each coefficient matrix and the known heating density vector. 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 uniform temperature distribution of the circuit, then the coupling solution of the corresponding temperature distribution is calculated, the Joule heating density distribution of the new circuit is calculated, 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.

The invention is based on the Fourier series temperature analytic solution and the finite volume numerical value discrete method (FVM), which realizes the discrete of the metal layer and the corresponding surface area of the insulating layer, but not the discrete of the whole PCB structure; the calculation of the current density distribution of the metal circuit and the joule heat distribution of the electric-thermal coupling is realized through numerical value discrete approximation of a current continuity equation of the metal layer; the numerical value dispersion approximation is carried out on the heat balance equation of the metal layer, so that the linear relation among the temperatures of the metal layer units is utilized to approximately represent the internal heat diffusion process, and the heat conduction between the metal layers through the through holes is calculated; by taking the thermal resistance parameters of the device into account, the heating rate of the device in the coupling solution is realized, the compensation and the correction of the assumed thermal boundary conditions related to the contact area of the device are realized, the thermal boundary conditions of the whole PCB surface are maintained, and finally the average temperature of the core of the device can be solved according to the temperature;

while the embodiments of the invention have been described in detail in connection with the accompanying drawings, it is not intended to limit 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 multilayer PCB structure is characterized by comprising the following steps:

s1, a steady-state thermal balance equation of the PCB substrate insulating layer is constructed, Fourier series analytic solution coefficients are solved based on boundary conditions of the steady-state thermal balance equation, and a matrix expression of analytic solutions of the substrate insulating layer is obtained, wherein the matrix expression comprises the following steps:

the fourier series corresponding to η is obtained by the same method as the fourier series corresponding to θ, and the fourier series solution is substituted into the lower surface z corresponding to θ, where z is 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 d is_qAfter discretizing surface heat source gridLength of grid cell, integral field Sq_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:

t is adopted because the temperature distribution analytical solution of the 2n-1 layer and the 2n layer metal layer is also included at the same time_2n-1And T_2nAs the temperature distribution vector sign of the corresponding surface; wherein q is_M,u,nAnd q is_M,d,nThe distribution vectors of the longitudinal heat flux density transmitted by the upper surface metal layer and the lower surface metal layer of the insulating layer of the nth substrate are obtained; r_M,2n-1、R_M,2nAre each q_M,u,n、q_M,d,nAt T_2n-1、T_2nA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_{M,2n_2n-1}Is q_M,d,nAt T_2n-1A thermal coefficient of resistance matrix corresponding to the thermal effect is generated in the limited area, namely the thermal effect generated by the longitudinal transfer heat flow between the 2 n-th metal layer and the n-th substrate insulating layer on the 2 n-1-th metal layer is related, so that the suffix is 2n _2 n-1; r_{M,2n-1_2n}Is q_M,u,nAt T_2nA thermal coefficient of resistance matrix corresponding to the thermal effect is generated in the limited area, namely the thermal effect generated by the longitudinal transmission heat flow between the 2n-1 th metal layer and the n-th substrate insulating layer on the 2 n-th metal layer is related, so that the lower patch is 2n-1_2 n; q. q.s_pre,n-1The heat flux density distribution vector is longitudinally transmitted between the nth layer substrate and the nth-1 layer substrate through the nth-1 layer prepreg layer; q. q.s_pre,nThe heat flux density distribution vector is longitudinally transmitted between the nth layer of substrate and the (n + 1) th layer of substrate through the nth layer of prepreg layer; r_I,2n-1And R_{I,2n_2n-1}Are each q_pre,n-1And q is_pre,nAt T_2n-1A thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; r_{I,2n-1_2n}And R_I,2nAre each q_pre,n-1And q is_pre,nAt T_2nA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area;

T₁＝R_M,1q_M,u,1+R_{M,2_1}q_M,d,1-R_{I,2_1}q_pre，1+R_iu,1q_u

T₂＝R_{M,1_2}q_M,u,1+R_M,2q_M,d,1-R_I,2q_pre,1+R_iu,2q_u

wherein q is_uIs the heat flux density distribution vector, q, transferred between the upper surface heating component and the covering region_uIncludes q_iu(x, y) a heat flow density distribution after discretization; r_iu,1And R_iu,2Are each q_uAt T₁And T₂The defined area produces a matrix of thermal resistivity corresponding to thermal effects due to q_uThe 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_iuHeat flow transferred between the associated top surface heat generating component and the surface area of the insulating layer directly covered thereby, R_iu,1And R_iu,2Will include a partial all-zero column vector; the definitions of other temperature vectors T, heat flux density distribution vectors q and matrix coefficients R can refer to the definitions of relevant vectors and matrix coefficients in the analytic solution matrix equation of the temperature distribution of the upper surface and the lower surface of the insulating layer of the nth layer of substrate;

T_2N-1＝R_M,2N-1q_M,u,N+R_{M,2N_2N-1}q_M,d,N+R_I,2N-1q_pre,N-1+R_id,2N-1q_d

T_2N＝R_{M,2N-1_2N}q_M,u,N+R_M,2Nq_M,d,N+R_{I,2N-1_2N}q_pre,N-1+R_id,2Nq_d

2. The steady-state thermal analysis method for the multi-layer PCB structure as recited in claim 1, wherein the step S2 includes calculating thermal effects of components in analytical solutions based on thermal resistance parameters, modifying substrate insulating layer boundary conditions, and calculating the heating rates of the components as known heat source conditions into a matrix expression of the analytical solutions, including:

q_h,u＝q_u+h_uT_J,u

the two formulas are combined to obtain:

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

C_Ju,1T₁＝R_iu,1C_T,uT_J,u

T₂＝R_{M,1_2}q_M,u,1+R_M,2q_M,d,1-R_I,2q_pre,1+R_iu,2C_P,uq_J,u-C_Ju,2T₁

3. The steady-state thermal analysis method for the multi-layer PCB structure of claim 1, wherein the step S3 approximates the heat balance equation of the discrete metal layer based on the finite volume method and takes into account the thermal conductivity of the metal via, comprising:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_c,2nIs the volume Joule heating rate, T, of a certain unit in the 2 n-th metal layer_c,2nIs the temperature of the metal unit, q_u,2nDensity of heat flow transferred longitudinally to the upper surface of the unit, q_d,2nThe heat flux density transmitted longitudinally to the lower surface of the unit can be transmitted to the heating deviceTransferring the density of the heat flow or transferring the density of the heat flow through a metal via in a prepreg layer contacted by the metal via; t is_E,2nIs the temperature of the metal cell adjacent to the east side of the cell, T_W,2nIs the temperature, T, of the west adjacent metal unit of the unit_N,2nIs the temperature of the adjacent metal cell on the north side of the cell, T_S,2nThe temperature of the adjacent metal cell on the south side of the cell.

4. The steady-state thermal analysis method for the multi-layer PCB structure of claim 1, wherein the step S3 generalizes the discrete thermal balance equation corresponding to the metal layer unit into a metal layer thermal diffusion lumped matrix equation, comprising:

M_MV,2n-1T_2n-1-M_{V,2n-1_2n}T_2n＝q_MJ，2n-1+C_M,2n-1q_pre,n-1-q_M,u,n

M_MV,2nT_2n-M_{V,2n_2n-1}T_2n-1＝q_MJ，2n-C_M,2nq_pre,n-q_M,d,n

wherein q is_MJ,2n-1And q is_MJ,2nRespectively corresponding unit joule heating vectors of the 2n-1 th metal layer and the 2n th metal layer; m_MV,2n-1The coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the 2n-1 layer metal layer unit is formed; m_{V,2n-1_2n}The coefficient corresponding to the temperature of the 2n layer metal layer unit in the discrete heat balance equation of the 2n-1 layer metal layer unit is formed; m_MV,2nThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the 2 n-th metal layer unit is formed; m_{V,2n_2n-1}The coefficient corresponding to the temperature of the 2n-1 layer metal layer unit in the discrete heat balance equation of the 2n layer metal layer unit is formed; q. q.s_pre,n-1And q is_pre,nThe heat flux density transmitted by the metal layer unit on the surface of the substrate and part of the heat flux insulated directly by the substrateThe heat flux density transmitted by the surface unit of the layer and the corresponding matrix coefficient C_M,2n-1And C_M,2nWill include a partial all-zero column vector; t is_2n-1And T_2nIncluding the temperature of all units of the 2n-1 th and 2n th metal layers and the temperature of the surface unit of the substrate insulating layer covered by partial non-metal layer, thereby M is_MV,2n-1、M_{V,2n-1_2n}、M_MV,2nAnd M_{V,2n_2n-1}Will include a partial all-zero column vector;

M_MV,1T₁-M_{V,1_2}T₂＝q_MJ，1+C_M,uq_h,u-q_M,u,1

M_MV,2T₂-M_{V,2_2}T₁＝q_MJ，2-C_M,2nq_pre,1-q_M,d,1

M_MV,2N-1T_2N-1-M_{V,2N-1_2N}T_2N＝q_MJ，2N-1+C_M,2N-1q_pre,N-1-q_M,u,N

M_MV,2NT_2N-M_{V,2N_2N-1}T_2N-1＝q_MJ，2N+C_M,dq_h,d-q_M,d,N

5. The steady-state thermal analysis method for the multilayer PCB structure according to claim 1, wherein the step S4 is implemented by constructing a thermal analysis coupling equation set of the multilayer PCB structure based on the assumption of the longitudinal discrete distributed thermal resistance of the prepreg layer and according to an analytic solution matrix of the substrate insulating layer and a lumped matrix equation of the thermal diffusion of the metal layer, and the method comprises the following steps:

T_2n-T_2n+1＝R_pre,nq_pre,n

wherein q is_pre,nLongitudinally and discretely distributing thermal resistance vectors for the nth prepreg layer;

thus, the system of thermal analysis coupling equations for an N-layer PCB substrate bonded by N-1 prepreg layers, i.e., a multilayer PCB structure with 2N metal layers, will have the following form:

6. The steady-state thermal analysis method for the multi-layer PCB structure of claim 5, wherein in step S5, the discrete metal layer current continuity equation is approximated based on the finite volume method, the lumped matrix equation describing the linear relationship between the unit potentials of the metal layer is obtained through calculation and induction, and the potential distribution in the metal layer is obtained, so as to further obtain the current density distribution and the Joule heating distribution of the metal layer, and the Joule heating distribution is included in the metal layer thermal diffusion matrix equation, and iterative calculation is performed between the Joule heating distribution and the temperature distribution, so as to obtain the temperature distribution of the surface where each layer of the PCB is located under the condition of the Joule heating of the circuit, comprising:

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 σ_PsConductance at Pe, Pw, Pn and Ps, respectivelyRate;

applying a finite volume method to the control volume with the vertex as the center, wherein the conductivity of the midpoint of the cross section is the conductivity corresponding to the equivalent resistance synthesized by the metal units at the left side, the right side, the upper side and the lower side, namely the reciprocal of the sum of the resistivity of the two units, sigma_PsConductivity σ Using C1 cell and C2 cell_C1And σ_C2Represents:

therefore, the conductivity of each unit of the metal layer is calculated according to the obtained temperature distribution;

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;

M_σ，2n-1V_2n-1-M_{σV，2n-1_2n-2}V_2n-2-M_{σV，2n-1_2n}V_2n＝K_p,2n-1

wherein, V_2n-2、V_2n-1And V_2nRespectively summarizing vectors of all unknown top potential of the 2n-2 nd layer, the 2n-1 th layer and the 2n-2 nd layer metal layer; m_σ,2n-1For V in the matrix equation_2n-1The coefficient matrix is formed by conductivity parameters; m_{σV,2n-1_2n-2}And M_{σV,2n-1_2n}Is V in the matrix equation_2n-2And V_2nThe coefficient matrix is composed of conductivity parameters corresponding to the 2n-2 layer and 2n layer metal layer units which are connected with the 2n-1 layer metal layer unit through the via hole, and for V which is not connected through the via hole_2n-2And V_2nCell top potential of (1), M_{σV,2n-1_2n-2}And M_{σV,2n-1_2n}The corresponding column in (1) is all-zero column, so the above matrix equation can also be included in the metal sheet without connecting via holeAn approximate discrete current continuity equation for the element; k_p,2n-1The parameter vector comprises known items in all matrix equations;

wherein, V₁～V_2NThe summary vectors of the top potential of the discrete units with unknown potentials respectively corresponding to the first to 2N metal layers can be obtained by simultaneously solving the 2N matrix equations; m_σ,1～M_σ,2NV in the current continuity matrix equation corresponding to the first to 2N metal layers respectively₁～V_2NCoefficient matrix of, M_{σV,1_2}～M_{σV,2N-1_2N}Respectively correspond to V in the first 2N-1 matrix equations₂～V_2NA coefficient matrix of (a); m_{σV,2_1}～M_{σV,2N_2N-1}Respectively correspond to V in the subsequent 2N-1 matrix equations₁～V_2N-1A coefficient matrix of (a); k_p,1～K_p,2NThe parameter vector comprises known items in all matrix equations;

Wherein, V_C1Potential at the center of cell C1;

then C1 corresponds to a total current density of:

wherein σ_ΔConductivity corresponding to the triangular region;

in obtaining the correspondence of metal unitsOn the premise of the current density of (1), the corresponding joule heating rate per unit volume can be obtained, and for the cell C1, the joule heating rate per unit volume q_C1Can be expressed as:

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