CN111859622B - Single-sided PCB structure steady-state thermal analysis method - Google Patents

Abstract

The invention discloses a steady-state thermal analysis method for a single-sided PCB structure, which comprises the following steps: s1, constructing a stable thermal balance Laplace equation of the PCB insulating layer to obtain an insulating layer temperature analytic solution matrix equation; s2, correcting boundary conditions of the insulating layer, and taking the heating rate of the component as a known heat source condition to be recorded into a matrix equation expression of an analytic solution; s3, summarizing to obtain a metal layer thermal diffusion lumped matrix equation, and recording the heat conduction between the metal layer and the surface heating device; s4, constructing a thermal analysis coupling equation set of the single-sided PCB structure; s5, calculating and summarizing to obtain a lumped matrix equation, solving the potential distribution in the metal layer, further solving the current density distribution and the Joule heating distribution of the metal layer, calculating the Joule heating distribution into the metal layer thermal diffusion matrix equation, and performing iterative calculation with the temperature distribution to obtain the temperature distribution of the PCB surface calculated by the Joule heating of the circuit; and S6, verifying the calculation accuracy of the thermal analysis method of the single-sided PCB structure by adopting a COMSOL modeling operation result.

Description

Single-sided PCB structure steady-state thermal analysis method

Technical Field

The invention belongs to the technical field of PCB thermal analysis, and particularly relates to a steady-state thermal analysis method for a single-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 approximately equivalent thermal conductivity related to the local copper coverage of the PCB in algorithmic principle, while the latter is based on 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.

However, no relevant software for PCB thermal analysis exists in China at present.

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 the PCB. 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, which approximates continuous unknown function to finite sub-functions (generally in the form of polynomial containing unknown variables which can be approximated to linear change in the discrete domain) effective in the discrete control domain in discrete finite space or time, then substitutes the finite sub-functions into the original continuous equation to establish the relationship between the corresponding variables of each discrete sub-domain, and then uses the given boundary condition to obtain the approximate distribution of the unknown variables. 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 by the students of the university of great continental engineering, Song Xueguan, et al, combined with the university of neckasei, lester, etc., on the experience of using multi-physics analysis software such as ANSYS, COMSOL, etc. in the industrial and academic circles show that over 60% of users are satisfied with the operational accuracy of such software, but that over 60% of users consider the operational speed of such software to be 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 single-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 single-sided PCB structure comprises the following steps:

s1, constructing a steady-state heat balance Laplace equation of the single-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, inducing to obtain a metal layer heat diffusion lumped matrix equation, and calculating the heat conduction between the metal layer and a surface heating device;

s4, constructing a thermal analysis coupling equation set of the single-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, thereby further calculating current density distribution and Joule heating distribution of the metal layer, calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with temperature distribution to calculate the temperature distribution of the PCB surface under the premise of calculating Joule heating of the circuit.

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

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

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

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

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

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

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

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

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

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

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

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

therefore, a single-sided PCB insulating layer temperature analytic solution matrix equation is obtained:

T_u＝R_Mq_M+R_Iq_u

wherein the heat flux density distribution vector of the metal layer transmitted into the insulating layer is q_M(ii) a And q is_uA heat flux density distribution vector, q, transferred between the surface heating component and its covered region_uComprising q_iu(x, y) a heat flow density distribution after discretization; r_MIs q_MA thermal resistance matrix corresponding to a thermal effect generated on the surface of the insulating layer; r_IIs q_uActing on insulating layerA corresponding thermal resistance matrix; due to q_uThe heat flux transferred through the surface of the metal layer is taken into account in the thermal diffusion matrix equation corresponding to the metal layer, and the analytical solution matrix equation only takes into account the heat flux transferred between the surface heating element and the surface region of the insulating layer directly covered by the surface heating element, so that R_IWill include a partial all zero column vector.

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

let the heat flux density vector of the covering area of the surface heating device be q_h,u：

q_h,u＝q_u+h_uT_J,u

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

the two formulas are combined to obtain:

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

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

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

Wherein R is_IC_P,uq_J,uThe term is a vector of known constants, and T_uIncluding T_J,uAll elements of the vector can be transformed by adding zero to the matrix coefficients in a rising dimension way to realize the use of T in the matrix equation_uInstead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_JT_u＝R_IC_T,uT_J,u

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

T_u＝R_Mq_M+R_IC_P,uq_J,u-C_J,uT_u。

preferably, the step S3 approximates the discrete metal layer thermal equilibrium equation based on the finite volume method, generalizes the lumped matrix equation for metal layer thermal diffusion, 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 a heat balance equation of the discrete metal layer:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_cIs the joule heating rate distribution per unit volume of the metal layer; t is_cIs the metal unit temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_c,uIs T_cThe density of heat flow transmitted from the upper surface of the corresponding unit can be the density of heat flow transmitted between the corresponding unit and the heating device; q. q.s_dIs T_cThe heat flux density of the lower surface of the corresponding unit; t is_EIs T_cThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_WIs T_cThe temperature of the corresponding unit west side adjacent metal unit; t is_NIs T_cThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_SIs T_cThe temperature of the adjacent metal unit on the south side of the corresponding unit;

and (3) inducing a metal layer thermal diffusion matrix equation based on numerical value dispersion:

M_M,uT_u＝q_MJ,u+C_M,uq_h,u-q_M

wherein q is_MJ,uThe unit joule heating vector corresponding to the metal layer; m_M,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the metal layer unit is formed; due to q_h,uThe compensated heat flow directly transmitted through the surface of the insulating layer is already counted in the analytic solution matrix equation, and only the heat flow transmitted between the heating device and the covering metal layer part of the heating device is counted in the lumped matrix equation of heat diffusion, so C_M,uWill also include partial all-zero column vectors;

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

wherein, C_Mu,uI.e. to realize T in the above equation_J,uAnd T_uThe correlation term is equivalent to the transformed coefficients.

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

wherein, T_u,q_MIf the temperature vector is an unknown vector, and other vector and matrix coefficients are known quantities, the unknown temperature vector and the unknown heat flow density vector can be obtained by combining the two 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 further calculate the joule heating distribution into a metal layer thermal diffusion matrix equation, and iteratively calculate the joule heating distribution with a temperature distribution to calculate a temperature distribution of a surface where the metal layer is located under the condition that joule heating of the circuit is included:

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

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

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

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

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

calculating the conductivity of each cell of the metal layer according to the temperature distribution obtained in the 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 positions, it may happen that the analyzed control volume cell is non-square, and the control cell number 2 surrounded by dark dotted lines has the following approximate discrete current continuity equation:

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

for a 45-degree inclined line in the PCB design process, a square grid is adopted for reconstructing the inclined line, and for a small square control volume unit with a quarter unit area and a Ptr2 vertex at the boundary of the upper left-hand line, the corresponding discrete approximate current continuity equation is as follows:

wherein σ_PmeAnd σ_PmnAre each P_meAnd P_mnThe electrical conductivity of (d); v_Ptr2、V_Pc2And V_Pc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation shows that the current continuity equation is the same as 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 boundary_cA/2, while the length in the actual line is d_cTherefore, such corners need to be individually processed in the analysis method to account for the approximate current continuity equation corresponding to their actual boundaries:

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

Finally, a summary matrix equation of an approximate discrete current continuity equation corresponding to the metal layer can be obtained:

M_σV_M＝K_p

wherein, V_MA summary vector of the discrete unit vertex potentials corresponding to all unknown potentials of the metal layer; m_σFor V in the matrix equation_MThe coefficient matrix of (a), is formed by conductivity parameters; k_pKnown terms in a matrix equation are included, including corresponding terms of known endpoint potentials or cross-sectional currents;

after the peak potential distribution in the line is found, the potentials corresponding to the centers of the square cells and the centers of the cross sections around the cells are calculated, and the potential V at the center of the cell C2 can be found by applying a finite volume method to a small diamond region surrounding the center of the cell C2_C2：

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

Wherein, V_C1Potential at the center of cell C1;

thus, the lateral and longitudinal current densities corresponding to each cell center can be approximated, e.g., for cell C1, the lateral current density J corresponding to its center_x,C1And longitudinal current density J_y,C1Can be expressed as:

then C1 corresponds to a total current density of:

for the triangular unit at the oblique line boundary, the electric potential can be assumed to be linearly distributed in the region, so as to obtain the current density of the corresponding triangular region, and the linear electric 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 obtained₀、a₁And a₂Then, the transverse current density J within the triangle is known_x,ΔAnd longitudinal current density J_y,ΔComprises the following steps:

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

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

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

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

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

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

the method for coupling the numerical value dispersion and the 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 area on the surface of the insulating layer of the substrate, can realize approximate representation of the heat diffusion in the metal layer based on the finite volume method, so as to count the boundary condition of the heat flow density distribution on the surface of the insulating layer of the substrate under the heat diffusion action of the metal layer, further realize the mutual coupling solution of the temperature analysis solution matrix equation and the metal layer heat diffusion matrix equation, and 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 schematic diagram of a discrete unit of the circuit 1.

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

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

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

Fig. 6 shows a general structure and flow of a corresponding procedure of the single-sided PCB structure coupled thermal analysis method.

FIG. 7 is a comparison of the results of thermal distribution simulation calculations for an example corresponding to the parameters of Table 5-1.

Figure 8 single line architecture.

Fig. 9 the structure of fig. 8 corresponds to the potential distribution calculation result.

Fig. 10 shows the calculation results of the current density (X direction) distribution corresponding to the structure of fig. 8.

Fig. 11 shows the calculation result of the current density (Y direction) distribution corresponding to the structure of fig. 8.

Fig. 12 is a calculation result of the total current density distribution corresponding to the structure of fig. 8.

Fig. 13 is a single line structure.

Fig. 14 is a calculation result of the total current density distribution corresponding to the structure of fig. 13.

FIG. 15 is a calculated Joule heat distribution for the structure of FIG. 13.

Fig. 16 is a calculation result of the line temperature distribution corresponding to the structure of fig. 13.

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 single-sided PCB structure comprises the following steps:

s1, constructing a steady-state heat balance Laplace equation of the single-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, inducing to obtain a metal layer heat diffusion lumped matrix equation, and calculating the heat conduction between the metal layer and a surface heating device;

s4, constructing a thermal analysis coupling equation set of the single-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, thereby further calculating current density distribution and Joule heating distribution of the metal layer, calculating the Joule heating distribution into a metal layer thermal diffusion matrix equation, and performing iterative calculation with temperature distribution to calculate the temperature distribution of the PCB surface under the premise of calculating Joule heating of the circuit.

And S6, verifying the calculation accuracy of the thermal analysis method of the single-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 assumption needs to be carried out on the structural characteristics of the single-sided PCB:

the single-sided PCB is mainly composed of a core insulating substrate layer and a metal layer covered on the surface of the core insulating substrate layer. A common substrate layer material is FR-4 glass epoxy resin, which has a thermal conductivity of 0.3W/mK and a crystallization temperature of about 120 ℃, and when the temperature is exceeded, both the mechanical properties and the insulating properties of FR-4 are significantly reduced. The thickness of the substrate insulating layer is generally much larger (>0.5mm) 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. The reporting algorithm mainly realizes thermal analysis of a single-sided PCB structure consisting of an insulating layer and a metal layer, and the influence of the silk-screen layer on heat transfer is ignored in the method because the area of the silk-screen layer on the surface of the PCB is low; in addition, since the solder resist coating is thin (typically 25.4 μm) and has low thermal conductivity (0.25W/mK), its equivalent thermal resistance can be taken into account in the heat transfer coefficient to the PCB surface and to 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 will also exist between the components and the PCB, and the component packaging 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 with the micrometer-sized thickness, which cannot be ignored in the thermal analysis. However, since the metal layer is thin and has high thermal conductivity, the temperature gradient between the upper and lower surfaces thereof can be ignored, and the insulating layer is closely attached to the metal layer, so that it can be assumed that the temperature distribution of the metal layer is the same as that of the surface region of the insulating layer in contact with the metal layer. The insulating layer is a main constituent material of the PCB, but since it is not a heat source, the temperature distribution inside the insulating layer is not a main solution target.

Therefore, in the method, the temperature distribution analytical solution of the surface of the insulating layer is used for giving 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 analytical solution; meanwhile, a numerical value dispersion method based on a finite volume method is used for calculating the thermal diffusion inside the metal layer, and an approximate steady-state thermal diffusion matrix equation is given, wherein the approximate steady-state thermal diffusion matrix equation contains the heat flow density distribution actually transmitted from the metal layer to the insulating layer, and the distribution is also used as one of the coupling variables; and finally solving the solution of the coupling variable through a simultaneous analytical solution matrix equation and a thermal diffusion matrix equation.

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

The following detailed description of the method is made under the above-mentioned assumption:

step S1, constructing a steady-state heat balance Laplace equation of the single-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:

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

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

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

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

substituting fourier series solution into the lower surface corresponding to T (z ═ L)_in) The relation between partial coefficient and characteristic value in solution is obtained under the condition of thermal boundary, and when h is reached_dWhen not zero, C₂And C_γnmThe corresponding expression is:

wherein H_uAnd H_dTo resolveSolving intermediate parameters of the derivation process; when h is generated_dIs zero, C₁,C_γn,mComprises the following steps:

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

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

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

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

therefore, a single-sided PCB insulating layer temperature analytic solution matrix equation is obtained:

T_u＝R_Mq_M+R_Iq_u (8)

wherein the heat flux density distribution vector of the metal layer transmitted into the insulating layer is q_M(ii) a And q is_uA heat flux density distribution vector, q, transferred between the surface heating component and its covered region_uComprising q_iu(x, y) a heat flow density distribution after discretization; r_MIs q_MA thermal resistance matrix corresponding to a thermal effect generated on the surface of the insulating layer; r_IIs q_uActing on the thermal resistance matrix corresponding to the insulating layer; due to q_uThe heat flux transferred through the surface of the metal layer is taken into account in the thermal diffusion matrix equation corresponding to the metal layer, and the analytical solution matrix equation only takes into account the heat flux transferred between the surface heating element and the surface region of the insulating layer directly covered by the surface heating element, so that R_IWill 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 boundary condition of the substrate insulating layer, and calculating the heating rate of the component as the known heat source condition into the matrix equation expression of the analytic solution, which specifically comprises:

covering with surface heating deviceThe heat flow density vector of the region is q_h,u：

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

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

the two formulas are combined to obtain:

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

with q_h,uAlternative surface temperature profile analytical solutionT_uQ in (1)_u：

Wherein R is_IC_P,uq_J,uThe term is a vector of known constants, and T_uIncluding T_J,uAll elements of the vector can be transformed by adding zero to the matrix coefficients in a rising dimension way to realize the use of T in the matrix equation_uInstead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_JT_u＝R_IC_T,uT_J,u (13)

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

T_u＝R_Mq_M+R_IC_P,uq_J,u-C_J,uT_u (14)

step S3, approximating the heat balance equation of the discrete metal layer based on the finite volume method, inducing to obtain the lumped matrix equation of the metal layer heat diffusion, and calculating the heat conduction between the metal layer and the surface heating device, which specifically comprises:

s3.1, approximating a heat balance equation of the discrete metal layer based on a finite volume method, which specifically comprises the following steps:

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

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

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

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

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

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

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

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

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_cIs the joule heating rate distribution per unit volume of the metal layer; t is_cIs the metal unit temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer; q. q.s_c,uIs T_cThe density of heat flow transmitted from the upper surface of the corresponding unit can be the density of heat flow transmitted between the corresponding unit and the heating device; q. q.s_dIs T_cThe heat flux density of the lower surface of the corresponding unit; t is_EIs T_cThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_WIs T_cThe temperature of the corresponding unit west side adjacent metal unit; t is_NIs T_cThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_SIs T_cThe temperature of the adjacent metal unit on the south side of the corresponding unit;

therefore, the temperature relationship between adjacent units is established by the equation after dispersion, so that the subsequent dispersion heat balance equations corresponding to all metal layer units can be conveniently summarized into a lumped matrix equation for representing the heat diffusion of the metal layer. 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, inducing a metal layer thermal diffusion matrix equation based on numerical value dispersion, which specifically comprises the following steps:

M_M,uT_u＝q_MJ,u+C_M,uq_h,u-q_M (21)

wherein q is_MJ,uThe unit joule heating vector corresponding to the metal layer; m_M,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the metal layer unit is formed; due to q_h,uThe compensated heat flow directly transmitted through the surface of the insulating layer is already counted in the analytic solution matrix equation, and only the heat flow transmitted between the heating device and the covering metal layer part of the heating device is counted in the lumped matrix equation of heat diffusion, so C_M,uWill also include partial all-zero column vectors;

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

wherein, C_Mu,uI.e. to realize T in the above equation_J,uAnd T_uCoefficients 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 single-sided PCB structure, which specifically comprises the following steps:

wherein, T_u,q_MIf the temperature vector is an unknown vector, and other vector and matrix coefficients are known quantities, the unknown temperature vector and the unknown heat flow density vector can be obtained by combining the two matrix equations.

Step S5, approximating a discrete metal layer current continuity equation based on a finite volume method, calculating and summarizing a lumped matrix equation describing a linear relationship between unit potentials of the metal layer, and calculating a potential distribution in the metal layer, thereby further calculating a current density distribution and a 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 a temperature distribution to calculate a temperature distribution of a surface where the metal layer is located under the condition of joule heating of the line, which specifically includes:

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

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

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

for a square area surrounded by the cross sections of Pn, Pw, Ps and Pe with the metal layer unit vertex Pc 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 area:

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

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

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

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

therefore, an approximate discrete equation corresponding to a current continuity equation of a square area surrounded by the cross section of the metal layer unit with the vertex Pc as the center and Pn, Pw, Ps and Pe can be obtained as shown in fig. 3:

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

the 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, ρ₀In the initial ringResistivity at ambient temperature, Δ T_MTemperature rise of discrete metal units compared to ambient temperature, α_TAlpha for copper, for a material dependent temperature coefficient_TIs positive, i.e. its resistivity increases with increasing temperature;

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

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

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

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

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

wherein σ_PmeAnd σ_PmnAre each P_meAnd P_mnThe electrical conductivity of (d); v_Ptr2、V_Pc2And V_Pc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation 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 jagged gaps of the boundary. That is, it is proved that the partial zigzag boundary formed discretely in the square grid and the boundary of the 45-degree inclined line are approximately equivalent in current continuity.

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

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

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:

finally, a summary matrix equation of an approximate discrete current continuity equation corresponding to the metal layer can be obtained:

M_σV_M＝K_p (36)

wherein, V_MA summary vector of the discrete unit vertex potentials corresponding to all unknown potentials of the metal layer; m_σFor V in the matrix equation_MThe coefficient matrix of (a), is formed by conductivity parameters; k_pKnown terms in the matrix equation, such as corresponding terms describing known endpoint potentials or cross-sectional currents, are included;

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

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

Wherein, V_C1Potential at the center of cell C1;

thus, the lateral and longitudinal current densities corresponding to each cell center can be approximated, e.g., for cell C1, the lateral current density J corresponding to its center_x,C1And longitudinal current density J_y,C1Can be expressed as:

then C1 corresponds to a total current density of:

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

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

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

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

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

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

joule (J)The heating rate is related to the conductivity, and the conductivity is influenced by the temperature, so when the joule heating of the line is counted, iterative operation between the temperature distribution and the joule heating rate distribution is required, when the error of the iterative operation is smaller than a set error, the temperature distribution and the joule heating rate distribution can be approximately calculated, and the error can be set as the maximum value of the allowable unit temperature change in two iterative calculations. q. q.s_C1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the metal layer_MJ，u。

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. 6.

The general flow structure for realizing thermal analysis of the PCB structure comprises the following steps:

the general structure and flow of the PCB structure thermal analysis program written by the method of the present invention are shown in FIG. 6 below, 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 generated, wherein the coefficient matrix comprises a coefficient matrix which is required by calculating, analyzing and solving a thermal resistance coefficient matrix and a numerical discrete equation and is included in the heat conduction of a metal layer and the 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 system by calling each coefficient matrix and a known heating density vector. If the joule heating of part of the metal layer circuit needs to be taken into account, the influence of the temperature on the metal conductivity needs to be taken into account, so that iterative operation between the joule heating and the temperature distribution is needed to obtain the steady-state circuit temperature distribution. In the iterative calculation part, the Joule heating density distribution of the circuit is calculated according to the assumed temperature distribution of the circuit, then the coupling solution of the corresponding temperature distribution is calculated, then the temperature distribution is used for calculating the new Joule heating density distribution of the circuit, and 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 set error limit value.

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

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. Previous studies 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 a multi-branch line, is that the sum of the known branch currents or the number of branch terminal potentials is equal to the number of line terminals, and that at least one terminal potential is known.

a. Examples of thermal diffusion of Metal layers

First, the importance of the thermal diffusion effect of the metal layer in the thermal analysis of the PCB structure is illustrated by way of example, and the following table shows the relevant parameters of the single-sided PCB simplified structure with two power amplifier chips.

Table 1 double power amplifier PCB structure parameter table

Base size (mm)	149x 98x 1.6
		Copper foil thickness (mum)	35
Heat generation rate (W)	1.5 (upper)/2 (lower)
		Upper and lower surface HTC (W/m)2K)	10

The above structure was subjected to thermal analysis calculations using the MATLAB program based on the method of the present invention, without/with the metal layer thermal diffusion being taken into account, to obtain the temperature distribution results at an ambient temperature of 0 ℃, as shown in fig. 7. It can be seen that the hot spot temperature would reach 220 ℃ without accounting for the thermal diffusion of the metal layer, whereas the hot spot temperature would be only 72 ℃ after accounting for the thermal diffusion of the metal layer. Although the thickness of the metal layer is only 35 μm, the effect of the metal layer on the heat diffusion of the auxiliary device is very remarkable, so that the heat diffusion of the metal layer is necessary to be added in the thermal analysis of the PCB.

b. Example of Metal line Current Density distribution calculation

When the current in the metal layer line is large, there is a possibility that high joule heat is generated in the line to cause a high temperature rise in the line. From the foregoing analysis, it can be seen that the joule heat distribution is obtained on the premise that the calculation of the current density distribution is completed, and in order to verify the rationality of the method for calculating the current density distribution in the method, a simulation calculation is performed on the circuit structure shown in fig. 8, the relevant parameters of the circuit are listed in table 2, and the relevant coupled thermal analysis simulation calculation results and the modeling analysis results of the circuit in the COMSOL software are presented in fig. 9 to 12, where fig. 9 is a potential distribution, fig. 10 is a current density distribution in the x direction, fig. 11 is a current density distribution in the y direction, and fig. 12 is a total current density distribution.

Table 2 circuit parameter table of fig. 8

Base size (mm)	239x 121x 1.6
		Copper foil thickness (mum)	18
Line current (A)	17.5A
		Number of discrete units of line	4892
Ambient temperature (. degree. C.)	20
		Number of discrete units of COMSOL model	61298

It can be seen from the simulation calculation results that the current density calculation method of the present invention realizes more accurate analysis of the current density of the double-ended line on the premise of less discrete unit number, and the simulation calculation results and the modeling analysis results of the COMSOL software according to the method present more consistent distribution trend and smaller numerical error (the error of the total current density peak value is within 1%).

c. Coupling solution simulation calculation example

Due to the influence of temperature on the resistivity of the line, in order to obtain more accurate temperature rise distribution caused by joule heat, analog calculation of electrothermal coupling needs to be carried out. To verify the rationality of the analytical solution and finite volume numerical discretization based methods used in the method, simulation calculations were performed on the circuit structure shown in fig. 13, with the relevant parameters of the circuit listed in table 3. The results of the associated coupled thermal analysis simulation calculations and the results of the modeling analysis of the line in the COMSOL software are presented in fig. 14-16. In which fig. 14 is the total current density distribution, fig. 15 is the joule heat density distribution per unit volume, and fig. 16 is the temperature distribution.

Table 3 circuit parameter table of fig. 13

Base size (mm)	77x 68x 1.2
		Copper foil thickness (mum)	35
Line voltage (V)	0.1V
		Number of discrete units of line	2839
Ambient temperature (. degree. C.)	20
		HTC (W/m) of PCB surface2K)	10
Number of discrete units of COMSOL model	16934

It can be seen that the simulation calculation results according to the proposed coupled thermal analysis method and the modeling analysis results of the COMSOL software show a more consistent distribution trend and a smaller numerical error (the error of the hot spot temperature is about 2.5%). The accuracy of the method and the effectiveness of the iterative operation mode are verified.

The thermal analysis method for the single-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; the calculation of the line current density distribution and the electric-thermal coupled joule heat distribution is realized through numerical discrete approximation of a current continuity equation of the metal layer; the heat diffusion process inside the metal layer is approximately represented by utilizing the linear relation among the temperatures of the metal layer units through carrying out numerical value discrete approximation on the heat balance equation of the metal layer; by taking the thermal resistance parameters of the device into account, the heating rate of the device in the coupling solution is realized, the compensation and the correction of the assumed thermal boundary conditions related to the contact area of the device are realized, the thermal boundary conditions of the whole PCB surface are maintained, and the average temperature of the core of the device and the temperature distribution of other hot spot areas can be finally calculated. Programming realization of the thermal analysis method is carried out in an MATLAB environment, and theoretical accuracy of the method is verified through comparison with COMSOL modeling operation results.

While the embodiments of the invention have been described in detail in connection with the accompanying drawings, it is not intended to limit the scope of the invention. Various modifications and changes may be made by those skilled in the art without inventive step within the scope of the appended claims.

Claims (5)

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

s1, constructing a steady-state heat balance Laplace equation of the single-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, inducing to obtain a metal layer heat diffusion lumped matrix equation, and calculating the heat conduction between the metal layer and a surface heating device;

s4, constructing a thermal analysis coupling equation set of the single-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, 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 PCB surface under the premise of calculating Joule heating of a circuit;

s6, verifying the calculation accuracy of the thermal analysis method of the single-sided PCB structure by adopting a COMSOL modeling operation result;

step S1 is to construct a single-sided PCB insulation layer steady-state heat balance Laplace equation, solve Fourier series analytic solution coefficients based on the boundary conditions of the PCB insulation layer, and obtain an insulation layer temperature analytic solution matrix equation, including:

constructing a steady-state heat balance Laplace equation of the single-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 as compared toTemperature change of ambient temperature or temperature at ambient temperature of 0 ℃, q_iu(x, y) is the heat flow density distribution function introduced by the upper surface, h_uAnd h_dHeat transfer coefficients of the upper and lower surfaces of the insulating layer, k_iIs the thermal conductivity of the insulating layer, L_x、L_y、L_inThe lengths of the insulating layer in the x direction and the y direction and the thickness of the insulating layer in the z direction are respectively measured under a Cartesian coordinate system;

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

wherein, C₁，C₂，C_n,m，

Is a Fourier series coefficient, beta_n、μ_m、γ_n,mCharacteristic values in a Fourier series;

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

the corresponding expression is:

wherein H_uAnd H_dDeriving intermediate parameters of the process for the analytic solution; when h is generated_dIs zero, C₁,

Comprises the following steps:

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

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

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

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

therefore, a single-sided PCB insulating layer temperature analytic solution matrix equation is obtained:

T_u＝R_Mq_M+R_Iq_u

wherein the heat flux density distribution vector of the metal layer transmitted into the insulating layer is q_M(ii) a And q is_uA heat flux density distribution vector, q, transferred between the surface heating component and its covered region_uComprising q_iu(x, y) a heat flow density distribution after discretization; r_MIs q_MA thermal resistance matrix corresponding to a thermal effect generated on the surface of the insulating layer; r_IIs q_uActing on the thermal resistance matrix corresponding to the insulating layer; due to q_uThe heat flux transferred through the surface of the metal layer is taken into account in the thermal diffusion matrix equation corresponding to the metal layer, and the analytical solution matrix equation only takes into account the heat flux transferred between the surface heating element and the surface region of the insulating layer directly covered by the surface heating element, so that R_IWill include a partial all zero column vector.

2. The steady-state thermal analysis method for the single-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 boundary condition of the substrate insulating layer, and calculating the heating rate of the component as the known heat source condition into the matrix equation expression of the analytic solution, including:

let the heat flux density vector of the covering area of the surface heating device be q_h,u：

q_h,u＝q_u+h_uT_J,u

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

the two formulas are combined to obtain:

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

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

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

Wherein R is_IC_P,uq_J,uThe term is a vector of known constants, and T_uIncluding T_J,uAll elements of the vector, and thus the matrix coefficients are transformed by zero-padding up-scaling to achieve T in the matrix equation_uInstead of T_J,uTo achieve the purpose of unifying vectors, the equivalent transformation equation is as follows:

C_JT_u＝R_IC_T,uT_J,u

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

T_u＝R_Mq_M+R_IC_P,uq_J,u-C_J,uT_u。

3. a steady-state thermal analysis method for a single-sided PCB structure as claimed in claim 1, wherein said step S3 is based on finite volume method to approximate the discrete metal layer thermal equilibrium equation, and to generalize the lumped matrix equation for metal layer thermal diffusion, and to take into account the thermal conduction between the metal layer and its surface heat generating device, and comprises:

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 a heat balance equation of the discrete metal layer:

wherein k is_mThe thermal conductivity of the metal layer; q. q.s_cIs the joule heating rate distribution per unit volume of the metal layer; t is_cIs the metal unit temperature; d_cThe unit side length of the structured and discrete metal unit is shown; d_mIs the thickness of the metal layer;q_c,uis T_cThe density of heat flow transmitted from the upper surface of the corresponding unit is the density of heat flow transmitted between the corresponding unit and the heating device; q. q.s_dIs T_cThe heat flux density of the lower surface of the corresponding unit; t is_EIs T_cThe temperature of the adjacent metal unit at the east side of the corresponding unit; t is_WIs T_cThe temperature of the corresponding unit west side adjacent metal unit; t is_NIs T_cThe temperature of the adjacent metal unit at the north side of the corresponding unit; t is_SIs T_cThe temperature of the adjacent metal unit on the south side of the corresponding unit;

and (3) inducing a metal layer thermal diffusion matrix equation based on numerical value dispersion:

M_M,uT_u＝q_MJ,u+C_M,uq_h,u-q_M

wherein q is_MJ,uThe unit joule heating vector corresponding to the metal layer; m_M,uThe coefficient corresponding to the temperature of the metal layer unit in the discrete heat balance equation of the metal layer unit is formed; due to q_h,uThe compensated heat flow directly transmitted through the surface of the insulating layer is already counted in the analytic solution matrix equation, and only the heat flow transmitted between the heating device and the covering metal layer part of the heating device is counted in the lumped matrix equation of heat diffusion, so C_M,uWill also include partial all-zero column vectors;

q is to be_h,uSubstituting the expression into the thermal diffusion matrix equation to obtain:

wherein, C_Mu,uI.e. to realize T in the above equation_J,uAnd T_uThe correlation term is equivalent to the transformed coefficients.

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

wherein, T_u,q_MIf the vector is unknown vector, and other vector and matrix coefficients are known quantities, the unknown temperature vector and the unknown heat flow density vector are obtained by combining the two matrix equations.

5. The steady-state thermal analysis method for the single-sided PCB structure of 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 through calculation and induction, and the potential distribution in the metal layer is obtained, so as to further obtain the current density distribution and the Joule heating distribution of the metal layer, and the Joule heating distribution is included in the metal layer thermal diffusion matrix equation, and iterative calculation is performed between the Joule heating distribution and the temperature distribution, so as to obtain the temperature distribution of the surface of the metal layer under the condition of the Joule heating of the circuit, and the method comprises:

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

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

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

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

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

calculating the conductivity of each cell of the metal layer according to the temperature distribution obtained in the step S4; the conductivity corresponding to the Pc is approximately calculated by using the conductivities corresponding to the four midpoints surrounding the Pc:

for the metal units at the line boundary, the analyzed control volume unit is non-square, and the number 2 control unit surrounded by dark dotted lines has the corresponding approximately discrete current continuity equation:

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

for a 45-degree inclined line in the PCB design process, a square grid is adopted for reconstructing the inclined line, and for a small square control volume unit with a quarter unit area and a Ptr2 vertex at the boundary of the upper left-hand line, the corresponding discrete approximate current continuity equation is as follows:

wherein σ_PmeAnd σ_PmnAre each P_meAnd P_mnThe electrical conductivity of (d); v_Ptr2、V_Pc2And V_Pc1Potentials at Ptr2, Pc2 and Pc1, respectively; the above equation shows that the current continuity equation is the same as 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 are 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 boundary_cA/2, while the length in the actual line is d_cTherefore, such corners need to be individually processed in the analysis method to account for the approximate current continuity equation corresponding to their actual boundaries:

wherein σ_Pmw3、σ_Pmn3And σ_Pms3Are each P_mw3、P_mn3And P_ms3The electrical conductivity of (d); v_Pb3、V_Pb4、V_Pc4And V_Pc5Pb3, Pb4, Pc4 and PThe potential at c 5;

finally, a summary matrix equation of an approximate discrete current continuity equation corresponding to the metal layer is obtained:

M_σV_M＝K_p

wherein, V_MA summary vector of the discrete unit vertex potentials corresponding to all unknown potentials of the metal layer; m_σFor V in the matrix equation_MThe coefficient matrix of (a), is formed by conductivity parameters; k_pKnown terms in a matrix equation are included, including corresponding terms of known endpoint potentials or cross-sectional currents;

after the peak potential distribution in the line is obtained, the potentials corresponding to the centers of the square cells and the centers of the cross sections around the cells are calculated, and the potential V at the center of the cell C2 is obtained by applying a finite volume method to a small diamond region surrounding the center of the cell C2_C2：

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

Wherein, V_C1Potential at the center of cell C1;

thus, the lateral current density and the vertical current density corresponding to the center of each cell are obtained by approximation, and the lateral current density J corresponding to the center of the cell C1_x,C1And longitudinal current density J_y,C1Expressed 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 is shown 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 is obtained₀、a₁And a₂Then, the transverse current density J within the triangle is known_x,ΔAnd longitudinal current density J_y,ΔComprises the following steps:

wherein σ_ΔConductivity corresponding to the triangular region;

on the premise of obtaining the current density corresponding to the metal cell, that is, the joule heating rate per unit volume corresponding to the metal cell, the joule heating rate per unit volume q is obtained for the cell C1_C1Expressed as:

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

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

the joule heating rate is related to the electrical conductivity, which is affected by the temperature, so that when the joule heating of the line is taken into account, an iterative operation between the temperature distribution and the joule heating rate distribution is required, and when the iterative operation is carried outWhen the calculated error is smaller than the set error, the temperature distribution and the Joule heating rate distribution are approximately calculated, and the error is set as the maximum value of the allowable unit temperature change in the two iterative calculations; q. q.s_C1Multiplying the thickness of the metal layer to obtain the heating surface density of the corresponding metal unit, and forming the corresponding Joule heating surface density vector q by the heating surface densities of all the units of the metal layer_MJ,u。

Images