CN111723486B - Double-sided PCB structure steady-state thermal analysis method - Google Patents
Double-sided PCB structure steady-state thermal analysis method Download PDFInfo
- Publication number
- CN111723486B CN111723486B CN202010572678.2A CN202010572678A CN111723486B CN 111723486 B CN111723486 B CN 111723486B CN 202010572678 A CN202010572678 A CN 202010572678A CN 111723486 B CN111723486 B CN 111723486B
- Authority
- CN
- China
- Prior art keywords
- equation
- metal layer
- unit
- temperature
- heat
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Operations Research (AREA)
- Investigating Or Analyzing Materials Using Thermal Means (AREA)
Abstract
The invention discloses a steady-state thermal analysis method for a double-sided PCB structure, which comprises the following steps: s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and obtaining an insulating layer temperature analytic solution matrix equation; s2, correcting the boundary condition of the substrate insulating layer, and taking the heating rate of the component as the known heat source condition to be recorded into a matrix equation expression of an analytic solution; s3, approximating a heat balance equation of the discrete metal layer based on a finite volume method, and summarizing to obtain a metal layer heat diffusion lumped matrix equation; s4, constructing a thermal analysis coupling equation set of the double-sided PCB structure; and S5, calculating and inducing to obtain a lumped matrix equation describing the linear relation between the unit potentials of the metal layer, calculating the Joule heating distribution into the metal layer thermal diffusion matrix equation, and performing iterative calculation with the temperature distribution to obtain the temperature distribution of the surface of each layer of the PCB under the premise of calculating the Joule heating of the circuit. And S6, verifying the calculation accuracy of the thermal analysis method of the double-sided PCB structure by adopting a COMSOL modeling operation result.
Description
Technical Field
The invention belongs to the technical field of PCB thermal analysis, and particularly relates to a steady-state thermal analysis method for a double-sided PCB structure.
Background
The PCB is a main carrier of electronic components, has become an important component of electronic equipment, and the reliability of the PCB is related to the reliable operation of the whole electronic equipment. One of the main factors affecting the reliable operation of the PCB and the components thereon is the temperature. The design of the PCB structure ensures that the heating of the components on the PCB structure can be effectively diffused and can exchange heat with the external environment, namely, the temperature rise of the components is ensured to be within a safe range, and the PCB also becomes a main medium for heat dissipation of the components. Therefore, under the condition of considering the actual PCB structure, the circuit layout and the thermal boundary condition, the thermal analysis of the PCB can realize the effect of assisting in optimizing the PCB design and analyzing the reliability of the electronic system.
The PCB thermal simulation software widely applied in the market at present is FlothERM of Mentor company and Icepak of ANSYS company. The former is based mainly on approximate equivalent thermal conductivity related to the local copper coverage of the PCB in the algorithm principle, and the latter is based on the finite element method. The COMSOL and Sigrity PowerDC based on the finite element method can also perform thermal simulation analysis on the basis of completing PCB and device modeling.
Byron Blackmore of Mentor corporation in 2009 disclosed a method adopted by FloTHERM to derive local equivalent thermal conductivity from the local copper coverage of PCBs. According to the data announced by the company website, the occupation rate of Flotherm in the PCB thermal analysis market in the year is up to 88%. In the past decade, Flotherm has also been used by many national scholars to develop relevant application studies. However, a white paper disclosed in 2018 by Mentor corporation reveals that the method based on local copper coating rate adopted previously may have the problems of overestimating local thermal conductivity and underestimating actual temperature rise. The white paper indicates that another empirical algorithm is added in the Flotherm XT v3.2 version of software to analyze the thermal distribution of the PCB, the new algorithm is a solution method based on the local transverse thermal conductivity summarized after the analysis of 140 PCB products, and the thermal conductivity is still related to the copper coating rate; however, the white paper also indicates that empirical algorithms may underestimate actual thermal conductivity and thus overestimate actual temperature rise. Furthermore, Flotherm is expensive to use, and has a key authorization (Lisence) cost of 60 to 200 million RMB in China.
Finite Element Method (FEM) is a numerical discrete method widely used, and is characterized in that a continuous unknown function is approximated to a finite sub-function (generally in a polynomial form containing unknown variables which can be approximated to linear change in the discrete domain) effective in the discrete control domain in a discrete finite space or time, the finite sub-functions are substituted into an original continuous equation to establish the relation between corresponding variables of each discrete sub-domain, and then the approximate distribution of the unknown variables can be solved by using a given boundary condition. The software such as ANSYS Icepak and COMSOL Multiphysics based algorithms are used by many scholars to perform thermal simulation analysis of power devices, PCBs, and electronic systems.
The results of questionnaires on the use experiences of multi-physical field analysis software such as ANSYS, COMSOL and the like in the industry and the academia of the university of Necklands, Lester university, UK, Song Xueguan and the like of the university of the big connective workers in combination with the university of Necklands, show that more than 60% of users are satisfied with the operation accuracy of the software, but the same more than 60% of users consider that the operation speed of the software is slow. The CFD software based on the finite element method generally needs to disperse the whole structure to be analyzed in the simulation operation process, which may be one of the main reasons for the unsatisfactory operation speed. Meanwhile, the price of the software is expensive, taking COMSOL as an example, the price of the software is between 10 and 40 thousands according to the selected functional module.
From the development trend of the thermal analysis method of the PCB, the method for establishing the overall transverse and longitudinal equivalent thermal conductivity of the PCB in the early stage cannot be used for thermal analysis of the PCB as a general method due to poor structural applicability, otherwise, larger calculation errors can be brought. Therefore, methods for deriving equivalent thermal conductivity and finding local discretization orthotropic thermal conductivity based on the local copper cladding rate of each layer of the PCB are developed subsequently, and common characteristics of the methods are that errors caused by using the local equivalent thermal conductivity are expected to be reduced by discretizing the PCB structure and enhancing analysis of the local line distribution and the laminated board structure of the PCB.
Disclosure of Invention
The invention aims to provide a steady-state thermal analysis method of a double-sided PCB structure aiming at the defects in the prior art so as to solve the problem of large thermal analysis error of the prior PCB structure.
In order to achieve the purpose, the invention adopts the technical scheme that:
a steady-state thermal analysis method for a double-sided PCB structure comprises the following steps:
s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation;
s2, calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the boundary condition of the substrate insulating layer, and calculating the heating rate of the component as the known heat source condition into a matrix equation expression of an analytic solution;
s3, approximating a heat balance equation of a discrete metal layer based on a finite volume method, summarizing to obtain a metal layer heat diffusion lumped matrix equation, and counting the heat conduction between the metal layer and a surface heating device thereof and the heat conduction of the metal via hole;
s4, constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the corrected insulating layer temperature analytic solution matrix equation and the metal layer thermal diffusion matrix equation based on numerical value dispersion;
s5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and calculating potential distribution in the metal layer, so as to further calculate current density distribution and Joule heating distribution of the metal layer, and calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with temperature distribution to calculate temperature distribution of the surface where each layer of the PCB is positioned under the condition of Joule heating of a circuit;
and S6, verifying the calculation accuracy of the thermal analysis method of the double-sided PCB structure by adopting a COMSOL modeling operation result.
Preferably, step S1 is to construct a steady-state thermal balance laplace equation of the double-sided PCB insulating layer, and solve a fourier series analytic solution coefficient based on the boundary condition of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation, including:
constructing a stable heat balance Laplace equation of the double-sided PCB insulating layer, and obtaining boundary conditions of the periphery, the upper surface and the lower surface of the PCB insulating layer:
wherein T is the amount of temperature change from ambient temperature or the temperature at ambient temperature of 0 ℃, qiu(x, y) and qid(x, y) is a distribution function of the heat flux density introduced by the upper and lower surfaces, huAnd hdHeat transfer coefficients of the upper and lower surfaces of the insulating layer, kiIs the thermal conductivity of the insulating layer, Lx、Ly、LinThe lengths of the insulating layer in the x direction and the y direction and the thickness of the insulating layer in the z direction are respectively measured under a Cartesian coordinate system;
the expression form of Fourier series analytic solution corresponding to the steady-state heat balance equation is as follows:
wherein, C1,C2,Cn,m,Cγn,mIs a Fourier series coefficient, betan、μm、γn,mCharacteristic values in a Fourier series;
decomposing a heat balance equation and a boundary condition by adopting a heat effect superposition principle:
wherein theta and eta are superposition component variables of the insulating layer heat balance equation solution;
substituting Fourier series solution into lower surface z corresponding to theta as LinThe relation between partial coefficient and characteristic value in solution is obtained under the condition of thermal boundary, and when h is reacheddWhen not zero, C2And Cγn,mThe corresponding expression is:
wherein HuAnd HdDeriving intermediate parameters of the process for the analytic solution; when h is generateddIs zero, C1,Cγn,mComprises the following steps:
substituting the series solution into the thermal boundary condition of the upper surface to obtain:
multiplying both sides of the equation by cos (. beta.)nx)cos(μmy) and the integral is carried out in the surface heating area, according to the characteristics of the trigonometric function, the integral after multiplying different series numbers is zero, so that the left side of the equation only contains cos2(βnx)cos2(μmy) is not zero after integration, and the right side of the equation is used for each heating unit q in the heating areaiu,iIntegrating and then summing, and calculating to obtain a coefficient C1And Cn,m:
Wherein d isqIntegration Domain Sq for grid cell Length after discretization of surface Heat Source gridiu,iI.e. representing an incoming heat flow density of qiu,iSurface discrete unit area of δnAnd deltamTo resolve dimensionless intermediate parameters of the derivation process, qiu,1~qiu,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, qiu,1~qiu,N1I.e. the column vector q of the heat flow density is constructediuThen, 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 Ru(x,y,z)Expanding into a thermal resistance matrix R according to the coordinates of the corresponding distribution points of the regionuWherein each row is corresponding to a thermal resistance row vector for solving the temperature of a certain point, and then R is addeduMultiplying the temperature distribution vector by the corresponding heat flow density vector q to obtain a temperature distribution vector of the specified area, namely the heat effect generated by the corresponding q in the specified area; finally, superposing Fourier series analytic solutions in a matrix form corresponding to theta and eta of the specified region to obtain the temperature distribution of the region;
thus, the matrix equation corresponding to the analytic solution temperature distribution vector of the upper and lower surfaces of the insulating layer in the double-sided PCB will have the following form:
wherein, the distribution vectors of the heat flux density transmitted into the insulating layer by the upper and lower surface metal layers are qM,u、qM,dWherein q isuAnd q isdThe distribution vector q of heat flux density transferred between the upper and lower surface heating components and the covering region thereofuComprising qiu(x, y) heat flux density distribution after dispersion, qdComprising qid(x, y) a heat flow density distribution after discretization; rMu,uIs qM,uAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMu,dIs qM,dAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMd,uIs qM,uAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMd,dqM,dAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; riu,uIs quAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; riu,dIs qdAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rid,uIs quAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rid,dIs qdAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; due to quAnd q isdThe 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 equationiu(x, y) and qid(x, y) heat flow transferred between the associated top and bottom surface heat generating components and the surface area of the insulating layer directly covered thereby, Riu,u、Riu,d、Rid,uAnd Rid,dWill include a partial all zero column vector.
Preferably, step S2 includes calculating the thermal effect of the device based on the thermal resistance parameter of the device, and modifying the thermal boundary condition of the substrate insulating layer to obtain a modified insulating layer temperature analytic solution matrix equation expression, which includes:
setting the heat flux density vector of the covering area of the upper surface heating device as qh,u:
qh,u=qu+huTJ,u
Wherein, TJ,uA temperature distribution vector for an area occupied by the upper surface device; because the covered area of the component does not exist direct connection with the external environmentIs in huThe calculated heat exchange under the condition, but the part of the heat exchange is already taken into account in the analytic solution derivation process of the step S1, so the occupied area of the component is required to be determined according to huCompensating for the multiple calculated heat exchanges; in addition, the thermal resistance parameters provided by the components can be used for deducing the heat transfer between the related components and the PCB surface, and the specific deduction process is as follows:
the two formulas are combined to obtain:
qh,u=CP,uqJ,u-CT,uTJ,u
wherein the temperature distribution of the upper surface element device surface is TtopThe heat transfer coefficient of the surface to the outside is htopThe equivalent heat conduction corresponding to the heat sink existing on the surface thereof is also taken into accounttop(ii) a Mean temperature of the core is Tj(ii) a The average thermal resistance between the surface and the core is psijt(ii) a The average thermal resistance between the junction center and the bottom of the device is psijc(ii) a Equivalent thermal resistance of solder mask coating is psic,u(ii) a The heat flux density transferred from the core to the surface of the component is qjt(ii) a The heating volume density vector corresponding to the upper surface element device is PJ,u(ii) a The diagonal matrix composed of height parameters of the elements is CJ,u(ii) a The density vector of the heating surface corresponding to the upper surface element device is qJ,u;CP,uFor taking into account the thermal resistance parameter of the deviceJ,uThe coefficient of (d); cT,uFor taking into account the corresponding T after the thermal resistance parameter of the deviceJ,uThe coefficient of (a);
with qh,uAlternative upper surface temperature distribution analytical solution TuQ in (1)u:
Wherein, TJ,dThe temperature distribution vector of the coverage area of the lower surface heating device is obtained; riu,uCP,uqJ,uThe term is a vector of known constants, and TuIncluding TJ,uAll elements of the vector, so the matrix coefficient can be transformed by adding zero to the ascending dimension to realize the T in the matrix equationuInstead of TJ,uTo achieve the goal of unifying vectors, the equivalent transformation equation is as follows:
CJu,uTu=Riu,uCT,uTJ,u
wherein, CJu,uTo realize TJ,uAnd TuAnd (3) equivalently transforming the coefficients of the related terms to further obtain:
Tu=RMu,uqM,u+RMd,uqM,d+Rid,uqd+Riu,uCP,uqJ,u-CJu,uTu
similarly, when there is a device on the lower surface, the heat flow transferred between the device on the lower surface and the surface unit, which takes into account the analytical solution, is also corrected, i.e. q in the above formuladItem replacement:
wherein the heat flux density vector of the compensated lower surface covered by the heating device is qh,d(ii) a The density vector of the heating surface corresponding to the lower surface component is qJ,d;CP,dFor taking into account the thermal resistance parameter of the deviceJ,dThe coefficient of (a); and CJu,uFunction of (C) is similar to the definitionJd,uFor realizing the temperature distribution T of the coverage area of the lower surface heating deviceJ,dAnd TdCoefficients of the correlation term equivalent transform;
in the same way, the T corrected by the thermal boundary condition can be obtaineddCorresponding analytic solution matrix equation, and after correction, TuAnd TdCorresponding coupletThe analytical solution matrix equation becomes:
wherein, with CJu,uAnd CJd,uFunction of (C) is similar to the definitionJu,d、CJd,dRespectively for realizing T in the lower surface analytical solution matrix equationJ,uAnd TuRelated item and TJ,dAnd TdThe correlation term is equivalent to the transformed coefficients.
Preferably, the step S3 approximates the discrete metal layer heat balance equation based on the finite volume method, generalizes the metal layer heat diffusion matrix equation, and includes the heat conduction between the metal layer and the surface heat generating device thereof, and the heat conduction of the metal via, and includes:
dispersing the metal layer by adopting a structured square grid, and dispersing a heat balance equation of heat conduction in the metal layer based on second-order Taylor series approximation to obtain an upper surface discrete metal layer heat balance equation:
wherein k ismThe thermal conductivity of the metal layer; q. q.sc,uThe joule heating rate distribution per unit volume of the upper surface metal layer; t isc,uIs the upper surface metal cell temperature; dcThe unit side length of the structured and discrete metal unit is shown; dmIs the thickness of the metal layer; q. q.su,uIs Tc,uThe longitudinal heat flow density of the upper surface of the corresponding unit can be the density of heat flow transferred between the corresponding unit and the heating device; q. q.sd,uIs Tc,uThe longitudinal heat flux density transmitted between the lower surface of the corresponding unit and the insulating layer; t isE,uIs Tc,uThe temperature of the metal unit adjacent to the east side of the corresponding unit; t isW,uIs Tc,uThe temperature of the corresponding unit west side adjacent metal unit; t isN,uIs Tc,uThe temperature of the adjacent metal unit at the north side of the corresponding unit; t isS,uIs Tc,uCorrespondence sheetThe temperature of the adjacent metal unit on the side of the Yunan; psiNv,zIs the unit thermal resistance of the metal via; t isc,dvThe temperature of the lower surface metal layer unit connected through the via hole is shown;
similarly, a discrete heat balance equation corresponding to the lower surface metal layer unit can be obtained:
wherein q isc,dThe joule heating rate distribution per unit volume of the lower surface metal layer; t isc,dIs the lower surface metal unit temperature; dcThe unit side length of the structured and discrete metal unit is shown; dmIs the thickness of the metal layer; q. q.sd,dIs Tc,dThe longitudinal heat flux density of the lower surface of the corresponding unit can be the density of heat flux transmitted between the corresponding unit and the heating device; q. q.su,dIs Tc,uThe longitudinal heat flux density transmitted between the upper surface of the corresponding unit and the insulating layer; t isE,dIs Tc,dThe temperature of the adjacent metal unit at the east side of the corresponding unit; t isW,dIs Tc,dThe temperature of the corresponding unit west side adjacent metal unit; t isN,dIs Tc,dThe temperature of the adjacent metal unit at the north side of the corresponding unit; t isS,dIs Tc,dThe temperature of the adjacent metal unit on the south side of the corresponding unit; t isc,uvThe temperature of the lower surface metal layer unit connected through the via hole is shown;
and (3) inducing a thermal diffusion matrix equation of the upper and lower surface metal layers based on numerical value dispersion:
wherein q isMJ,uAnd q isMJ,dRespectively corresponding unit joule heating vectors of the upper surface metal layer and the lower surface metal layer; h istopThe heat transfer coefficient between the surface of the component and the outside is shown; mMV,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; mMV,dBy detachment of the metal layer unit from the lower surfaceThe coefficient corresponding to the temperature of the metal layer unit in the heat dissipation balance equation is formed; mV,dThe coefficient of the temperature of the lower surface metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; mV,uThe coefficient corresponding to the temperature of the upper surface metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; due to qh,u、qh,dThe compensated heat flow directly transmitted through the surface of the insulating layer is counted in an analytic solution matrix equation, and only the heat flow transmitted between the heating device and the part of the covering metal layer is counted in a lumped matrix equation of heat diffusion, so that CM,u、CM,dWill also include partial all-zero column vectors;
q is to beh,uAnd q ish,dSubstituting the expression of (a) into the thermal diffusion matrix equation to obtain:
wherein, CMu,uAnd CMd,dRespectively for realizing T in the above equationJ,uAnd TuRelated item and TJ,dAnd TdThe correlation term is equivalent to the transformed coefficients.
Preferably, in step S4, the method includes constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the modified insulating layer temperature analytic solution matrix equation and the metal layer numerical value discrete-based thermal diffusion matrix equation, where the method includes:
the thermal analysis coupling equation set of the double-sided PCB structure is as follows:
wherein, Tu,Td,qM,u,qM,dIf the vector is unknown, and other vector and matrix coefficients are known, the unknown temperature vector and the unknown heat flux density vector can be obtained by combining the four matrix equations.
Preferably, in step S5, the method approximates a discrete metal layer current continuity equation based on a finite volume method, calculates and summarizes a lumped matrix equation describing a linear relationship between unit potentials of the metal layer, and calculates a potential distribution in the metal layer, so as to further calculate a current density distribution and a joule heating distribution of the metal layer, and calculate the joule heating distribution into a metal layer thermal diffusion matrix equation, and perform iterative calculation with a temperature distribution to calculate a temperature distribution of a surface where each layer of the PCB is located under the condition of joule heating of the circuit, including:
based on second-order Taylor series approximation, obtaining an approximate discrete equation corresponding to a current continuity equation of a square area which is surrounded by the sections of Pn, Pw, Ps and Pe and takes the vertex Pc of the metal layer unit as the center:
wherein, VPc、VPE、VPW、VPNAnd VPSPc and P are respectivelyE、PW、PNAnd the potential at PS; wherein σPe、σPw、σPnAnd σPsRespectively, respectively
Electrical conductivities at Pe, Pw, Pn and Ps;
the above is to apply the finite volume method to the control volume centered at the vertex;
the electrical conductivity of a metal varies with the temperature of the material, and its relationship to temperature can be represented by the following equation:
where ρ is the resistivity, ρ0For resistivity at initial ambient temperature, Δ TMTemperature rise of discrete metal units compared to ambient temperature, αTAlpha for copper, for a material dependent temperature coefficientTIs positive, i.e. its resistivity increases with increasing temperature;
the conductivity of the midpoint of the cross section is considered to be right and left or upper and lower sides thereofThe conductivity corresponding to the equivalent resistance after the metal unit synthesis is the reciprocal of the sum of the two unit resistivities, sigmaPsConductivity σ Using C1 cell and C2 cellC1And σC2Represents:
therefore, the electrical conductivity of each unit of the metal layer is calculated according to the temperature distribution obtained in step S4; the conductivity corresponding to Pc can be approximately found by applying the conductivities corresponding to the four midpoints surrounding Pc:
for metal cells at line boundary locations, it may happen that the analyzed control volume cells are non-square, and the corresponding approximately discrete current continuity equation is:
wherein σPm1、σPm2、σPm3And σPm4Are each Pm1、Pm2、Pm3And Pm4The electrical conductivity of (d); vPvt1、VPvt2、VPvt3And VPvt4Potentials at Pvt1, Pvt2, Pvt2 and Pvt2, respectively;
in the PCB design process, 45 degree diagonal lines are designed, but they can be reconstructed by a square grid, and the discrete approximate current continuity equation of a small square control volume cell of quarter cell area with Ptr2 as the vertex at the line boundary is:
wherein σPmeAnd σPmnAre each PmeAnd PmnThe electrical conductivity of (d); vPtr2、VPc2And VPc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation shows that the current continuity equation is the same as that corresponding to the triangular control unit surrounding Ptr2 after filling the sawtooth gaps of the boundary, i.e. it proves that the partial sawtooth boundaries formed discretely by the square grid and the boundary of the 45-degree inclined line can be approximately equivalent in current continuity;
while the control volume units at the line corners are not directly equivalent to the actual line corners, the control volume units with the cross sections of Pmn3, Pmw3 and Pms3 surrounding the boundary vertex Pb3 have the length d of the cross section of the lower side Pms3 in the sawtooth boundarycA/2, while the length in the actual line is dcTherefore, such corners need to be individually processed in the algorithm to account for the approximate current continuity equation corresponding to their actual boundaries:
wherein σPmw3、σPmn3And σPms3Are each Pmw3、Pmn3And Pms3The electrical conductivity of (d); vPb3、VPb4、VPc4And VPc5Potentials at Pb3, Pb4, Pc4 and Pc5, respectively;
for the metal layer units connected by the via holes, the corresponding approximately discrete current continuity equation can also be obtained according to the above vertex-centered finite volume method, and assuming that vertices of adjacent metal layer units, of which the metal layer unit vertex Pc is connected by the via holes, are Pu and Pd, the corresponding approximately discrete current continuity equation becomes:
wherein, VPuAnd VPdAre each PuAnd PdThe potential of (d); sigmaPuAnd σPdElectrical conductivity at Pu and Pd, respectively;
finally, a summary matrix equation of approximate discrete current continuity equations corresponding to the upper and lower surface metal layers can be obtained:
wherein, V1And V2Respectively corresponding to the summary vectors of all the discrete unit top potential with unknown potential of the upper surface metal layer and the lower surface metal layer; mσ,1And Mσ,2Respectively as V in the upper surface matrix equation1With V in the lower surface matrix equation2The coefficient matrix of (a), is formed by conductivity parameters; kp,1And Kp,2Known terms in two matrix equations, such as corresponding terms describing known endpoint potentials or cross-sectional currents, are included respectively; mσV,1_2Is V in the upper surface matrix equation2The coefficient matrix of (A) is composed of conductivity parameters corresponding to the lower surface metal layer unit connected by the via hole and the upper surface metal unit discrete heat balance equation, for V2In cells not connected by vias, MσV,1_2The corresponding column in (1) is an all-zero column; mσV,2_1For V in the lower surface matrix equation1The coefficient matrix of (a) is formed by conductivity parameters corresponding to the upper surface metal layer unit connected by the via hole in the discrete thermal equilibrium equation of the lower surface metal unit, for V1In cells not connected by vias, MσV,2_1The corresponding column in (1) is an all-zero column; the above matrix equation can therefore also incorporate the approximate discrete current continuity equation for metal cells not connected to a via; v can be obtained by simultaneously solving the two matrix equations1And V2;
After the potential distribution of the top point in the line is obtained, the potentials corresponding to the center of each square cell and the center of the cross section around the cell are obtained, and the potential V at the center of the cell C2 is obtained by applying the finite volume method to the small diamond region around the center of the cell C2C2:
Wherein, VPvt4Potential at cell C2 apex Pvt 4; similarly, the potential V of the point Ps at the center of the cross section can be obtainedPs:
Wherein, VC1Potential at the center of cell C1;
thus, the lateral current density and the vertical current density corresponding to each cell center can be approximated, and the lateral current density J corresponding to the cell center can be approximated for the cell C1x,C1And longitudinal current density Jy,C1Can be expressed as:
then C1 corresponds to a total current density of:
for the triangular unit at the oblique line boundary, assuming that the potential is linearly distributed in the region, the current density of the corresponding triangular region is obtained, and the linear potential distribution corresponding to the triangle formed by three points Ptr2, Ptr1 and Pc1 can be expressed as follows:
wherein, the third-order square matrix in the equation is composed of operation formulas corresponding to the coordinates of three vertexes under a Cartesian coordinate system, SΔIs a triangular area, so that the coefficient a corresponding to the potential linear distribution expression can be obtained0、a1And a2Then, the transverse current density J within the triangle is knownx,ΔAnd longitudinal current density Jy,ΔComprises the following steps:
wherein σx,ΔConductivity corresponding to the triangular region;
on the premise of obtaining the current density corresponding to the metal cell, the joule heating rate per unit volume can be obtained, for example, for cell C1, the joule heating rate per unit volume q can be obtainedC1Can be expressed as:
qC1=EC1·JC1=(Jx,C1 2+Jy,C1 2)/σC1=JC1 2/σC1
wherein E isC1The electric field strength at the center of cell C1;
the joule heating rate is related to the electrical conductivity, and the electrical conductivity is affected by the temperature, so when the joule heating of the line is counted, iterative operation between the temperature distribution and the joule heating rate distribution is required to be carried out, when the error of the iterative operation is smaller than a set error, the temperature distribution and the joule heating rate distribution can be approximately calculated, and the error can be set as the maximum value of the allowable unit temperature change in two iterative calculations. q. q.sC1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the upper and lower surface metal layersMJ,uAnd q isMJ,d。
The steady-state thermal analysis method for the double-sided PCB structure provided by the invention has the following beneficial effects:
the method for coupling numerical value dispersion and temperature analysis solution can avoid integral dispersion and solution of the integral PCB structure by solving the temperature distribution temperature analysis solution of the designated areas on the upper surface and the lower surface of the substrate insulating layer, can realize approximate representation of the thermal diffusion between the metal circuit layers through the metal through holes and the metal layers on the basis of a finite volume method, can approximately correct the homogeneous structure hypothesis of the temperature analysis solution on the substrate insulating layer containing the through holes, can simultaneously count the boundary condition that the surface heat flow density distribution of the substrate insulating layer is subjected to the thermal diffusion effect of the metal layers, and further realizes the mutual coupling solution of the temperature analysis solution matrix equation and the metal layer thermal diffusion matrix equation to obtain the surface temperature distribution of the PCB.
Drawings
Fig. 1 is a schematic diagram of structured discrete metal layer unit heat conduction.
FIG. 2 is a cross-sectional view of a via.
Fig. 3 is a schematic diagram of a discrete unit of the circuit 1.
Fig. 4 is a schematic diagram of a discrete element of the circuit 2.
FIG. 5 is a schematic diagram of a zigzag boundary equivalent to a 45 degree diagonal boundary.
Fig. 6 is a schematic diagram of a corner case where a 45 degree boundary may occur.
FIG. 7 shows a general structure and flow of a corresponding procedure of a dual-sided PCB structure coupled thermal analysis modeling method.
Fig. 8 shows the layout and temperature distribution of the power MOSFET pad leads.
Fig. 9 is a circuit layout of upper and lower surface pad pins and upper surface temperature distribution after the addition of thermally conductive vias.
Fig. 10 is a diagram of a double-sided wiring structure.
Fig. 11 is a top surface line joule heating rate distribution.
Fig. 12 is a lower surface line joule heating rate distribution.
Fig. 13 is a top surface line temperature distribution.
Fig. 14 is a bottom surface line temperature distribution.
Detailed Description
The following description of the embodiments of the present invention is provided to facilitate the understanding of the present invention by those skilled in the art, but it should be understood that the present invention is not limited to the scope of the embodiments, and it will be apparent to those skilled in the art that various changes may be made without departing from the spirit and scope of the invention as defined and defined in the appended claims, and all matters produced by the invention using the inventive concept are protected.
According to one embodiment of the application, the steady-state thermal analysis method for the double-sided PCB structure comprises the following steps:
s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation;
s2, calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the boundary condition of the substrate insulating layer, and calculating the heating rate of the component as the known heat source condition into a matrix equation expression of an analytic solution;
s3, approximating a heat balance equation of a discrete metal layer based on a finite volume method, summarizing to obtain a metal layer heat diffusion lumped matrix equation, and counting the heat conduction between the metal layer and a surface heating device thereof and the heat conduction of the metal via hole;
s4, constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the corrected insulating layer temperature analytic solution matrix equation and the metal layer thermal diffusion matrix equation based on numerical value dispersion;
s5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and calculating potential distribution in the metal layer, so as to further calculate current density distribution and Joule heating distribution of the metal layer, and calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with temperature distribution to calculate temperature distribution of the surface where each layer of the PCB is positioned under the condition of Joule heating of a circuit;
and S6, verifying the calculation accuracy of the thermal analysis method of the double-sided PCB structure by adopting a COMSOL modeling operation result.
The above steps will be described in detail below according to one embodiment of the present application.
Firstly, thermal analysis premise hypothesis is required to be carried out on the structural characteristics of the double-sided PCB:
the double-sided PCB structure is mainly composed of a core substrate insulating layer and metal circuit layers covered on the upper and lower surfaces of the core substrate insulating layer. A common substrate layer material is FR-4 glass epoxy, which has a thermal conductivity of about 0.3W/mK and a crystallization temperature of about 120 ℃, and when this temperature is exceeded, both the mechanical and insulating properties of FR-4 will be significantly reduced. The metal layer is mainly made of copper, and the thermal conductivity is about 390W/mK. While the thickness of the substrate insulating layer (>0.5mm) is much larger than the thickness of the metal layer (generally less than 0.1 mm). The PCB surface is covered with a soldemask solder resist coating and a silk screen printing layer for marking devices and printing trademarks and other information. Because the area of the silk-screen layer on the surface of the PCB is low, the influence of the silk-screen layer on the heat transfer is ignored. In addition, because the solder resist coating is thin (typically 1mil) and has low thermal conductivity (0.25W/mK), its equivalent thermal resistance is taken into account in the heat transfer coefficient between the PCB surface and the external environment. The Heat Transfer Coefficient (HTC) is also one of the important thermal boundary conditions for deriving a temperature-resolved solution.
However, the existence of the components will destroy the assumed condition that the upper and lower surfaces of the PCB and the external environment exchange heat under a certain HTC, and the PCB and the external environment also exchange heat through the surfaces of the components, so that heat flow transmission also exists between the components and the PCB, and the package of the components will also affect the assumed thermal boundary condition.
Although the insulating material occupies most of the PCB structure, the large difference in thermal conductivity between the metal layer and the insulating layer results in the metal layer wiring having much higher thermal diffusivity than the insulating material even in the micrometer-sized thickness, which cannot be ignored in the thermal analysis. However, since the metal layer is thin and has high thermal conductivity, the temperature gradient between the upper and lower surfaces thereof can be ignored, and the insulating layer is closely attached to the metal layer, so that it can be assumed that the temperature distribution of the metal layer is the same as that of the surface region of the insulating layer in contact with the metal layer. The substrate insulating layer is a main constituent material of the PCB, but since it is not a heat source, the temperature distribution inside the substrate insulating layer is not a main solution target.
Therefore, in the modeling method, the temperature distribution analytic solution of the surface of the insulating layer is used for providing the temperature distribution expression of the metal layer and the heating device which are in contact with the insulating layer, namely the temperature distribution expression is used as one of the coupling variables of the numerical discrete analysis and the analytic solution; meanwhile, a numerical value dispersion method based on finite volume is used for counting the heat diffusion inside the metal layer and the heat conduction process between the metal layers through the metal through holes, and an approximate steady-state heat diffusion matrix equation of the metal layer is given, wherein the heat diffusion matrix equation contains the heat flow density distribution actually transmitted by the metal layer to the insulating layer, and the distribution is also used as one of coupling variables; and finally solving the coupling variable through a simultaneous analytic solution matrix equation and a thermal diffusion matrix equation. And finally, obtaining the core temperature of the heating element by using the thermal resistance parameters of the heating element.
In addition, because the heat dissipation areas of the four sides of the PCB are generally much smaller than the heat dissipation areas of the upper and lower surfaces of the PCB, and the main heat sources are both on the upper and lower surfaces, the heat exchange between the four sides of the PCB and the outside is neglected. The method is therefore not suitable for thermal analysis of PCBs of the type in which heat transfer is achieved primarily through the side and where the side has a plug section, etc. in an application.
The following detailed description of the algorithm is made under the above-described assumption:
step S1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer, and obtaining an insulating layer temperature analytic solution matrix equation, wherein the equation specifically comprises the following steps:
the corresponding steady state heat balance Laplace equation of the insulating layer in the double-sided PCB and the boundary conditions of the periphery, the upper surface and the lower surface are as follows:
wherein T is the temperature variation compared with the ambient temperature or the temperature when the ambient temperature is 0 ℃, and the density distribution function of the heat transmitted from the upper surface and the lower surface is qiu(x, y) and qid(x,y),hu(W/m2K) And hd(W/m2K)Heat Transfer Coefficients (HTC), k, of the upper and lower surfaces of the insulating layer, respectivelyi(W/mK) is insulation layer thermal conductivity, Lx、Ly、 LinThe length of the insulating layer in the x and y directions and the thickness of the insulating layer in the z direction in a cartesian coordinate system, respectively.
The expression form of Fourier series analytic solution corresponding to the steady-state heat balance equation is as follows:
wherein, C1,C2,Cn,m,Cγn,mIs a Fourier series coefficient, betan、μm、γn,mCharacteristic values in a Fourier series;
according to the principle of thermal effect superposition, the thermal effect generated by the upper surface device and the lower surface device or other heat sources in the insulating layer can be regarded as the result of superposition of the thermal effects after the upper surface heat source and the lower surface heat source act independently, and at the moment, the equation (1) is changed into:
therefore, only the analytic solution corresponding to the surface temperature distribution of the insulating layer with the single-sided structure needs to be deduced, and the analytic solution corresponding to the double-sided structure can be obtained through superposition. Only the solving process of the thermal effect corresponding to the analytic solution theta generated by the upper surface heat source will be described, and the solving process of eta is similar to the solving process. Substituting the order solution 2 into the lower surface (z ═ L)in) Corresponding thermal boundary conditions can be obtained, when h is reached, the corresponding relation between partial coefficient and characteristic value in solution can be obtaineddWhen not zero, C2And Cγn,mThe corresponding expression is:
when h isdIs zero, C1,Cγn,mCan be expressed as:
discussion of h onlydGeneral derivation of other coefficients in an example analytic solution corresponding to a single-sided structure when non-zero, hdThe derivation process for zero time is similar and will not be described here. The number of stages is substituted into the thermal boundary condition of the upper surface (z is 0), and the following results are obtained:
multiplying both sides of the above equation by cos (. beta.)nx)cos(μmy) and integration is performed in the surface heating area, the integration after multiplication of different numbers of stages is zero according to the characteristics of trigonometric function, so that the left side of the equation only contains cos2(βnx)cos2(μmy) is not zero after integration, and the right side of the equation is used for each heating unit q in the heating areaiu,iIntegrating and then summing, and calculating to obtain coefficient C1And Cn,m:
Wherein d isqIntegration Domain Sq for grid cell Length after discretization of surface Heat Source gridiu,iI.e. representing an incoming heat flow density of qiu,i(W/m2) Surface discrete unit area of (d)nAnd deltamTo resolve dimensionless intermediate parameters of the derivation process, qiu,1~qiu,N1In order to discretize the unit heat flux density of the heat source grids on the surface of the insulating layer, if the heat flux density transmitted by a certain heat source is uniform, grid discretization is not carried out on the longitudinal heat flux density transmitted by the area where the heat source is located; c1The molecule in the expression is the heat quantity transmitted in the whole unit time, Cn,mThe double integral summation term in the expression can be directlyIntegration is performed in the area where the heat source is located.
Further, the summation operation of Fourier series can be expressed by means of vector operation, and for the above analytic solution, qiu,1~qiu,N1Namely, a column vector q of heat flux density is constructediuThen, 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 Ru(x,y,z)Expanding into a thermal resistance matrix R according to the coordinates of the corresponding distribution points of the regionuWherein each row is corresponding to a thermal resistance row vector for solving the temperature at a certain point, and then R is addeduMultiplying the temperature distribution vector by the corresponding heat flow density vector q to obtain a temperature distribution vector of the specified area, namely the heat effect generated by the corresponding q in the specified area; and finally, superposing Fourier series analytic solutions in a matrix form corresponding to theta and eta of the specified region to obtain the temperature distribution of the region.
Thus, the matrix equation corresponding to the analytic solution temperature distribution vector of the upper and lower surfaces of the insulating layer in the double-sided PCB will have the following form:
wherein, the distribution vectors of the heat flux density transmitted into the insulating layer by the upper and lower surface metal layers are qM,u、qM,dWherein q isuAnd q isdThe distribution vector q of heat flux density transferred between the upper and lower surface heating components and the covering region thereofuComprising qiu(x, y) heat flux density distribution after dispersion, qdComprising qid(x, y) a heat flow density distribution after discretization; rMu,uIs qM,uAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMu,dIs qM,dAt TuThe defined area generating a thermal effectA corresponding thermal resistivity matrix; rMd,uIs qM,uAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMd,dqM,dAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; riu,uIs quAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; riu,dIs qdAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rid,uIs quAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rid,dIs qdAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; due to quAnd q isdThe 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 equationiu(x, y) and qid(x, y) heat flow transferred between the associated top and bottom surface heat generating components and the surface area of the insulating layer directly covered thereby, Riu,u、Riu,d、Rid,uAnd Rid,dWill include a partial all-zero column vector;
step S2, calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the thermal boundary condition of the substrate insulating layer, and obtaining a corrected insulating layer temperature analytic solution matrix equation expression, which specifically comprises the following steps:
the thermal boundary conditions according to which the analytical solution of the surface temperature distribution of the insulating layer of the substrate is solved are mainly the heat flux density distribution transmitted by the upper surface and the lower surface of the insulating layer of the substrate and the heat transfer coefficient for representing the heat exchange between the surface and the external environment. However, the thermal boundary conditions ignore the fact that the partial region covered by the heat generating component is only in heat exchange with the component, rather than directly with the external environment, and do not account for heat exchange between the surface of the component and the external environment. To correct the assumed thermal boundary conditions, the area occupied by the components must be based on huThe multi-calculated heat exchange heat flow density is compensated, and on the other hand, the thermal resistance parameters provided by the components can be used for deducing the position between the related components and the PCB surfaceTo heat transfer.
Setting the heat flux density vector of the covering area of the upper surface heating device as qh,uIt can be expressed as follows:
qh,u=qu+huTJ,u (10)
wherein, TJ,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 environmentuThe 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 huCompensating for the multi-calculated heat exchange; in addition, the thermal resistance parameters provided by the components can be used for deducing the heat transfer between the related components and the PCB surface, and the specific deduction process is as follows:
the two formulas are combined to obtain:
qh,u=CP,uqJ,u-CT,uTJ,u
wherein the temperature distribution of the upper surface element device surface is TtopThe heat transfer coefficient of the surface to the outside is htopThe equivalent heat conduction corresponding to the heat sink existing on the surface thereof is also taken into accounttop(ii) a Mean temperature of the core is Tj(ii) a The average thermal resistance between the surface and the core is psijt(ii) a The average thermal resistance between the junction center and the bottom of the device is psijc(ii) a Equivalent thermal resistance of solder mask coating is psic,u(ii) a The heat flux density transferred from the core to the surface of the component is qjt(ii) a The heating volume density vector corresponding to the upper surface element device is PJ,u(ii) a The diagonal matrix composed of height parameters of the components is CJ,u(ii) a The density vector of the heating surface corresponding to the upper surface element device is qJ,u;CP,uFor taking into account the thermal resistance parameter of the deviceJ,uThe coefficient of (a); cT,uFor taking into account the corresponding T after the thermal resistance parameter of the deviceJ,uThe coefficient of (a);
with qh,uAlternative upper surface temperature distribution analytical solution TuQ in (1)u:
Wherein, TJ,dThe temperature distribution vector of the coverage area of the lower surface heating device is obtained; riu,uCP,uqJ,uThe term is a vector of known constants, and TuIncluding TJ,uAll elements of the vector can be transformed by adding zero to the matrix coefficients in a rising-dimension manner to realize the use of T in the matrix equationuInstead of TJ,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:
CJu,uTu=Riu,uCT,uTJ,u (14)
wherein, CJu,uTo realize TJ,uAnd TuAnd (3) equivalently transforming the coefficients of the related terms to further obtain:
Tu=RMu,uqM,u+RMd,uqM,d+Rid,uqd+Riu,uCP,uqJ,u-CJu,uTu (15)
similarly, when there is a device on the lower surface, the heat flow transferred between the device on the lower surface and the surface unit, which takes into account the analytical solution, is also corrected, i.e. q in the above formuladItem replacement:
in the same way, the T corrected by the thermal boundary condition can be obtaineddCorresponding analytic solution matrix equation, and after correction, TuAnd TdThe corresponding simultaneous analytical solution matrix equation becomes:
wherein, with CJu,uAnd CJd,uFunction of (C) is similar to the definitionJu,d、CJd,dRespectively for realizing T in the lower surface analytical solution matrix equationJ,uAnd TuRelated item and TJ,dAnd TdCoefficients of the correlation term equivalent transform;
s3, approximating a discrete metal layer heat balance equation based on a finite volume method, summarizing to obtain a metal layer heat diffusion matrix equation, and calculating the heat conduction between the metal layer and a surface heating device thereof and the heat conduction of the metal via hole therein, wherein the method specifically comprises the following steps:
s3.1, approximating a heat balance equation of a discrete metal layer based on a finite volume method;
the steady state thermal equilibrium equation for the metal layer has the form of poisson's equation:
wherein k ism(W/mK) represents the metal layer thermal conductivity; q. q.sc(x,y)(W/m3) Is a joule heating rate distribution function per unit volume of the metal layer.
To calculate the thermal diffusion in the metal layer, the above thermal equilibrium equation will be discretized mainly using a finite volume method, and the metal layer will be discretized using a structured square grid, as shown in fig. 1. Wherein the contribution of the metal layer to the heat diffusion is realized by the heat conduction among the metal layer units, and the heat flow density conducted between the central unit and the four surrounding units is q as shown in fig. 1n,qs,qwAnd q ise(W/m2) And the heat exchange between the top of the device and the component is expressed as heat flow density qu(W/m2) The heat conduction with the bottom insulating layer is formed by qd(i.e., constituting q in the aforementioned analytical solution equation)M,uElements of a vector) (W/m)2) Represents, but onlyJoule heat in the interior of the element is generated by the volume heat generation rate qc(W/m3) And (4) representing. When Joule heating of certain lines or elements is not considered in the thermal analysis, the corresponding qcI.e. zero.
Performing triple volume integration on the above equation in a unit, and applying the divergence theorem, we can obtain:
wherein q iscIs the unit volumetric heat rate, dcIs the unit side length of the structured and discrete metal unit, dmMultiplying both ends of the equation by k for the thickness of the metal layermAnd 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 unitcAnd its east cell temperature TEThe corresponding Taylor-series expansion is:
wherein, the position (x)e,ye) The midpoint of the connecting line of the corresponding central unit and the right unit; o (d)c 3) Is the sum of all terms containing partial derivatives of the third order and above, and d contained thereincThe power exponent of the power term is greater than or equal to 3, since d is generallycSmaller value of (a) can be represented by O (d)c 3) Are omitted. And then subtracting the two expressions to approximately obtain the partial derivative of the temperature on the corresponding section:
similarly, other partial derivatives in equation 20 can be processed approximately linearly, and finally, the heat dissipation balance equation can be realized, and an approximate linear relationship between the temperatures of adjacent units is obtained:
therefore, the temperature relationship between adjacent units is established by the equation after dispersion, so that the subsequent dispersion heat balance equation corresponding to all the metal layer units can be conveniently induced into a lumped matrix equation representing the metal layer heat diffusion, and when the units are connected with other adjacent metal layer units through metal via holes, the q-axis heat resistance can be based on the metal via hole heat resistanceuOr qdAnd (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 Lin(m) via radius rv(m) via thermal conductivity kv(W/mK),dpw(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 isin(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 rangecu: kFR-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 rv<2000dpw(dpwTypically 1-4 mil), the thermal resistance of a thin-walled metal via will be lower than the thermal resistance of the same via volume filled with FR-4, and this condition is met in most PCB via designs; on the other hand, when the via hole is a solid via hole filled with copper completely, the via hole is generally used as a heat conducting via hole due to low thermal resistance; therefore, the thermal conduction contribution of the vias must be accounted for in the thermal analysis of the PCB structure.
Therefore, on the premise of taking the heat conduction of the metal via hole into account, the heat balance equation of the discrete metal layer on the upper surface can be obtained:
wherein k ismThe thermal conductivity of the metal layer; q. q.sc,uThe joule heating rate distribution per unit volume of the upper surface metal layer; t isc,uIs the upper surface metal cell temperature; dcThe unit side length of the structured and discrete metal unit is shown; dmIs the thickness of the metal layer; q. q.su,uIs Tc,uThe longitudinal heat flow density of the upper surface of the corresponding unit can be the density of heat flow transferred between the corresponding unit and the heating device; q. q.sd,uIs Tc,uThe longitudinal heat flux density transmitted between the lower surface of the corresponding unit and the insulating layer; t isE,uIs Tc,uThe temperature of the adjacent metal unit at the east side of the corresponding unit; t isW,uIs Tc,uThe temperature of the corresponding unit west side adjacent metal unit; t isN,uIs Tc,uThe temperature of the adjacent metal unit at the north side of the corresponding unit; t isS,uIs Tc,uThe temperature of adjacent metal units on the south side of the corresponding unit; psiNv,zIs the unit thermal resistance of the metal via; t isc,dvThe temperature of the lower surface metal layer unit connected through the via hole is shown;
similarly, a discrete heat balance equation corresponding to the lower surface metal layer unit can also be obtained:
wherein q isc,dThe joule heating rate distribution per unit volume of the lower surface metal layer; t isc,dIs the lower surface metal unit temperature; dcThe unit side length of the structured and discrete metal unit is shown; dmIs the thickness of the metal layer; q. q.sd,dIs Tc,dThe longitudinal heat flux density of the lower surface of the corresponding unit can be the density of heat flux transferred between the corresponding unit and the heating device; q. q.su,dIs Tc,uThe longitudinal heat flux density transmitted between the upper surface of the corresponding unit and the insulating layer; t isE,dIs Tc,dThe temperature of the adjacent metal unit at the east side of the corresponding unit; t isW,dIs Tc,dCorresponding to the temperature of the adjacent metal unit on the west side of the unit; t isN,dIs Tc,dThe temperature of the adjacent metal unit at the north side of the corresponding unit; t isS,dIs Tc,dThe temperature of the adjacent metal unit on the south side of the corresponding unit; t isc,uvThe temperature of the lower surface metal layer unit connected through the via hole is shown;
step S3.3, the upper and lower surface metal layer thermal diffusion matrix equation based on numerical value dispersion is summarized:
wherein q isMJ,uAnd q isMJ,dRespectively corresponding unit joule heating vectors of the upper surface metal layer and the lower surface metal layer; h istopHeat of the surface of the component and the outsideA transfer coefficient; mMV,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; mMV,dThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; mV,dThe coefficient of the discrete heat balance equation of the upper surface metal layer unit corresponding to the temperature of the lower surface metal layer unit is formed; mV,uThe coefficient corresponding to the temperature of the upper surface metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; due to qh,u、 qh,dThe compensated heat flow directly transmitted through the surface of the insulating layer is counted in an analytic solution matrix equation, and only the heat flow transmitted by the heating device and the part of the covering metal layer of the heating device is counted in a lumped matrix equation of heat diffusion, so that CM,u、CM,dWill also include partial all-zero column vectors;
q is to beh,uAnd q ish,dSubstituting the expression of (a) into the thermal diffusion matrix equation to obtain:
wherein, CMu,uAnd CMd,dRespectively for realizing T in the above equationJ,uAnd TuRelated item and TJ,dAnd TdCoefficients of the correlation term equivalent transform;
step S4, according to the corrected insulation layer temperature analytic solution matrix equation and the metal layer thermal diffusion matrix equation based on numerical value dispersion, constructing a thermal analysis coupling equation set of the double-sided PCB structure, which specifically comprises the following steps:
from the temperature-resolved solution matrix equation 17, and the thermal-analysis coupling equation set corresponding to the numerical-dispersion-based thermal-diffusion matrix equation 30, the double-sided PCB structure can be expressed as:
wherein, Tu,Td,qM,u,qM,dIf the temperature vector is unknown, and other vector and matrix coefficients are known, the unknown temperature vector and the unknown heat flow density vector can be obtained by combining the four matrix equations.
S5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and obtaining potential distribution in the metal layer, so that current density distribution and Joule heating distribution of the metal layer can be further obtained, the Joule heating distribution can be included in a metal layer thermal diffusion matrix equation, iterative calculation is carried out between the Joule heating distribution and temperature distribution, and temperature distribution of the surface where each layer of the PCB is included under the premise of Joule heating of a circuit is obtained, and the method specifically comprises the following steps:
s5.1, approximating a discrete metal layer current continuity equation based on a vertex-centered finite volume method, which specifically comprises:
the current continuity equation in the metal layer can be expressed as:
wherein, J (A/m)2) 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 squarePe、JPwThe transverse current densities at Pe and Pw, respectively; j. the design is a squarePnAnd JPsThe 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 PEThe potential of (c) is developed using a taylor series:
wherein, VPeIs PeThe potential of (d); (x)Pe,yPe) Is PeCoordinates of (c); dcIs the unit length of the metal layer unit; further subtracting the two to obtain a partial differential expression based on the second-order taylor series approximation:
therefore, an approximate discrete equation corresponding to a current continuity equation of a square region surrounded by the cross section where Pn, Pw, Ps and Pe are located can be obtained by using the metal layer unit vertex Pc shown in fig. 3 as the center (not connected with the metal via):
wherein, VPc、VPE、VPW、VPNAnd VPSPc and P are respectivelyE、PW、PNAnd the potential at PS;
the above is to apply the finite volume method to the control volume centered at the vertex;
the electrical conductivity of a metal varies with the temperature of the material, and its relationship to temperature can be represented by the following equation:
where ρ is the resistivity, ρ0For resistivity at initial ambient temperature, Δ TMTemperature rise of discrete metal units compared to ambient temperature, αTAlpha for copper, for a material dependent temperature coefficientTIs positive, i.e. its resistivity increases with increasing temperature;
the conductivity of the midpoint of the cross section is considered to be the conductivity corresponding to the equivalent resistance of the synthesized metal units 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. sigmaPsThe conductivities of the C1 cell and the C2 cell can be usedC1And σC2Represents:
the electrical conductivity of each cell of the metal layer is therefore calculated from the temperature distribution obtained in claim 5; the conductivity corresponding to Pc can be approximately found by applying the conductivities corresponding to the four midpoints surrounding Pc:
for metal cells at line boundary positions, it may happen that the control volume cell to be analyzed is non-square, such as the control cell No. 2 (which is the smaller square control cell No. 1 below) surrounded by dark black multi-point lines shown in fig. 4, and its corresponding approximately discrete current continuity equation is:
wherein σPm1、σPm2、σPm3And σPm4Are each Pm1、Pm2、Pm3And Pm4The electrical conductivity of (d); vPvt1、VPvt2、VPvt3And VPvt4Potentials at Pvt1, Pvt2, Pvt2 and Pvt2, respectively;
in the PCB design process, 45 degree inclined lines are often designed, but the square grid can be better reconstructed, for example, for a small square control volume unit with a quarter cell area with Ptr2 as the vertex at the boundary of the upper left corner line of fig. 5, the discrete approximate current continuity equation is:
wherein σPmeAnd σPmnAre each PmeAnd PmnThe electrical conductivity of (d); vPtr2、VPc2And VPc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation indicates that the current continuity equation is the same as the current continuity equation corresponding to the triangle control cell surrounding Ptr2 after filling the boundary sawtooth gap. That is, it is proved that the partial zigzag boundary formed discretely in the square grid and the boundary of the 45-degree inclined line are approximately equivalent in current continuity.
The control volume units at the line corners are not directly equivalent to the actual line corners, for example, the control volume units with the cross sections of Pmn3, Pmw3 and Pms3 around the boundary vertex Pb3 have the length d of the cross section of the lower side Pms3 in the sawtooth boundarycA/2, while the length in the actual line is dcTherefore, 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 Pmw3、Pmn3And Pms3The electrical conductivity of (d); vPb3、VPb4、VPc4And VPc5Potentials 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:
wherein, VPuAnd VPdAre each PuAnd PdThe potential of (d); sigmaPuAnd σPdElectrical conductivity at Pu and Pd, respectively;
s5.2, summarizing to obtain a lumped matrix equation for describing the linear relation among the electric potentials of the metal layer units, and solving the electric potential distribution of the top points of the metal layer units, wherein the lumped matrix equation specifically comprises the following steps:
summarizing and summarizing the discrete current continuity equations corresponding to the metal layer units of various types to obtain a matrix equation of an approximate discrete current continuity equation corresponding to each metal layer, and finally obtaining a summarized matrix equation of the approximate discrete current continuity equations corresponding to the upper and lower surface metal layers:
wherein, V1And V2Respectively corresponding to the summary vectors of all the discrete unit top potential with unknown potential of the upper surface metal layer and the lower surface metal layer; mσ,1And Mσ,2Respectively V in the upper surface matrix equation1With the lower surface matrix equationV2The coefficient matrix of (a), is formed by conductivity parameters; kp,1And Kp,2Known terms in two matrix equations, such as corresponding terms describing known endpoint potentials or cross-sectional currents, are included respectively; mσV,1_2Is V in the upper surface matrix equation2The coefficient matrix of (a) is formed by conductivity parameters corresponding to the lower surface metal layer unit connected by the via hole and the upper surface metal unit discrete heat balance equation, for V2In cells not connected by vias, MσV,1_2The corresponding column in (1) is an all-zero column; mσV,2_1For V in the lower surface matrix equation1The coefficient matrix of (a) is formed by the conductivity parameters corresponding to the upper surface metal layer unit connected by the via hole in the discrete thermal equilibrium equation of the lower surface metal unit, for V1In cells not connected by vias, MσV,2_1The corresponding column in (1) is an all-zero column; the above matrix equation can therefore also incorporate the approximate discrete current continuity equation for metal cells not connected to a via; v can be obtained by simultaneously solving the two matrix equations1And V2;
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 C2C2:
Wherein, VPvt4Potential at cell C2 apex Pvt 4; similarly, the potential V of the point Ps at the center of the cross section can be obtainedPs:
Wherein, VC1Potential at the center of cell C1;
step S5.4, further solving the current density distribution and the Joule heating distribution of the metal layer, and calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with the temperature distribution to obtain the temperature distribution of the surface where each layer of the PCB is located under the premise of calculating the Joule heating of the circuit, which specifically comprises the following steps:
on the basis of the potential distribution of the metal layer, the transverse current density and the longitudinal current density corresponding to the center of each cell can be approximately obtained, for example, for the cell C1, the transverse current density J corresponding to the center of the cellx,C1And longitudinal current density Jy,C1Can be expressed as:
then C1 corresponds to a total current density of:
for the triangular cells at the oblique line boundary, the electric potential is assumed to be linearly distributed in the region, so as to obtain the current density of the corresponding triangular region, for example, the linear electric potential distribution corresponding to the triangle enclosed by three points Ptr2, Ptr1 and Pc1 can be expressed as follows:
wherein, the third-order square matrix in the equation is composed of operation formulas corresponding to the coordinates of three vertexes under a Cartesian coordinate system, SΔIs a triangular area, so that the coefficient a corresponding to the potential linear distribution expression can be obtained0、a1And a2Then, the transverse current density J within the triangle is knownx,ΔAnd longitudinal current density Jy,ΔComprises the following steps:
wherein σx,ΔConductivity corresponding to the triangular region;
on the premise of obtaining the current density corresponding to the metal cell, the joule heating rate per unit volume can be obtained, for example, for cell C1, the joule heating rate per unit volume q can be obtainedC1Can be expressed as:
qC1=EC1·JC1=(Jx,C1 2+Jy,C1 2)/σC1=JC1 2/σC1 (52)
wherein E isC1The electric field strength at the center of cell C1; q. q.sC1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the upper and lower surface metal layersMJ,uAnd q isMJ,d。
The joule heating rate is related to the electrical conductivity, and the electrical conductivity is affected by the temperature, so when the joule heating of the line is counted, the iterative operation between the temperature distribution and the joule heating rate distribution needs to be performed, when the error of the iterative operation is smaller than the set error, the temperature distribution and the joule heating rate distribution can be approximately calculated, and the error can be set as the maximum value of the unit temperature change allowed in the two iterative calculations, as shown in fig. 7.
The general program structure for realizing the thermal analysis modeling of the double-sided PCB structure is as follows:
the general structure and flow of the PCB structure thermal analysis program which can be written by using the method are shown in FIG. 7, and the program can be divided into an initialization program and a main program. If an analysis model is newly established in an initialization program, a circuit diagram is firstly imported, all components such as a circuit, a device and a bonding pad need to be identified and position confirmed in the process of analyzing the circuit by the program, and then a coefficient matrix required by summarizing a coupling equation set can be calculated, wherein the coefficient matrix comprises coefficient matrixes required by analyzing an antipyretic resistance coefficient matrix and a numerical value discrete equation and taken into account of heat conduction of a metal layer and thermal resistance parameters of the device and a radiator. If the existing model is modified, the step of setting/modifying system parameters can be directly carried out, and finally, the data is stored.
The main program comprises two parts of direct coupling calculation temperature distribution and iterative calculation temperature distribution, and the direct coupling calculation calculates a thermal analysis result according to a coupling equation set by calling each coefficient matrix and a heating density vector of a known device. If the joule heating of part of the metal layer circuit needs to be taken into account, the influence of the temperature on the metal conductivity needs to be taken into account, so that iterative operation between the joule heating and the temperature distribution is needed to obtain the steady-state circuit temperature distribution. In the iterative calculation part, the Joule heating density distribution of the circuit is calculated according to the assumed temperature distribution of the circuit, then the coupling solution of the corresponding temperature distribution is calculated, the new Joule heating density distribution of the circuit is calculated according to the temperature distribution, then the new coupling solution is calculated, and the iterative calculation is carried out in sequence until the maximum difference value of the unit temperatures in the two temperature distributions is lower than the given error limit value.
According to one embodiment of the present application, a specific PCB structure is analyzed to verify the computational accuracy of the present algorithm:
and the MATLAB is utilized to complete the double-sided PCB structure thermal analysis modeling program capable of analyzing the two-dimensional line structure diagram based on the method. In the program, the thermal parameters such as the size of the PCB, the thickness of the metal layer and the insulating layer, the ambient temperature, the thermal conductivity and the surface heat transfer coefficient of the material, the thermal resistance of surface coatings such as a solder mask, the thermal resistance of the device and the like can be set, and the functions of setting the electrical parameters such as the power consumption of the device, the current and the voltage of the double-end line, the branch current of the multi-branch line, the terminal potential of. The previous studies of the inventor have concluded that a sufficient requirement for the known electrical parameters of the line, which are required for the potential distribution and current density analysis of the multi-branch line, is that the sum of the known branch currents or the number of branch terminal potentials needs to be equal to the number of line terminals, and at least one terminal potential is known.
a. Power MOSFET thermal analysis modeling example
An example of thermal analysis of a local MOS pin pad structure on a PCB and further optimization of the structure is first given. A power MOS with the model number of STB20N95K5 is used as a main circuit switch on a double-sided PCB of a certain vehicle-mounted DC-DC converter, and specific parameters of a pin pad structure of the MOS are shown in a table below.
TABLE 1 MOSFET pad Pin line parameters
Size of substrate | Thickness of copper foil | Via hole wall thickness | Pore diameter | Heat generation rate | Bottom cold plate temperature | Type of package |
22x25x2(mm) | 35μm | 25.4μm(1mil) | 1mm | 3.4 |
60℃ | D2PAK |
The original structure and the thermal analysis temperature distribution result of the MOS pin pad are shown in fig. 8, which is the most important heating device in the converter PCB, and here, the contribution of other parts of the PCB to the MOS heat conduction is ignored in the analysis process, and only the heat conduction from the PCB to the bottom cold plate is taken into account, so as to simulate the theoretical temperature rise extreme value of the MOS under the working condition. And simultaneously modeling the structure in COMSOL and carrying out thermal analysis so as to check the reasonability and the accuracy of the modeling method through comparison. As can be seen from the analysis results, the thermal distribution of the analysis results obtained based on the thermal analysis method is basically consistent with the COMSOL analysis results, and the temperature difference of the two hot spots accounts for less than 1% of the total hot spot temperature rise. However, the hot spot temperature of the lead bonding pad region exceeds 195 ℃, and the MOS cannot work normally, so that the lead bonding pad structure needs to be improved.
Because the bottom is a constant temperature cold plate, in order to enhance the heat conduction from the MOS to the bottom, some metal through holes are added in the pin pad area, and a layer of copper-clad metal layer is also added at the bottom, the structure and the thermal analysis temperature distribution result are shown in fig. 9, and the related analysis parameters and the structure discrete parameters corresponding to the coupled thermal analysis program and the COMSOL model are given at the same time.
TABLE 2 MOSFET pad pin line simulation calculation parameters
Therefore, under the condition of adding the metal heat conduction through hole, the temperature rise of the hot spot of the welding pad area is reduced by over 90 ℃, the core temperature of the welding pad area does not exceed 106 ℃ according to the calculation of the related thermal resistance parameters of the MOS package, and the welding pad area is in a safe working range. Compared with simulation results of COMSOL, the simulation results show that on one hand, the operation principle of the coupling thermal analysis is reasonable and accurate, and on the other hand, on the premise of having similar precision, the total number of numerical discrete units of the PCB in the coupling thermal analysis related by the invention is only 7495 of the pad area, which is far lower than the number of units of the PCB which is wholly discrete by COMSOL based on a finite element analysis method, so that the method has obvious advantages in numerical analysis efficiency.
b. Modeling example for joule heating thermal analysis of double-sided PCB (printed Circuit Board) circuit
In the following, the joule heating and the temperature distribution of a double-sided line having a via hole are simulated, calculated and compared, and the analysis of joule heating will involve the iterative operation mentioned in the foregoing program structure. The basic initial parameters of the circuit are listed in tables 3 and 4, wherein the negative current represents the current flowing into the corresponding end point, the circuit structure diagrams and the corresponding simulation calculation results are shown in fig. 10 to 14, fig. 11 and 12 are the joule heating rates per unit volume of the upper and lower surfaces, respectively, and fig. 13 and 14 are the temperature distributions of the upper and lower surfaces, respectively.
TABLE 3 two-sided line basic parameters
Size of substrate | Thickness of copper foil | Via hole wall thickness | Pore diameter | Pad1 | Pad2 | Pad3 | Via1(2nd layer) | Via2 |
70x34x1.6(mm) | 18μm | 25.4μm(1mil) | 2mm | -3A | -5A | -3A | -3A | 12V |
TABLE 4 double-sided line simulation calculation parameters
Therefore, the coupling thermal analysis operation result based on the method keeps higher consistency with the COMSOL operation result, and the effectiveness of an iterative operation mode related to joule heating of a line is verified.
The thermal analysis method for the double-sided PCB structure is mainly based on a Fourier series temperature analytical solution and a finite volume numerical value dispersion method, and can realize that only the metal layer and the corresponding surface area of the insulating layer are dispersed according to a two-dimensional circuit diagram, but not the whole PCB structure is dispersed; by carrying out numerical value discrete approximation on the heat balance equation of the metal layer, the linear relation between the temperatures of the metal layer units is utilized to approximately represent the heat diffusion inside the metal layer units and the heat conduction process between the metal layer layers through the metal via holes; and the heating rate of the device is calculated in the coupling solution by counting the thermal resistance parameters of the device.
The method can compensate and correct the thermal boundary condition related to the assumption of the contact area of the device, correct the homogeneity assumption of the insulating layer of the PCB substrate and the thermal boundary condition including the surface heat flux density distribution, and finally obtain the average temperature of the core of the device and the temperature distribution of other hot spot areas. Programming realization and simulation calculation of the method are carried out in an MATLAB environment, and the theoretical accuracy of the modeling method is verified by comparing with the COMSOL modeling operation result.
While the embodiments of the invention have been described in detail in connection with the drawings, it should not be construed as limiting the scope of the invention. Various modifications and changes may be made by those skilled in the art without inventive step within the scope of the appended claims.
Claims (5)
1. A steady-state thermal analysis method for a double-sided PCB structure is characterized by comprising the following steps:
s1, constructing a stable thermal balance Laplace equation of the double-sided PCB insulating layer, and solving Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation, wherein the equation comprises the following steps:
constructing a stable heat balance Laplace equation of the double-sided PCB insulating layer, and obtaining boundary conditions of the periphery, the upper surface and the lower surface of the PCB insulating layer:
wherein T is the amount of temperature change from ambient temperature or the temperature at ambient temperature of 0 ℃, qiu(x, y) and qid(x, y) is a heat flow density distribution function introduced by the upper and lower surfaces, huAnd hdHeat transfer coefficients of the upper and lower surfaces of the insulating layer, kiIs the thermal conductivity of the insulating layer, Lx、Ly、LinThe lengths of the insulating layer in the x direction and the y direction and the thickness of the insulating layer in the z direction are respectively measured under a Cartesian coordinate system;
the expression form of Fourier series analytic solution corresponding to the steady-state heat balance equation is as follows:
wherein, C1,C2,Cn,m,Cγn,mIs a Fourier series coefficient, betan、μm、γn,mCharacteristic values in a Fourier series;
decomposing a heat balance equation and a boundary condition by adopting a heat effect superposition principle:
wherein theta and eta are superposition component variables of the insulating layer heat balance equation solution;
substituting Fourier series solution into lower surface z corresponding to theta as LinThe relation between partial coefficient and characteristic value in solution is obtained under the condition of thermal boundary, and when h is reacheddWhen not zero, C2And Cγn,mThe corresponding expression is:
wherein HuAnd HdDeriving intermediate parameters of the process for the analytic solution; when h is generateddIs zero, C1,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 cos2(βnx)cos2(μmy) is not zero after integration, and the right side of the equation is used for each heating unit q in the heating areaiu,iIntegrating and then summing, and obtaining after operationCoefficient C1And Cn,m:
Wherein d isqIntegration Domain Sq for grid cell Length after discretization of surface Heat Source gridiu,iI.e. representing an incoming heat flow density of qiu,iSurface discrete unit area of δnAnd deltamTo resolve dimensionless intermediate parameters of the derivation process, qiu,1~qiu,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, qiu,1~qiu,N1I.e. the column vector q of the heat flow density is constructediuThen, 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 Ru(x,y,z)Expanding into a thermal resistance matrix R according to the coordinates of the corresponding distribution points of the regionuWherein each row is corresponding to a thermal resistance row vector for solving the temperature at a certain point, and then R is addeduMultiplying the temperature distribution vector by the corresponding heat flow density vector q to obtain a temperature distribution vector of the specified area, namely the heat effect generated by the corresponding q in the specified area; finally, superposing Fourier series analytic solutions in a matrix form corresponding to theta and eta of the specified region to obtain the temperature distribution of the region;
thus, the matrix equation corresponding to the analytic solution temperature distribution vector of the upper and lower surfaces of the insulating layer in the double-sided PCB will have the following form:
wherein, the distribution vectors of the heat flux density transmitted into the insulating layer by the upper and lower surface metal layers are qM,u、qM,dWherein q isuAnd q isdRespectively the heat flux density distribution vector q transmitted between the upper and lower surface heating components and the covering region thereofuComprising qiu(x, y) heat flux density distribution after dispersion, qdComprising qid(x, y) a heat flow density distribution after discretization; rMu,uIs qM,uAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMu,dIs qM,dAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMd,uIs qM,uAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rMd,dqM,dAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; riu,uIs quAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; riu,dIs qdAt TuA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rid,uIs quAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; rid,dIs qdAt TdA thermal resistance coefficient matrix corresponding to the thermal effect is generated in the limited area; due to quAnd q isdThe 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 equationiu(x, y) and qid(x, y) heat flow transferred between the upper and lower surface heat-generating components associated therewith and the surface areas of the insulating layer directly covered thereby, Riu,u、Riu,d、Rid,uAnd Rid,dWill include a partial all-zero column vector;
s2, calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the boundary condition of the substrate insulating layer, and calculating the heating rate of the component as the known heat source condition into a matrix equation expression of an analytic solution;
s3, approximating a heat balance equation of a discrete metal layer based on a finite volume method, summarizing to obtain a metal layer heat diffusion lumped matrix equation, and counting the heat conduction between the metal layer and a surface heating device thereof and the heat conduction of the metal via hole;
s4, constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the corrected insulating layer temperature analytic solution matrix equation and the metal layer thermal diffusion matrix equation based on numerical value dispersion;
s5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and inducing to obtain a lumped matrix equation describing a linear relation between unit potentials of the metal layer, and calculating potential distribution in the metal layer, so as to further calculate current density distribution and Joule heating distribution of the metal layer, and calculate the Joule heating distribution into a metal layer thermal diffusion matrix equation, and perform iterative calculation with temperature distribution to calculate the temperature distribution of the surface of each layer of the PCB under the condition of calculating Joule heating of the circuit;
and S6, verifying the calculation accuracy of the thermal analysis method of the double-sided PCB structure by adopting a COMSOL modeling operation result.
2. The steady-state thermal analysis method for the double-sided PCB structure as recited in claim 1, wherein the step S2 includes the steps of calculating the thermal effect of the component based on the thermal resistance parameter of the component, correcting the thermal boundary condition of the substrate insulating layer, and obtaining the corrected insulating layer temperature analytic solution matrix equation expression, which includes:
setting the heat flux density vector of the covering area of the upper surface heating device as qh,u:
qh,u=qu+huTJ,u
Wherein, TJ,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 environmentuThe calculated heat exchange under the condition, but the part of the heat exchange is taken into account in the analytic solution derivation process of the step S1, so the occupied area of the component is determined according to huCompensating for the multi-calculated heat exchange; in addition, can benefitThe thermal resistance parameter provided by the component is used for deducing the heat transfer between the relevant component and the PCB surface, and the specific deduction process is as follows:
the two formulas are combined to obtain:
qh,u=CP,uqJ,u-CT,uTJ,u
wherein the temperature distribution of the upper surface element device surface is TtopThe heat transfer coefficient of the surface to the outside is htopThe equivalent heat conduction corresponding to the heat sink existing on the surface thereof is also taken into accounttop(ii) a Mean temperature of the core is Tj(ii) a The average thermal resistance between the surface and the core is psijt(ii) a The average thermal resistance between the junction center and the bottom of the device is psijc(ii) a Equivalent thermal resistance of solder mask coating is psic,u(ii) a The heat flux density transferred from the core to the surface of the component is qjt(ii) a The heating volume density vector corresponding to the upper surface element device is PJ,u(ii) a The diagonal matrix composed of height parameters of the components is CJ,u(ii) a The density vector of the heating surface corresponding to the upper surface element device is qJ,u;CP,uFor taking into account the thermal resistance parameter of the deviceJ,uThe coefficient of (a); cT,uFor taking into account the corresponding T after the thermal resistance parameter of the deviceJ,uThe coefficient of (a);
with qh,uAlternative upper surface temperature distribution analytical solution TuQ in (1)u:
Wherein, TJ,dThe temperature distribution vector of the coverage area of the lower surface heating device is obtained; riu,uCP,uqJ,uThe term is a vector of known constants, and TuIncluding TJ,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 equationuInstead of TJ,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:
CJu,uTu=Riu,uCT,uTJ,u
wherein, CJu,uTo realize TJ,uAnd TuAnd (3) equivalently transforming the coefficients of the related terms to further obtain:
Tu=RMu,uqM,u+RMd,uqM,d+Rid,uqd+Riu,uCP,uqJ,u-CJu,uTu
similarly, when there is a device on the lower surface, the heat flow transferred between the device on the lower surface and the surface unit, which takes into account the analytical solution, is also corrected, i.e. q in the above formuladItem replacement:
wherein the heat flux density vector of the compensated lower surface covered by the heating device is qh,d(ii) a The density vector of the heating surface corresponding to the lower surface component is qJ,d;CP,dFor taking into account the thermal resistance parameter of the deviceJ,dThe coefficient of (a); cJd,uFor realizing the temperature distribution T of the covering area of the lower surface heating deviceJ,dAnd TdCoefficients of the correlation term equivalent transform;
in the same way, the T corrected by the thermal boundary condition can be obtaineddCorresponding analytic solution matrix equation, and after correction, TuAnd TdThe corresponding simultaneous analytical solution matrix equation becomes:
wherein, CJu,d、CJd,dRespectively for realizing T in the lower surface analytical solution matrix equationJ,uAnd TuRelated item and TJ,dAnd TdThe correlation term is equivalent to the transformed coefficients.
3. The steady-state thermal analysis method for the double-sided PCB structure of claim 1, wherein the step S3 approximates the discrete metal layer thermal equilibrium equation based on the finite volume method, generalizes the metal layer thermal diffusion matrix equation, and includes the following steps:
dispersing the metal layer by adopting a structured square grid, and dispersing a heat balance equation of heat conduction in the metal layer based on second-order Taylor series approximation to obtain an upper surface discrete metal layer heat balance equation:
wherein k ismThe thermal conductivity of the metal layer; q. q.sc,uThe joule heating rate distribution per unit volume of the upper surface metal layer; t isc,uIs the upper surface metal cell temperature; dcThe unit side length of the structured and discrete metal unit is shown; dmIs the thickness of the metal layer; q. q.su,uIs Tc,uThe longitudinal heat flow density of the upper surface of the corresponding unit can be the density of heat flow transferred between the corresponding unit and the heating device; q. q.sd,uIs Tc,uThe longitudinal heat flux density transmitted between the lower surface of the corresponding unit and the insulating layer; t isE,uIs Tc,uThe temperature of the adjacent metal unit at the east side of the corresponding unit; t isW,uIs Tc,uThe temperature of the corresponding unit west side adjacent metal unit; t isN,uIs Tc,uThe temperature of the adjacent metal unit at the north side of the corresponding unit; t isS,uIs Tc,uThe temperature of the adjacent metal unit on the south side of the corresponding unit; psiNv,zIs the unit thermal resistance of the metal via; t isc,dvIs the temperature of the lower surface metal layer unit connected by the via holeDegree;
similarly, a discrete heat balance equation corresponding to the lower surface metal layer unit can be obtained:
wherein q isc,dThe joule heating rate distribution per unit volume of the lower surface metal layer; t isc,dIs the lower surface metal unit temperature; dcThe unit side length of the structured and discrete metal unit is shown; dmIs the thickness of the metal layer; q. q.sd,dIs Tc,dThe longitudinal heat flux density of the lower surface of the corresponding unit can be the density of heat flux transferred between the corresponding unit and the heating device; q. q.su,dIs Tc,uThe longitudinal heat flux density transmitted between the upper surface of the corresponding unit and the insulating layer; t isE,dIs Tc,dThe temperature of the adjacent metal unit at the east side of the corresponding unit; t isW,dIs Tc,dThe temperature of the corresponding unit west side adjacent metal unit; t isN,dIs Tc,dThe temperature of the adjacent metal unit at the north side of the corresponding unit; t isS,dIs Tc,dThe temperature of the adjacent metal unit on the south side of the corresponding unit; t isc,uvThe temperature of the lower surface metal layer unit connected through the via hole is shown;
and (3) inducing a thermal diffusion matrix equation of the upper and lower surface metal layers based on numerical value dispersion:
wherein q isMJ,uAnd q isMJ,dRespectively corresponding unit joule heating vectors of the upper surface metal layer and the lower surface metal layer; h istopThe heat transfer coefficient between the surface of the component and the outside is taken as the coefficient; mMV,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the upper surface metal layer unit is formed; mMV,dThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; mV,dDiscrete thermal equilibrium equation from top surface metal layer unitThe temperature coefficient of the middle corresponding lower surface metal layer unit is formed; mV,uThe coefficient corresponding to the temperature of the upper surface metal layer unit in the discrete heat balance equation of the lower surface metal layer unit is formed; due to qh,u、qh,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 CM,u、CM,dWill also include partial all-zero column vectors;
q is to beh,uAnd q ish,dSubstituting the expression of (a) into the thermal diffusion matrix equation to obtain:
wherein, CMu,uAnd CMd,dRespectively for realizing T in the above equationJ,uAnd TuRelated item and TJ,dAnd TdThe correlation term is equivalent to the transformed coefficients.
4. The steady-state thermal analysis method for the double-sided PCB structure of claim 1, wherein the step S4 is implemented for constructing a thermal analysis coupling equation set of the double-sided PCB structure according to the modified insulation layer temperature analytic solution matrix equation and the metal layer numerical value discrete-based thermal diffusion matrix equation, and comprises:
the thermal analysis coupling equation set of the double-sided PCB structure is as follows:
wherein, Tu,Td,qM,u,qM,dIf the vector is unknown, and other vector and matrix coefficients are known, the unknown temperature vector and the unknown heat flux density vector can be obtained by combining the four matrix equations.
5. A method for steady-state thermal analysis of a dual-sided PCB structure according to claim 1, wherein in step S5, a discrete metal layer current continuity equation is approximated based on a finite volume method, a lumped matrix equation describing a linear relationship between unit potentials of the metal layer is obtained by calculation and induction, and a potential distribution in the metal layer is obtained, so that a current density distribution and a joule heating distribution of the metal layer can be further obtained, and the joule heating distribution can be included in a metal layer thermal diffusion matrix equation, and an iterative calculation can be performed between the joule heating distribution and the temperature distribution, so as to obtain a temperature distribution of a surface where each layer of the PCB is located under the condition of the joule heating of the circuit, including:
based on second-order Taylor series approximation, obtaining an approximate discrete equation corresponding to a current continuity equation of a square area which is surrounded by the sections of Pn, Pw, Ps and Pe and takes the vertex Pc of the metal layer unit as the center:
wherein, VPc、VPE、VPW、VPNAnd VPSPc and P are respectivelyE、PW、PNAnd the potential at PS; wherein σPe、σPw、σPnAnd σPsElectrical conductivities at Pe, Pw, Pn and Ps, respectively;
the above is to apply the finite volume method to the control volume centered at the vertex;
the electrical conductivity of a metal varies with the temperature of the material, and its relationship to temperature can be represented by the following equation:
where ρ is the resistivity, ρ0For resistivity at initial ambient temperature, Δ TMTemperature rise of discrete metal units compared to ambient temperature, αTAlpha for copper, for a material dependent temperature coefficientTIs positive, i.e. its resistanceThe rate increases with increasing temperature;
the conductivity of the midpoint of the cross section is considered to be the conductivity corresponding to the equivalent resistance of the synthesized metal units at the left side, the right side, the upper side and the lower side, namely the reciprocal of the sum of the resistivities of the two units, sigmaPsConductivity σ Using C1 cell and C2 cellC1And σC2Represents:
therefore, the electrical conductivity of each unit of the metal layer is calculated according to the temperature distribution obtained in step S4; the conductivity corresponding to Pc can be approximately found by applying the conductivities corresponding to the four midpoints surrounding Pc:
for 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 Pm1、Pm2、Pm3And Pm4The electrical conductivity of (d); vPvt1、VPvt2、VPvt3And VPvt4Potentials at Pvt1, Pvt2, Pvt2 and Pvt2, respectively;
during the PCB design process, 45 degree inclined lines are designed, but a square grid can be reconstructed, and the discrete approximate current continuity equation of a square control volume unit with a quarter cell area and Ptr2 as a vertex at the line boundary is as follows:
wherein σPmeAnd σPmnAre each PmeAnd PmnThe electrical conductivity of (d); vPtr2、VPc2And VPc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation shows that the current continuity equation is the same as that corresponding to the triangular control unit surrounding Ptr2 after filling the sawtooth gaps of the boundary, i.e. it proves that the partial sawtooth boundaries formed discretely by the square grid and the boundary of the 45-degree inclined line can be approximately equivalent in current continuity;
while the control volume units at the line corners are not directly equivalent to the actual line corners, the control volume units with the cross sections of Pmn3, Pmw3 and Pms3 surrounding the boundary vertex Pb3 have the length d of the cross section of the lower Pms3 in the sawtooth boundarycA/2, while the length in the actual line is dcTherefore, 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 Pmw3、Pmn3And Pms3The electrical conductivity of (d); vPb3、VPb4、VPc4And VPc5Potentials at Pb3, Pb4, Pc4 and Pc5, respectively;
for the metal layer units connected by the via holes, the corresponding approximately discrete current continuity equation can also be obtained according to the above vertex-centered finite volume method, and assuming that vertices of adjacent metal layer units, of which the metal layer unit vertex Pc is connected by the via holes, are Pu and Pd, the corresponding approximately discrete current continuity equation becomes:
wherein, VPuAnd VPdAre each PuAnd PdThe potential of (d); sigmaPuAnd σPdElectrical conductivity at Pu and Pd, respectively;
finally, a summary matrix equation of approximate discrete current continuity equations corresponding to the upper and lower surface metal layers can be obtained:
wherein, V1And V2Respectively corresponding to the summary vectors of all the discrete unit top potential with unknown potential of the upper surface metal layer and the lower surface metal layer; mσ,1And Mσ,2Respectively V in the upper surface matrix equation1With V in the lower surface matrix equation2The coefficient matrix of (a), is formed by conductivity parameters; kp,1And Kp,2Known terms in two matrix equations are respectively included, and corresponding terms of known endpoint potential or section current are included; mσV,1_2Is V in the upper surface matrix equation2The coefficient matrix of (a) is formed by conductivity parameters corresponding to the lower surface metal layer unit connected by the via hole and the upper surface metal unit discrete heat balance equation, for V2In cells not connected by vias, MσV,1_2The corresponding column in (1) is an all-zero column; mσV,2_1For V in the lower surface matrix equation1The coefficient matrix of (a) is formed by conductivity parameters corresponding to the upper surface metal layer unit connected by the via in the discrete thermal equilibrium equation of the lower surface metal unit, for V1In cells not connected by vias, MσV,2_1The corresponding column in (1) is an all-zero column; the above matrix equation can therefore also incorporate the approximate discrete current continuity equation for metal cells not connected to a via; v can be obtained by simultaneously solving the two matrix equations1And V2;
After the peak potential distribution in the line is found, each square can be foundPotentials at the center of the cell and at the center of the cross-section around the cell were determined by finite volume method using a small diamond region around the center of the cell C2 to determine the potential V at the center of the cell C2C2:
Wherein, VPvt4Potential at cell C2 apex Pvt 4; similarly, the potential V of the point Ps at the center of the cross section can be obtainedPs:
Wherein, VC1Potential at the center of cell C1;
thus, the lateral current density and the vertical current density corresponding to each cell center can be approximated, and the lateral current density J corresponding to the cell center can be approximated for the cell C1x,C1And longitudinal current density Jy,C1Can be expressed as:
then C1 corresponds to a total current density of:
for the triangular unit at the oblique line boundary, assuming that the potential is linearly distributed in the region, the current density of the corresponding triangular region is obtained, and the linear potential distribution corresponding to the triangle formed by three points Ptr2, Ptr1 and Pc1 can be expressed as follows:
wherein, the third-order square matrix in the equation is composed of operation formulas corresponding to the coordinates of three vertexes under a Cartesian coordinate system, SΔIs a triangular area, and a coefficient a corresponding to the potential linear distribution expression can be obtained0、a1And a2Then, the transverse current density J within the triangle is knownx,ΔAnd longitudinal current density Jy,ΔComprises the following steps:
wherein σΔConductivity corresponding to the triangular region; on the premise of obtaining the current density corresponding to the metal cell, the joule heating rate per unit volume can be obtained, for example, for cell C1, the joule heating rate per unit volume q can be obtainedC1Can be expressed as:
qC1=EC1·JC1=(Jx,C1 2+Jy,C1 2)/σC1=JC1 2/σC1
wherein E isC1The 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.sC1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the upper and lower surface metal layersMJ,uAnd q isMJ,d。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010572678.2A CN111723486B (en) | 2020-06-22 | 2020-06-22 | Double-sided PCB structure steady-state thermal analysis method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010572678.2A CN111723486B (en) | 2020-06-22 | 2020-06-22 | Double-sided PCB structure steady-state thermal analysis method |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111723486A CN111723486A (en) | 2020-09-29 |
CN111723486B true CN111723486B (en) | 2021-05-04 |
Family
ID=72569916
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010572678.2A Active CN111723486B (en) | 2020-06-22 | 2020-06-22 | Double-sided PCB structure steady-state thermal analysis method |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111723486B (en) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117094281B (en) * | 2023-10-19 | 2024-02-13 | 杭州行芯科技有限公司 | Method for obtaining thermodynamic parameter, electronic equipment and storage medium |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102955477A (en) * | 2012-10-26 | 2013-03-06 | 南京信息工程大学 | Attitude control system and control method of four-rotor aircraft |
CN104237300A (en) * | 2014-08-27 | 2014-12-24 | 北京时代民芯科技有限公司 | Fixture and method for testing steady state thermal resistance of glass sealed surface-mount diode |
WO2015104715A1 (en) * | 2014-01-13 | 2015-07-16 | Corning Optical Communications Wireless Ltd. | Dissipating heat from electronic devices |
CN107024838A (en) * | 2017-06-05 | 2017-08-08 | 西华大学 | A kind of PCB exposure machines lamp box, PCB exposure machines and reset roll adjustment method |
CN107887669A (en) * | 2017-11-07 | 2018-04-06 | 大连理工大学 | A kind of heat dissipation metal electrokinetic cell pack arrangement design method and battery bag |
CN109726777A (en) * | 2017-10-30 | 2019-05-07 | 东莞理工学院 | PCB appearance detection system and detection method Internet-based |
-
2020
- 2020-06-22 CN CN202010572678.2A patent/CN111723486B/en active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102955477A (en) * | 2012-10-26 | 2013-03-06 | 南京信息工程大学 | Attitude control system and control method of four-rotor aircraft |
WO2015104715A1 (en) * | 2014-01-13 | 2015-07-16 | Corning Optical Communications Wireless Ltd. | Dissipating heat from electronic devices |
CN104237300A (en) * | 2014-08-27 | 2014-12-24 | 北京时代民芯科技有限公司 | Fixture and method for testing steady state thermal resistance of glass sealed surface-mount diode |
CN107024838A (en) * | 2017-06-05 | 2017-08-08 | 西华大学 | A kind of PCB exposure machines lamp box, PCB exposure machines and reset roll adjustment method |
CN109726777A (en) * | 2017-10-30 | 2019-05-07 | 东莞理工学院 | PCB appearance detection system and detection method Internet-based |
CN107887669A (en) * | 2017-11-07 | 2018-04-06 | 大连理工大学 | A kind of heat dissipation metal electrokinetic cell pack arrangement design method and battery bag |
Non-Patent Citations (2)
Title |
---|
Generalized Heat Transition Matrix for Arbitrarily Shaped Thermal Media and Its Applications to Steady-State Heat Conduction Problems in Large-Scale Systems;Gaobiao Xiao 等;《Numerical Heat Transfer》;20111231;第319-338页 * |
PCB及元件的温度场有限元分析;李晓明 等;《北京航空航天大学》;20000229;第26卷(第1期);第5-7页 * |
Also Published As
Publication number | Publication date |
---|---|
CN111723486A (en) | 2020-09-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Lu et al. | Transient electrical-thermal analysis of 3-D power distribution network with FETI-enabled parallel computing | |
Adam | New correlations between electrical current and temperature rise in PCB traces | |
Zhang | Improved numerical-analytical thermal modeling method of the PCB with considering radiation heat transfer and calculation of components’ temperature | |
Zhang et al. | A modeling methodology for thermal analysis of the PCB structure | |
CN111723486B (en) | Double-sided PCB structure steady-state thermal analysis method | |
Monier-Vinard et al. | Analytical modeling of multi-layered Printed Circuit Board dedicated to electronic component thermal characterization | |
CN111859622B (en) | Single-sided PCB structure steady-state thermal analysis method | |
Košel et al. | FEM simulation approach to investigate electro-thermal behavior of power transistors in 3-D | |
Li et al. | DC IR-drop analysis of multilayered power distribution network by discontinuous Galerkin method with thermal effects incorporated | |
CN111723507B (en) | Multilayer PCB structure steady-state thermal analysis method | |
Lall et al. | Thermal design rules for electronic components on conducting boards in passively cooled enclosures | |
Bashkirov et al. | The influence of a 3D model of a radio electronic component on thermal simulation | |
US7596482B2 (en) | System and method to analyze and determine ampacity risks on PCB interconnections | |
Rogié et al. | Practical analytical modeling of 3D multi-layer Printed Wired Board with buried volumetric heating sources | |
Shao et al. | Thermal-aware DC IR-drop co-analysis using non-conformal domain decomposition methods | |
Nowak et al. | A thermal network approach for a quick and accurate study of multiple connected switchgears | |
Khairnasov | Modeling and thermal analysis of heat sink layers of multilayer board | |
KR et al. | Testing of current carrying capacity of conducting tracks in high power flexible printed circuit boards | |
Zhang et al. | A thermal model for the PCB structure including thermal contribution of traces and vias | |
Guofeng et al. | The research of thermal design for vehicle controller based on simulation | |
Dogruoz et al. | Spatial variation of temperature on printed circuit boards: Effects of anisotropic thermal conductivity and Joule heating | |
Summers et al. | Design of data centre rack arrangements using open source software | |
Weiss et al. | Modeling air flow in electronic packages | |
Beckermann et al. | Use of a two-dimensional simulation model in the thermal analysis of a multi-board electronic module | |
Felczak et al. | Optimal placement of electronic devices in forced convective cooling conditions |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |