A modeling method for thermal steady-state simulation of the four-layer printed circuit board

Thermal simulation of a Printed Circuit Board (PCB) can help identify potential overheating risks in the circuit. The proposed modeling method combines analytical temperature solutions and numerical approximations. Only Fourier-series analytical solutions related to the prepreg-layer surfaces need to be calculated, rather than the entire structure. Heat transfer through the lateral sides of a PCB is approximately considered as part of the compensated heat flux of the insulating-layer surface boundaries. Heat diffusion within or between metal layers is numerically approximated using the finite volume method. The core layer is treated as “thermally-thick”. Temperature-dependent boundary conditions are considered through iterations. A test solver was developed based on the method. The modeling accuracy was validated by comparison with COMSOL Multiphysics for a four-layer structure with a moderate degree of discretization. Additionally, a PCB for generating DC 3.3V was designed, tested, and modeled, with the modeling results confirmed by the thermal images. The electro-thermal analysis of the distribution of electric potential and Joule heating in traces and vias was integrated into the PCB model. The layout maps of the PCB were further adjusted to reduce Joule heating in the output circuit, and the improvement on reducing the IR drop and hotspot temperature was examined.


Introduction
Thermal simulation of Printed Circuit Boards (PCBs) during the design stage can help optimize layouts [1][2][3][4][5][6][7] and identify potential risks of thermal runaway, ultimately improving circuit reliability [8].Localized overstress may also be found if combined with mechanical analysis [9].One way to perform such simulations is through mathematical modeling based on analyzing the discrete thermal resistances [2,6,10] or discrete effective thermal conductivities [11][12][13] of a PCB, which usually enables the construction of a thermal resistance network for the PCB.A simplified four-layer PCB model was developed to predict SiC MOSFET temperatures for cooling design in an integrated motor drive [10].Only the small PCBs with thermal vias attached to MOSFETS underwent thermal resistance discretization, while the main board was simplified into a single vertical thermal resistance for calculating the equivalent heat transfer coefficient (HTC) [10].This approach facilitated temperature calculations under forced air cooling.However, if lateral thermal conduction in metal layers and other components was considered, the accuracy of the predicted temperatures may be further improved.Considering the trace direction of each discrete region, a continuity check approach was also developed [11].The trace widths at four sub-area boundaries were checked and compared to verify the trace continuity [11].But the typical pattern with two separate bending sub-traces of the same width may be treated as two intersecting traces, leading to an over-estimated effective thermal conductivity.A neural network algorithm has also been employed to predict the effective thermal conductivity map of semiconductor package substrates [12], but the accuracy may depend on the quality of the learning data and the level of discreteness.
Although the computational efficiency of these methods can be higher compared to the detailed modeling using FEM-based method, the accuracy may be compromised because each trace layout pattern may not always be evenly distributed or fully considered in each discrete region.Moreover, in these models, heat transfer through the lateral sides of a PCB or substrate is usually not considered, thus this ignorance is also a potential factor that influences modeling accuracy.
Certainly, several commercial software, such as FloTHERM, Ansys Icepak, COMSOL Multiphysics and Celsius Thermal Solver, can also be utilized to simulate the thermal behaviour of PCBs [1][2][3][4][5][6][7][13][14][15].FloTHERM calculates the regional effective thermal conductivities of a PCB approximately by analyzing the regional cover rates of metal-foil layers [13].But the accuracy of this strategy may heavily rely on the discreteness of the structure.Other software is mainly based on the finite element method (FEM).Results of FEM-based software are usually considered to be agreeable.However, they usually require the complete discretization of the structure, which may reduce computational efficiency [15].
As shown in Fig 1, the proposed modeling method is based on combining the analytical solution of temperature with numerical approximations.Each line indicates a direct relationship between two connected blocks.The brown blocks represent the two main categories of approaches on which the modeling method is based.The green blocks represent the conventional methods involved.The blue blocks represent the calculations specifically related to numerical discretization approximations.The coffee block "Iterations" specifically indicates that when temperature-related thermal boundary conditions are present, the coupled thermal equations will undergo iterative calculations.The benefit of identifying analytical solutions is avoiding complete discretization of the PCB structure, but only the metal layers and their attached insulating surfaces [16,17].The analytical solutions take the form of Fourier series.Other similar analytical approaches have been utilized to establish thermal models of power devices [18][19][20][21], an LED module [22,23] and a heat spreader [24] as well as a cooling mechanism for microprocessor chips [25].
Since there are both metal layers and insulating layers in a PCB, the inhomogeneous structure prevents the use of analytical solutions alone for thermal analysis calculations and necessitates the inclusion of numerical approximations.Numerical approximations involve the compensated heat flux method used to correct assumed unrealistic average HTCs, the finite volume method for approximating heat diffusion and electric potential distribution within metal layers, the multigrid method for enhancing operational efficiency, considering discrete thermal resistance between ICs or components and the board, the "thermallythick" approximation of the core layer, and iterations to address temperature-dependent boundary conditions and calculate Joule heating in the circuit.Most of the approaches are discussed in Section 2, along with the introduction of a test solver developed according to the method.A four-layer structure was modeled using both the test solver and COMSOL Multiphysics to assess the method's feasibility and accuracy.The corresponding comparison is presented in Section 3. A PCB for generating DC 3.3V was also modeled and tested, and the PCB layouts were further improved.The corresponding model and experiment are provided in Section 4.

Fundamentals of the modeling method
As depicted in Fig 2, the four-layer stack consists of two prepreg layers, and one core layer, with each prepreg layer attached to two metal layers.Thus, the structure can also be viewed as two "double-sided PCB" units bonded together through the core.In addition, the temperature distribution of the lateral sides of each thin insulating layer is approximated as the average of the boundary temperature distributions of both its upper and lower surfaces.In this way, heat transfer through the lateral sides of a PCB can thus be considered as one part of the surface compensated-heat-flux.
Since the thermal conductivity of the FR-4 core is significantly lower than that of the attached metal layer, the core layer is assumed as a "thermally-thick" layer.Thus, only vertical heat conduction through the core is considered.A similar approximation method can also be found in COMSOL Multiphysics for modeling a poor thermal conductor with much lower thermal conductivity compared to the adjacent geometry.In this way, the two bonded "double-sided PCB" can be modeled together.

Analytical solutions of the prepreg layer
To model each double-sided structure attached to the core, the analytical solution of the prepreg layer is represented using Fourier series, and heat diffusion in the two metal-foil layers is numerically approximated.These representations are then converted into matrix equations.
Despite the non-uniformity of the prepreg layer caused by the metal vias, their impact on the lateral conduction of heat can be ignored due to the presence of insulation material around them [16].Additionally, vias' contribution to heat transfer between two connected metal layers can be accounted for using the approximation of discrete vertical thermal resistance, which will be explained in Section 2.3.Therefore, the governing equation for steady-state heat diffusion in the prepreg layer under Cartesian coordinates (illustrated in Fig 3) can be approximated by a Laplace's equation: where T denotes the increase in temperature.The corresponding thermal boundary conditions can be assumed firstly as follows: where L x and L y denote the dimensions along X-and Y-axis, respectively; L in denotes the thickness; k i denotes the thermal conductivity of the prepreg layer; q iu (x,y) denotes the vertical heat flux at the upper surface (z = 0) and q id (x,y) denotes that at the lower surface (z = L in ); h u and h d denote the initial assumed average HTCs of the upper and lower surfaces, respectively.In the four-layer PCB, obviously, h d of the top prepreg layer is zero, and h u of the bottom prepreg layer is zero.Consideration of heat transfer through the four lateral sides (x = 0, x = L x , y = 0, and y = L y ) will be explained in the next section.By applying the superposition theorem of thermal effect, Eq (1) can be divided into two sub-equations: where θ and η are sub-variables of T related to q iu (x,y) and q id (x,y), respectively.
To the extent that each group of differential equations is homogenous, their analytical solutions can be expressed as the product of three sub-solutions that are related to only one variable each.For example, the general solution of θ can be expressed as follows: Next, by substituting each group of sub-solutions back into the heat diffusion equation for θ and dividing the corresponding group, the equation is converted to the following form: Since each term of (5) only contains only one independent variable, they must remain constant when any group of (x,y,z) is substituted into the equation: where β n , μ m , and γ n,m are eigenvalues.Taking the boundary conditions into consideration, the relation -β n 2 -μ m 2 +γ n,m 2 = 0 appears to be more appropriate.The Fourier series form of the general solution can then be derived: where C 1 , C 2 , C n,m , and C γn,m are coefficients; n and m are not simultaneously equal to zero in Fourier series.By substituting the general solution into the last thermal boundary condition of θ, C 2 and C γn,m can be derived [17].Then, we can substitute the general solution into the boundary condition of θ at the top surface (z = 0).By further multiplying one set of cosðb n xÞcosðm m yÞ on both sides of the equation each time, and then calculating the integrals along the surface (z = 0), we can obtain C 1 and C n,m [17].This is done based on the integration property of trigonometric functions.
Subsequently, the analytical solution of θ at any position within the prepreg layer can be expressed as an array product that is shown in (8), in which C 1 ' denotes the z-dependent coefficient that includes both C 1 and C 2 ; C qiu,i denotes the (x, y)-dependent coefficient that includes C n,m ; C ' qiu,i denotes the cell coefficient of R u (x,y,z); C n,m,qiu,i denotes the cell coefficient used to derive C ' qiu,i ; q Iu denotes the array form of q iu (x,y); and N e denotes the number of β n or μ m ; and N 1 denotes the number of discrete surface regions; q iu,i denotes the approximate uniform heat flux in the i th individual square region; d q denotes the edge length of each region; and δ m and δ n denote constant coefficients.The analytical solution of η can be obtained in a similar manner.yðx; y; zÞ ¼ C 1 0 ð1 1 � � � 1Þ q iu;1 q iu;2 . . .
Þq Iu ¼ R y;uðx;y;zÞ q Iu R y;uðx;y;zÞ ¼ ðC 0 q iu;1 C 0 q iu;2 � � � C 0 Further, the analytical solutions for both the upper and lower surfaces of each prepreg layer can be expressed as matrix equations: where T u and T d denote the temperature arrays of the upper and lower surfaces, respectively; q Id denotes the array form of q id (x,y); and R u,u , R u,d , R d,u , and R d,d denote the coefficient matrices of thermal resistance related to thermal effects.Therefore, both surfaces of the prepreg layer are divided into discrete surface cells, and in each cell region the heat flux transferred is considered uniform.Particularly, the surface region attached to the metal layer has the same discretization layout as the metal layer to facilitate the consideration of the heat flux transferred from the metal layer.Moreover, the assumed HTCs h u and h d shown in ( 2) and (3), may not accurately represent the actual conditions of some surface regions.For instance, direct thermal convection between covered regions of the PCB's top surface and the environment is almost non-existent due to component presence.Moreover, if radiation heat transfer is taken into account, the equivalent HTCs of some surface regions may greatly deviate from the assumed average value.These unrealistic assumptions can be corrected by taking into account the respective compensated heat flux, which is related to the HTC difference and can be considered based on the superposition theorem of thermal effect.
When considering the compensated heat flux of each discrete region, the expression of the analytical solutions can be further expanded into the following form: where M Hdif denotes the diagonal matrix representing the difference between h u or h d and the actual HTC of each discrete surface cell; q ex denotes the array of compensated heat flux; the subscripts NMu and NMd refer to the upper and lower surface regions not attached to the metal layer, respectively; and the four matrices, R u,NMu , R u,NMd , R d,NMu , and R d,NMd , also denote the coefficient matrices related to thermal effect.For example, R u,NMu q ex,NMu represents the thermal effect of q ex,NMu on T u , and other coefficient matrices related to thermal effect function similarly.The compensated heat flux related to the region attached to the metal layer is going to be considered in the numerical approximation of heat diffusion in the metal layer.
When radiation heat transfer is considered, M Hdif becomes dependent on the temperature distribution of the corresponding surface.In this case, iterations between the temperature distribution and the temperature-dependent HTC of each surface cell must be executed to search for the radiation-equivalent HTC of each cell and eventually identify the solutions.
To incorporate the numerical approximation of heat diffusion in the metal layer, the thermal effects of the metal-attached and non-metal-attached regions must be separately considered.Thus, the expressions of the analytical solution presented in (10) are transformed to the following form.
where, q Mu and q Md denote the heat flux transferred to the upper and lower insulating surface regions attached to the metal layer, respectively; q NMu and q NMd denote the heat flux related to the regions not attached to the metal layer.Thus q Iu is divided into q Mu and q NMu , whereas q Id is divided into q Md and q NMd .Similar to the case of R u,NMu , the four matrices, R u,Mu , R u,Md , R d, Mu , and R d,Md , are also related to thermal effect.

Numerical approximation of heat transfer through the lateral sides of a PCB
Since the thickness of the prepreg layer or core layer is usually less than 2 mm and much smaller than the length and width of a PCB, the lateral sides of the insulating layer can only be discretized horizontally, as depicted in  (2), heat transfer between the lateral sides of a PCB and the environment can also be approximately calculated by multiplying the HTC with the temperature of each individual lateral cell.For example, if a PCB is under the boundary condition of natural convection (still air) with an average HTC of h a , then the array of the heat flux on the lateral sides of the top prepreg layer, q le_tp , can be expressed as follows.
where T le denotes the array consisting of the temperatures of all individual cells on the lateral sides of the layer; the subscript tp refers to the top prepreg layer.
To take q le into consideration, its thermal effect should also be included in the analytical solutions shown in (11).But, depending on the boundary conditions given in (2), only the heat flux at the upper and lower surface of the prepreg layer can be taken into account in the analytical solutions.Thus, q le or T le must be related to the surface heat flux or temperature distributions of the corresponding insulating layer.One option is to average the temperatures of surface boundary cells to approximate those of the lateral cells: where T bc,u and T bc,d denote the temperature arrays of the boundary cells of the upper and lower surfaces of an insulating layer, respectively; M bc,u and M bc,d denote the coefficient matrices that facilitate the transformations between the temperature arrays; T bc denotes the average of T bc,u and T bc,d , and can further be used to derive T le : where T le has four more elements compared to T bc , because of four more corner cells at the lateral sides compared to the surface boundary.M b denotes the coefficient matrix that facilitates the transformation between T le and T bc .M b,u and M b,d also denote the coefficient matrices that facilitate the transformation between these temperature arrays.Thus, (12) can be expressed as follows: If there is no uniform HTC available for the lateral sides, the specified HTC for each lateral cell can also be included in a coefficient matrix like M Hdif to derive q le : Similarly, heat transfer through the lateral sides of the core layer and the bottom prepreg layer can also be approximately described as follows: where the subscripts bp and core refer to the bottom prepreg layer and middle core, respectively.
Further, each product term in the expression of q le can be viewed as part of the compensated heat flux and included in q ex of the corresponding layer-surface.For example, those two terms related to T d_tp can be included in q ex,NMd_tp : where q ex,NMd and M Hdif,NMd maintain their respective meanings, as defined previously for (11); M Hdif,NMd_tp 'T d_tp refers to other compensated heat flux related to T d_tp except the heat flux through the lateral sides.The two terms inside parentheses are subtracted, because q le is considered positive when heat is transferred from the lateral sides of the insulating layer to the environment.Therefore, q le can be approximated as an element of q ex , and its thermal effect can thus be taken into account.

Numerical approximation of heat diffusion in metal-foil layers
In (11), all the heat flux passing through each prepreg-layer surface is depending on both the heat flux transferred to the metal layer and the following process of heat diffusion within the metal layer.A metal layer may comprise multiple metal traces of various patterns.Thus, the aforementioned Fourier-series based analytical solution is unsuitable for modeling the metal layer.Instead, Finite volume method (FVM) [26] can be used to discretize the metal layers and approximate internal heat diffusion.FVM is typically based on simplifying the integration of the differential equation that corresponds to each discrete element [26].The iterative electrothermal modeling method for analyzing IR-drop and Joule heating in metal layers was also partly based on the FVM [16].
As the thickness of the metal layer usually falls in the tens of microns, the vertical temperature variation through the layer can be ignored.Thus, the metal layer can be discretized only in the lateral direction.Consequently, the temperature distribution of the metal layer also mirrors that of the attached surface region of the prepreg layer.Square discrete cells, as shown in Fig 5, can be easily obtained by analyzing the pixel information of the layout map, and are also convenient to generate the multigrid.The governing equation for steady-state heat diffusion of a discrete metal cell has the form of Poisson's equation: where q J denotes the Joule heating rate, and k m denotes the thermal conductivity of the metal layer.Using Gauss's divergence theorem, the volume integration of ( 20) can be converted into closed-surface integration.For example, when considering the central discrete cell depicted in Fig 5, the integration over its six surfaces can be expressed as follows: ð where d c denotes the cell length; A denotes the area variable; d m denotes the thickness; q c denotes the Joule heating rate; and q e , q w , q n , and q s denote the heat flux normal to the four lateral surfaces; q u and q d denote the heat flux transferred into the upper side and out from the lower side of the metal cell, respectively.Thus, corresponding to the first metal layer at the top side of the PCB, q u can represent the compensated heat flux due to the possible non-realistic assumption of h u or the heat flux transferred from the component.By employing Taylor series truncation technique with 2 nd order precision [26], ( 21) can be simplified further to (22).The temperatures of the four neighboring cells depicted in The upper approximate discrete equation of heat diffusion for each group of adjacent metal cells can be gathered in a matrix equation.For example, corresponding to each pair of metal layers attached to a prepreg layer, the approximate matrix equations of heat diffusion can be derived as follows: where the subscripts, Mu and Md , refer to the upper and lower metal layers, respectively; the matrix, M Mu or M Md , is composed of the coefficients of T E , T W , T N , T S , and T c in each discrete equation of the corresponding metal layer; and the matrix, M Vu or M Vd , is composed of the corresponding coefficient of T C,dv , in each discrete equation of the upper or lower metal layer; the terms, M Hdif,Mu T u and M Hdif,Md T d , depend on the compensated heat flux of the metal-layer regions; q c , Mu and q u , Mu include q c d m and q u of each metal cell in the upper metal layer, respectively, whereas q c , Md and q d , Md are defined for the lower layer; and q Mu or q Md includes q d or q u of the corresponding metal layer, thereby representing the heat flux flowing to the attached prepreg layer and determining the analytical solutions as shown previously in (11).
To reduce the number of variables, T Mu and T Md are expressed using T u and T d , respectively: ( where C M,u and C M,d denote the coefficient matrices related to the conversion.Eventually, the matrix equations that describe heat diffusion in two metal layers can be linked to the analytical solution of the attached prepreg layer as follows: Therefore, heat diffusion in metal layers can also be considered as an additional thermal boundary condition of the attached prepreg layer.

"Thermally-thick" approximation of the core layer and the test solver
For easily coupling the two "double-sided PCB" units in a four-layer PCB, the middle core layer is approximated as a "thermally thick" layer.Because there is a significant difference in thermal conductivity between the epoxy-based core layer and the attached copper layers.In a typical FR-4 PCB, the electrodeposited (ED) copper-foil layers of 20~50μm thick have a thermal conductivity of around 316 W/mK [27], much higher than that of the epoxy-based core layer, 0.3 W/mK (according to the material library of COMSOL).Hence, the heat flow through the core can be derived using the surface temperature distributions and discrete thermal resistances of the core.
In this way, the double-sided heat-transfer matrix equations shown in [25] for each group of two metal layers and one prepreg layer can be interconnected, forming the group of equations as depicted in [27] for the four-layer PCB structure.The first four equations correspond to the double-sided structure of the top prepreg layer that is denoted by the subscript tp , and the final four equations with the subscript bp correspond to that of the bottom prepreg layer.The heat flux transferred between two double-sided structures denoted by q u,Md_tp , q NMd_tp , q u, Mu_bp , and q NMu_bp , are approximately related to T u_bp and T d_tp , which also represent the temperature distributions of the corresponding core surfaces, respectively.Further, T u_bp and T d_tp are also used to approximately calculate the temperature distributions of the lateral sides of the core layer.Thus, heat transfer through the lateral sides of the core layer can also be taken into account using the approach of compensated heat flux.C Md_tp , C NMd_tp , C Mu_bp , and C NMu_bp denote the coefficient matrices corresponding to the discrete vertical thermal resistances of the core layer.The remaining arrays, coefficient matrices, and subscripts have been defined previously.
Based on the modeling method, a test solver was developed in MATLAB.The test solver consists of two main parts, as shown in Fig 6 .In the "Model construction program", all the coefficient matrices presented in ( 27) are generated.When modifying a model, the coefficient matrices should only be recalculated if there are changes in the structural parameters, such as the dimensions, or thermal parameters, like the assumed HTCs associated with the analytical solutions.
In the "Model calculation program", the upper equations are solved using a single left-division operation.If there are boundary conditions of temperature-dependent heat flux defined in the model, such as radiation heat transfer, a new group of heat-flux equivalent HTCs can be derived and used for a new round of calculations.Such iterations continue until the preset acceptable errors in both the temperature maps and temperature-dependent HTCs are reached.The total Joule heating in traces and vias can also be iteratively calculated based on the temperature-dependent distribution of electric potential, and is then used to update the power dissipation of components based on the total power dissipation of the circuit (P d ).

Modeling parameters
A four-layer stack-Structure A-as illustrated in Fig 7 was modeled using both the test solver and COMSOL Multiphysics.The stack has dimensions of 29 mm × 40 mm.Each prepreg layer was set to be 0.639 mm thick.The copper thickness of the top and bottom layers was set to 35μm, while that of two middle layers was set to 30μm.This setting was based on the structure parameter of a PCB for generating DC 3.3V, which was also thermally modeled and will be introduced in the next section.The thickness of the core layer, Th core , was set to either 0.15 mm or 0.55 mm during modeling to study its impact on accuracy.
The thermal conductivity of the prepreg material was assumed to be equal to that of the epoxy-based core.This is because "prepreg" only indicates that the epoxy resin has dried but not yet hardened.The stack was assumed to be a common FR-4 type.Thus, as mentioned in the previous section, the thermal conductivity of the insulating layers was taken to be 0.3 W/ mK, and that of the copper foil was taken to be 316 W/mK [27].
There are three copper vias in the structure.The via in the middle only connects the first and second metal-foil layers, whereas the other two vias run through all the layers.The thickness of each via wall was taken to be 18 μm.The heat was assumed to be uniformly transferred into the grey region in the first layer, and P d was set to from 0.2W to 1 W to investigate the stability of the modeling accuracy.The thermal boundary conditions were assumed to be both natural convection (still air) with an equivalent HTC of 10 W/m 2 K [1] and radiation heat transfer.The radiation emissivity of all the surfaces was set to 0.9, an approximate value corresponding to the surfaces of the PCB and ICs [2,30,31].The environment was assumed to be a blackbody at 20˚C.Two modeling cases as shown in Table 1 were compared to investigate both the influence on the accuracy when heat transfer through the lateral sides of the PCB was not included, and that of the numerical approximation of the lateral sides.
Using the test solver, 2D layout maps shown in Fig 7 can be analyzed by first considering each pixel as a discrete cell.The pixel distribution with a resolution of 116 × 160 can be considered as a one-level grid with the unit cell length of 0.25 mm.In order to model the stack more efficiently, a three-level multigrid was used to discretize metal layers and surfaces.Fig 8 shows a schematic diagram of the discrete cells generated under the three-level multigrid.The generation of the multi-level discrete cells can be realized simply by unifying the neighboring cells together [17].To guarantee the accuracy of the hot-spot region, the coarser grids were not used for the heating region.Finally, the number of discrete cells of each layer surface was reduced to less than 5000, thus decreasing the total operation burden significantly: The number of the eigenvalue, β n , is equal to that of μm and is denoted by Ne.Ne was finally set to 300, resulting in about 9 x 10 4 Fourier series in each analytical solution.Further increasing Ne did not have an appreciable impact on the modeling accuracy.During the iterations for calculating radiation heat transfer, the accepted errors of the equivalent HTC and temperature of each discrete cell were set to 0.001W/m 2 K and 0.001˚C, respectively.
In order to verify the modeling results, two 3D models were constructed using COMSOL 5.5, corresponding to the two different Th core values.The mesh view and parameters of the COMSOL model with Th core = 0.15 mm are shown in Fig 9 .The "Heat Transfer in Solids" module and the "Surface-to-Surface Radiation" module were selected to identify the solution.

Modeling results
In Fig 10, temperature maps of the first layer surface from both the Test Solver (TS) and COM-SOL models with Th core = 0.15 mm and P d = 1W are displayed, along with a comparison of temperature profiles along the marking line (x = 2.625mm, y = 13.125mm~28.375mm)in the  heat transfer of the lateral sides in Case-1 produces significant errors.This is because the amount of heat transfer from the lateral sides of the PCB reached about 7.1% and 8.9% of the total heating power for the two structures, respectively.Comparisons under Case-2 with other heating power were also conducted, and the corresponding temperature curves of the marking hot region are depicted in Fig 12 .For models with Th core = 0.15 mm, the average modeling error rates of the marking hot region were around 1.8%, whereas those for Th core = 0.55 mm were around 2.5%.Hence, the "thermally thick" approximation of the core layer didn't exert a significant impact on the accuracy, which can be attributed to the much more effective lateral-heat-diffusion within the copper layers attached to core.
In Fig 13, the comparison of temperature maps of other three layer surfaces are illustrated for the Case-2 model with Th core = 0.15 mm and P d = 1W.Therefore, the modeling results of the test solver were testified to be relatively consistent with those obtained using COMSOL.
The models were run on a computer composed of an Intel 6138 CPU with 40 cores, 256GB memory, and 1TB solid-state drive.Both the preset accepted errors under Case-1 and Case-2 were reached around 10 iterations.The operation time of each TS model was less than 9 minutes, whereas the COMSOL model in each case required around 50 minutes to achieve an approximately converged hotspot temperature with at least 7.6 x 10 5 degrees of freedom.This was based on the results from varying different numbers of mesh elements.The much longer running time of the COMSOL models may be attributed to the combined computation and potential iterations between the "Heat Transfer in Solids" module and the "Surface-to-Surface Radiation" module.Secondly, in the COMSOL models with approximately converged hotspot temperatures, the total number of mesh elements was higher than 5 x 10 5 , whereas that of metal layers was only less than 2 x 10 5 .Hence, probably most of the computing resources were used to determine heat transfer and temperature distribution related to the insulating layers, rather than the metal layers which served as the primary heat diffusion path.

PCB parameters
A four-layer PCB for generating DC 3.3V was also modeled using the test solver.To assess the modeling accuracy, the board was also thermally imaged under natural convection (still air) To achieve a compact circuit, LMZ20501SIL was used as the control IC, in which an inductor is integrated.The initial circuit schematic and PCB layouts were generated using TI's WEBENCH POWER DESIGNER, and the layouts were then redesigned to obtain a smaller size.The copper planes were nearly fully coated for each layer to improve shielding against external interference, create plate capacitors at both the input and output of the power supply, and enhance heat diffusion.
The board is a standard FR-4 type with surface dimensions of 25 mm x 30 mm.According to the manufacturer's specifications, each prepreg layer is approximately 0.639 mm thick, and the core layer is 0.15 mm thick.The top and bottom copper foils are 35μm thick, whereas the other two middle copper layers are 30μm thick.Each via wall is approximately 18 μm thick.

Modeling and measurement
In the model built in the test solver, firstly it was necessary to determine the dissipated power of the components.Through simulation using the WEBENCH POWER DESIGNER, it was determined that with a 5V input and 3.3V/1A output, the control IC U1 resulted in nearly 99.9% of the total heating power loss of all the components.Hence, in the model U1 was approximated as the only heating component.In the iterations for the analysis of temperaturedependent radiation heat transfer, the temperature distributions were thus also used for electro-thermal Joule heating analysis of metal layers.Electrical analysis was separately conducted for the input (V in ) traces, output (V out ) traces and ground (GND) traces based on the connections between layers.
The thermal resistance between a component's junction and each pad cell, R pc , is approximately considered equal to each other.R θJB represents the average thermal resistance between the board and a component's junction, and is more appropriate than R θJC(bot) for indicating the thermal resistance between the component's pads and its junction.This choice was made because R θJC(bot) is measured with the thermal conductive glue filled under the components https://doi.org/10.1371/journal.pone.0310237.g015[28], and using this parameter in the model would lead to an overestimation of heat conduction between the components and board.Therefore, R pc can be approximately derived based on R θJB , and can be used to include heat transfer between the component and the board in the model.Additionally, R θJC(top) can also be used to include heat transfer between the component's surface and the environment [17].The thermal-resistance parameters of ceramic capacitors were derived based on an approximate relation [17] between the equivalent thermal conductivity of the capacitor and the capacitance.The thermal-resistance parameters of those resistors were already specified by the manufacturer [29].Table 2 lists components' thermalresistance parameters, radiation areas and emissivities considered in the model.
The layout maps used for modeling have a resolution of 250 x 300 and are illustrated in Fig 16 .Hence, the first-level discrete cell had a length of 0.1 mm.The emissivities of the PCB surface and IC were set to 0.9 [2,30,31], except for the regions of the vias and soldering, whereas the emissivity of ceramic capacitors was set to 0.94 [31] and that of resistors was set to 0.88 [31].The environment was considered a blackbody.The average HTC of all the sides of the stack was set to 10 W/m 2 K, corresponding to the boundary condition of natural convection (still air) [1].As previously mentioned, the thermal conductivity of the insulating layers was set to 0.3 W/mK, and that of the copper layers was set to 316 W/mK [26].The ambient temperature was set to 10.2˚C based on the measurement.A FLIR T640 thermal camera, as shown in Fig 17, was used to record the temperature.The board was supported by a nylon stick that was about one meter long and 8.6mm in diameter.Heat transfer through the region attached to the nylon stick was neglected.The 5V input was supplied by a programmable DC power supply of RIGOL DP832, and the input current was also measured by the power supply.The output current was sampled by a 100 mΩ (measured resistance: 106 mΩ) resistor and measured by a GDS-2202E oscilloscope.Since there were losses through the lines, the input and output voltages were directly measured at the corresponding terminals of the board.Several cement resistors were used as the load.
Initially, the experiment was conducted with a load of about 1A (Load A).After about half an hour from the start, the board reached thermal equilibrium.The average output was about 0.959A and 3.302V, while the input was approximately 4.902V and 0.69A.Hence, the total dissipating power in the PCB (P td ) was about 0.216W.Fig 19 shows the temperature maps of four layer surfaces obtained from the TS model of the PCB.Based on the thermal parameters of U1, T U1,Top was further derived to be about 28.48˚C, overestimated by about 0.78˚C, so that the simulation error rate of T U1,Top was about 4.46%.Moreover, the total Joule heating loss in the circuit was derived to be about 0.004W, thus the heating rate of U1 was approximately equal to 0.212W, but only about 0.0036W was dissipated from the topside of U1.Because of the iterations for electro-thermal analysis of Joule heating, the operation time of the model was about 2 hours.
The experiment was also conducted with a little heavier load (Load B).In the thermal steady state, the average output was about 1.384A and 3.310V, while the input current and voltage were approximately 1.01A and 4.892V.Hence, P td was about 0.36W.considering the results for both Load A and Load B, the modeling error rate of T U1,Top was between 4% and 5%.Such error rate can be attributed in part to the approximations used in the modeling method and ignorance of thermal conduction through the connecting wires and nylon stick.Furthermore, the actual power loss from the topside of U1 under Load B was estimated to be only about 0.006W, resulting in an estimated U1's junction temperature of about 39.72˚C, only slightly higher than that of the topside.Hence, most heating power of U1 was flowing downwards to the board and dissipated from the PCB surfaces.The temperature maps of other three layer surfaces also indicated thermal diffusion through the PCB, highlighting the key role of the PCB structure in the heat diffusion process.

Improving the layout
The electric-potential maps of the PCB traces with Load B are shown in Fig 22.The Joule heating power in V in /V out /GND traces and the connected vias was estimated to be 8.79 x 10 −4 W/ 7.2x 10 −3 W/8.86 x 10 −4 W, respectively.Hence, the total Joule heating loss in traces and vias was about 0.009W, thus the heating rate of U1 was approximately equal to 0.351W.The IR drop in the traces and vias connected to the input terminals was only 1.1 x 10 −3 V, but that corresponding to the output terminals was about 6.07 x 10 -3 V.As shown in Fig 22b, half of the output drop was attributed to the relatively thin trace connected to the output pin of U1.
Hence, as shown in Fig 23, the layout maps were further modified to reduce Joule heating in the V out traces.In the corresponding model built using the test solver, the heating rate of U1 was set to 0.351W, then iterations including the electro-thermal Joule heating analysis were running in the same way like the previous model.Finally, the IR drop of the V out trace in the top layer has decreased to 1.42 x 10 −3 V as shown in Fig 24, and Joule heating in V out traces and corresponding vias decreased to 1.9 x 10 −3 W, which is more than 70% less.Additionally, as shown in Fig 25, the hotspot temperature of the first layer surface has also decreased to 38.57˚C, about 1.3˚C less.Hence, heat spreading through the traces has also been improved.

Conclusions
Based on the modeling and experiment results shown in previous sections, the proposed method was verified to be feasible.For the model of a four-layer structure, the average modeling error rate using the test solver was less than 3% based on comparison with COMSOL Multiphysics, while the corresponding running time was less than 9 minutes, much lower than that of the COMSOL models.Hence, utilizing the analytical solution permits a moderate degree of discretization.For the PCB model, the modeling error rate of the surface temperature of the hottest IC was between 4% and 5%.Therefore, both the "thermally-thick" approximation of the core layer and the numerical approximation of heat transfer through the lateral sides of the PCB were testified to be reasonable.It can be predicted that the method may also be suitable for other multilayer PCB structures with six or more metal-foil layers, as well as for some semiconductor power modules with a laminate structure comprising multiple layers.

Fig 5 .
Fig 5. Schematic diagram of heat conduction through discrete cells in the metal layer.https://doi.org/10.1371/journal.pone.0310237.g005 Fig 5 are denoted by T E , T W , T N , and T S , respectively.T c denotes the temperature of the central cell.ψ Nv, z denotes the discrete vertical thermal resistance of the via in the insulating layer, and T c,dv denotes the temperature of the connected via cell in the other metal layer that is attached to the same insulating layer.

Fig 22 .Fig 23 .
Fig 22. Simulated electric-potential maps of the main-circuit traces of the PCB with Load B. (a) V in (V) trace of the top layer.(b) V out (V) trace of the top layer.(c) V GND (10 -4 V) trace of the top layer.(d) V GND (10 -4 V) trace of the 2 nd layer.(e) V in (V) trace of the 3 rd layer.(f) V out (V) trace of the 3 rd layer.(g) V GND (10 -4 V) trace of the 4 th layer.https://doi.org/10.1371/journal.pone.0310237.g022

Fig 24 .
Fig 24.Simulated electric-potential (V) maps of the V out traces of the new layout of the PCB with Load B. (a) the top layer.(b) the 3 rd layer.https://doi.org/10.1371/journal.pone.0310237.g024

Table 1 . Modeling condition definition related to the model of Structure A. Modeling case in the Test Solver Case-1 Case-2
Heat transfer through the lateral sides of the PCB Not included included https://doi.org/10.1371/journal.pone.0310237.t001