Circuit board heat dissipation structure optimization design method based on multi-physical field coupling

Abstract

The application discloses a circuit board heat dissipation structure optimization design method based on multi-physical field coupling and concretely relates to the field of circuit board heat dissipation design; after obtaining lamination parameters, power consumption time sequence, installation constraints and environmental boundaries, an electric-thermal-flow-force multi-physical field coupling twin simulation model is established; micron-level gap grids are generated between the copper layer and the dielectric layer at the step junction, the gas heat conduction contribution is set to zero when the air pressure is lower than 1 kPa, and the near-field radiation heat transfer thermal conductivity is calculated according to the gap width and the surface emissivity, which is loaded to the thermal field equation as an additional interfacial thermal resistance; the local failure risk area is determined based on the temperature field, the flow field and the deformation stress field obtained through simulation, the via density, the heat block size, the interface layer thickness and the heat dissipation fin size are extracted for iterative optimization, and the optimal variable is output to select a heat dissipation structure optimization scheme; the application improves the accuracy of the circuit board thermal simulation and the reliability of the heat dissipation structure optimization under a vacuum environment.

Description

Technical Field

[0001] This invention relates to the field of circuit board heat dissipation design, and more specifically to a method for optimizing the design of circuit board heat dissipation structures based on multi-physics coupling. Background Technology

[0002] As high-power-density electronic devices develop towards miniaturization, sealing, and complex installation boundaries, sandwich stepped circuit boards, which can take into account partitioning, structural avoidance, and local load-bearing, have been gradually applied to airborne radio frequency components, spaceborne control units, and high-reliability power modules.

[0003] However, in actual operation, the copper layer distribution, dielectric thickness, thermal via density, and interface contact state of such circuit boards vary significantly in different stepped regions. Especially under the combined effects of device pulse power consumption timing changes, confined flow channel heat dissipation, and mounting pre-tightening, local heat accumulation and stress coupling are easily formed at stepped corners, adjacent areas of via groups, and transition connections of heat sinks. This can lead to problems such as slight warping, interface delamination, abnormal local voltage drops, and hot spot migration. Existing technologies mostly adopt heat dissipation structure design methods based on empirical parameters, or only perform single simulation analysis on temperature field and flow field. It is difficult to uniformly correlate device power consumption timing, environmental boundaries, structural deformation, and flow state. It is also difficult to accurately identify local failure risk areas during the design stage and further complete the coordinated optimization of via density, heat spreader size, interface layer thickness, and heat sink fin size. Summary of the Invention

[0004] The purpose of this invention is to provide an optimized design method for heat dissipation structure of circuit boards based on multi-physics coupling, so as to overcome the shortcomings of the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a circuit board heat dissipation structure optimization design method based on multi-physics field coupling, comprising: Obtain the stack-up parameters, device power consumption timing, installation constraints, and environmental boundaries of the target sandwich stepped circuit board to obtain the initial parameter set P; Based on the initial parameter set P, establish an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W, which includes a copper layer, a dielectric layer, thermally conductive vias, interface materials, and heat dissipation components; Based on the power consumption timing and environmental boundary in the initial parameter set P, the twin simulation is solved using the multiphysics coupled twin simulation model W to obtain the corresponding temperature field T, flow field V and deformation stress field M. Based on the temperature field T, flow field V, and deformation stress field M, the local failure risk zone R at the step corner is determined, and the through hole density, heat spreader size, interface layer thickness, and heat dissipation fin size corresponding to the risk zone R are extracted as variables to be optimized X. The variable X to be optimized is mapped to a multiphysics coupled twin simulation model W for iterative twin simulation to obtain the objective function Y composed of peak temperature, local temperature difference, warpage and flow pressure drop; Based on the objective function Y, the optimal variable to be optimized is selected under preset constraints. The output is the same as the optimal variable to be optimized. Corresponding circuit board heat dissipation structure optimization scheme.

[0006] Preferably, the specific methods for establishing a coupled electro-thermal-fluid-mechanical multiphysics twin simulation model W include: A micron-scale gap grid is generated between the copper layer and the dielectric layer at the step junction, and the actual width of each gap unit is obtained by optical profile measurement. When the environmental boundary in the initial parameter set P contains a gas pressure below 1 kPa, the gas thermal conductivity contribution in the gap element is set to zero, and the near-field radiation heat transfer thermal conductivity is calculated based on the gap width and the emissivity of the copper layer surface. The near-field radiation heat transfer thermal conductivity is used as an additional interface thermal resistance boundary condition at the interface between the copper layer and the dielectric layer, and is loaded into the thermal field equation of the electro-thermal-fluid-mechanical multiphysics coupled simulation to form an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W.

[0007] Preferably, the specific methods for solving twin simulations using the electro-thermal-fluid-force multiphysics coupled twin simulation model W include: The power consumption timing is divided into continuous time steps, and the normal contact pressure at the step corner in the deformation stress field M is read at the beginning of each time step. The actual contact gap width at the interface between the copper layer and the dielectric layer is updated based on the normal contact pressure, and the near-field radiative heat transfer thermal conductivity is recalculated based on the updated gap width. The recalculated thermal conductivity is substituted into the thermal field equation to solve for the temperature field T and flow field V at the current time step, and then the deformation stress field M at the next time step is updated based on the temperature field T.

[0008] Preferably, updating the actual contact gap width at the interface between the copper layer and the dielectric layer based on the normal contact pressure includes: Extract the roughness amplitude and micro-protrusion radius of curvature of the copper layer surface and the dielectric layer surface from the initial parameter set P; Substituting the normal contact pressure into the Greenwood-Williamson contact model, the actual contact area ratio and the average flattening amount of the micro-protrusions were calculated. Subtract the average flattening amount of the micro-protrusions from the nominal gap width, and then multiply by the proportion of the non-contact area to obtain the updated actual contact gap width.

[0009] Preferably, the specific methods for determining the local failure risk zone R include: The temperature gradient amplitude of each step angle is extracted from the temperature field T, and the equivalent stress amplitude of each step angle is extracted from the deformation stress field M. The coupling risk index is obtained by multiplying the temperature gradient magnitude by the equivalent stress magnitude, and the rate of change of flow pressure drop near the step angle is extracted from the flow field V. When the coupling risk index exceeds the preset threshold and the rate of change of flow pressure drop is positive, the step corner is marked as a local failure risk zone R.

[0010] Preferably, the temperature gradient magnitude at each step turning angle is extracted from the temperature field T, including: At each step corner, select two adjacent nodes along the copper layer extension direction, the dielectric layer normal direction, and the bisector of the angle between the two, and record the temperature difference and straight-line distance between the three sets of nodes. Divide each temperature difference by the corresponding straight-line distance to obtain the component gradient values ​​in three directions. Then, square the three component gradient values ​​respectively and sum them to obtain the temperature gradient amplitude.

[0011] Preferably, the variable X to be optimized is mapped to an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W for iterative twin simulation, including: The density of through holes, the size of the heat spreader, the thickness of the interface layer and the size of the heat dissipation fins are treated as independent design variables, and the initial variable combination is generated by Latin hypercube sampling within their respective allowable ranges. Each set of variables is written into the corresponding position of the electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W by parameter substitution, and transient simulation is performed sequentially to obtain peak temperature, local temperature difference, warpage and flow pressure drop; With the goal of minimizing peak temperature and ensuring that local temperature difference, warpage, and flow pressure drop do not exceed their respective upper limits, radial basis function interpolation is used to construct the response relationship between variables and the target. After optimization on the response surface, the variable combination is updated for the next iteration.

[0012] Preferably, the optimal variable to be optimized is selected based on the objective function Y under preset constraints. ,include: Subtract the preset upper limit of the local temperature difference from the measured value of the local temperature difference in the objective function Y. If the difference is positive, divide the difference by the preset upper limit to obtain the local temperature difference excess ratio. Similarly, calculate the warpage excess ratio and the flow pressure drop excess ratio. If the difference is negative or zero, the excess ratio is zero. The weighted peak temperature is obtained by weighting and summing the measured peak temperature, the local temperature difference excess ratio, the warpage excess ratio, and the flow pressure drop excess ratio. The combination of variables with the smallest weighted peak temperature is selected as the optimal variable to be optimized. .

[0013] The technical effects and advantages provided by the present invention in the above technical solution are as follows: 1. This invention overcomes the fundamental shortcomings of existing empirical contact thermal resistance models at ambient pressure under vacuum or low-pressure environments. During the construction of the multiphysics coupled twin simulation model W, a micrometer-scale gap mesh is explicitly generated between the copper layer and the dielectric layer at the stepped interface. The gas thermal conductivity contribution is set to zero based on ambient pressures below 1000 Pa or 100 Pa. Simultaneously, near-field radiative heat transfer thermal conductivity based on gap width, surface emissivity, and dielectric function is introduced as an additional interface thermal resistance boundary condition, replacing the original empirical contact thermal resistance at ambient pressure. Furthermore, in each transient simulation time step, the normal contact pressure fed back from the deformation stress field M is read, and the actual contact gap width is dynamically updated using the Greenwood-Williamson contact model. Based on this, the near-field radiative heat transfer thermal conductivity is recalculated and then substituted into the thermal field equations for solution. These methods enable the twin simulation model W to realistically reflect the microscopic heat transfer mechanism in the submicron-level "virtual contact" gap under vacuum conditions, solving the technical problem that the temperature prediction is seriously understated due to the neglect of near-field radiation in existing technologies, thus providing a high-fidelity simulation basis for the optimization of heat dissipation structures.

[0014] 2. This invention achieves accurate identification and multivariate collaborative optimization of local failure risk zones based on the coupling results of multiple physics fields (thermal-fluid-mechanical). After obtaining the temperature field T, flow field V, and deformation stress field M through twin simulation, the coupling risk index is obtained by multiplying the temperature gradient amplitude at the stepped corner with the equivalent stress amplitude. Combined with the flow pressure drop rate in the flow field, the local failure risk zone R is jointly determined, avoiding the limitations of single-physics field analysis. Then, the through-hole density, heat spreader size, interface layer thickness, and heat dissipation fin size corresponding to the risk zone R are used as variables to be optimized, and iterative simulation is performed in the twin simulation model W. The objective function Y is constructed with peak temperature, local temperature difference, warpage, and flow pressure drop, and the optimal variable combination is selected by minimizing the weighted peak temperature. This process unifies power consumption timing, structural deformation, flow state, and thermal effects, enabling coordinated control of heat accumulation, warpage deformation, and flow resistance during the design phase. It effectively prevents failure modes such as interface delamination, hot spot migration, and abnormal voltage drop at stepped corners. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0016] Figure 1 This is a flowchart of the circuit board heat dissipation structure optimization design method based on multi-physics field coupling of the present invention.

[0017] Figure 2 This is a flowchart of the method for determining the local failure risk zone of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] For examples, please refer to Figure 1 As shown in this embodiment, the circuit board heat dissipation structure optimization design method based on multi-physics coupling includes: The stack-up parameters, device power consumption timings, mounting constraints, and environmental boundaries of the target sandwich-type stepped circuit board are obtained to obtain an initial parameter set P. The initial parameter set P includes four main categories of parameters: stack-up parameters, device power consumption timings, mounting constraints, and environmental boundaries. The specific acquisition method is as follows: The target sandwich stepped circuit board refers to a printed circuit board containing at least two stepped regions of different thicknesses, and internally having a conductive copper layer, an insulating dielectric layer, thermally conductive vias, and a sandwich layer (optional). Layer stack-up parameters are obtained through the following methods: from circuit board design drawings (such as Gerber files, PCB design source files), process stack-up tables, bills of materials, or direct measurement. Specifically, this includes: Copper layer parameters: copper layer thickness in each step area (e.g., 17μm, 35μm, 70μm), number of copper layers, copper layer distribution location (top layer, bottom layer, inner layer), and copper coverage.

[0020] Dielectric layer parameters: material type (e.g., FR4, polyimide, ceramic-filled composite material), thickness, thermal conductivity, density, specific heat capacity, Young's modulus, Poisson's ratio, and coefficient of thermal expansion of each dielectric layer.

[0021] Sandwich layer parameters: If the circuit board contains a metal core (such as aluminum or copper) or a ceramic core, obtain the sandwich layer material, thickness, location (which step area it is located in), and its connection method with the copper layer and dielectric layer.

[0022] Thermal via parameters: diameter of vias in each area (e.g., 0.2mm to 0.5mm), copper thickness of via wall, filling material (resin, copper paste or thermally conductive adhesive), distribution density of vias on the plane (number per unit area), and thermal layer connected to the via.

[0023] Interface material parameters: the initial thickness, compressive modulus, thermal conductivity and hardness of the thermal interface materials (such as thermal pads, thermal grease, phase change materials, thermal gels) used between the chip and the circuit board, and between the circuit board and the heat sink.

[0024] Heat sink parameters: geometric dimensions, material properties, and installation location of the heat sink (finned or pin-finned), heat spreader, heat pipe, or liquid cooling plate.

[0025] Device power consumption timing data describes the relationship between the power consumption of various heat-generating devices (such as power amplifiers, processors, power management chips, and RF transceiver components) mounted on a circuit board and the time required for operation. Acquisition methods include: consulting typical power consumption curves in the device datasheet, extracting transient power waveforms through circuit simulation (such as SPICE), or measuring the power consumption of an actual prototype under typical operating conditions. For pulsed power consumption, the pulse width, duty cycle, rise / fall time, and peak power must be recorded; for periodically varying power consumption, the power sequence within a complete cycle must be recorded. Finally, the power consumption timing data of all devices is stored in time series format (e.g., one sampling point per second or millisecond) and associated with the device's coordinate position on the circuit board.

[0026] Mounting constraints refer to the mechanical fixing conditions of the target circuit board in the whole machine or module, including: Installation methods: screw fixing, clip clamping, guide rail sliding, or welding fixing.

[0027] Preload / Preload torque: The range of preload at each fixing point (e.g., 2N to 5N) or the screw tightening torque (e.g., 0.3Nm to 0.6Nm).

[0028] Fixed point location and number: Mark the coordinates of all mounting points in the circuit board coordinate system, as well as the corresponding constraint degrees of freedom (such as full constraint, normal constraint only).

[0029] Mounting clearance and tolerance: The nominal clearance and manufacturing tolerance between the circuit board and the housing or supporting structure.

[0030] Adjacent component constraints: Whether there are other components (such as connectors, wiring harnesses, shielding covers) in the vicinity that physically interfere with the size of the heat dissipation structure.

[0031] The above constraints can be obtained through the machine's mechanical drawings, assembly process documents, or actual measurements.

[0032] The environmental boundary defines the thermal and fluid environment conditions in which the circuit board operates, including: Atmospheric pressure: normal pressure (101325Pa), low pressure (such as 1000Pa, 100Pa) or vacuum (<1Pa). For aerospace applications, the orbital vacuum level needs to be specified.

[0033] Ambient temperature: The temperature of the surrounding air or the casing (which can be a constant value or a temperature curve over time).

[0034] Cooling medium and flow conditions: natural convection, forced air cooling (wind speed, wind direction, turbulence intensity), liquid cooling (coolant type, flow rate, inlet temperature).

[0035] Radiation boundary: emissivity, temperature, and viewing angle coefficient of the surrounding wall surface.

[0036] Gravity direction: Determine the direction of the gravitational acceleration vector in order to calculate the natural convection buoyancy force.

[0037] The above environmental boundary parameters are determined based on the design specifications or measured data of the target application scenario (such as ground cabinets, high-altitude aircraft, and spaceborne equipment).

[0038] All the above parameters are organized according to a unified naming rule and data structure to form an initial parameter set P.

[0039] Based on the initial parameter set P, an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W is established, which includes a copper layer, a dielectric layer, thermally conductive vias, interface materials, and heat dissipation components.

[0040] Based on the stacking parameters in the initial parameter set P, the positions of all stepped corners on the target sandwich stepped circuit board are determined. At each stepped corner, the interface between the copper layer and the dielectric layer is selected as the key area of ​​interest. Using a 3D solid mesh generation tool, thin mesh cells with a thickness scale of micrometers are generated in the normal direction of this interface, called gap cells.

[0041] The planar dimensions of the gap cells (i.e., the length and width along the interface) are set to 10 micrometers by 10 micrometers to adequately resolve submicrometer-level contact gap variations. The gap cells consist of at least two layers, corresponding to the air or vacuum region near the copper layer surface and the region near the dielectric layer surface, respectively. After mesh generation, a white light interferometric profilometer is used to measure the copper layer surface and dielectric layer surface at the corresponding step corners of the actual circuit board sample.

[0042] During measurement, 10 different locations were selected, and the arithmetic mean of the surface height data at each location was recorded as the surface profile at that location. The height difference between the copper layer surface and the dielectric layer surface in the same plane coordinate system was defined as the actual width of the gap element. For each gap element, its actual width value was recorded in nanometers. This actual width will be used to subsequently determine whether the gas thermal conductivity contribution is negligible and to calculate near-field radiative heat transfer.

[0043] When the environmental boundary in the initial parameter set P contains a gas pressure below 1000 Pa, the gas thermal conductivity contribution within the gap unit is set to zero, and the near-field radiative heat transfer thermal conductivity is calculated based on the gap width and the emissivity of the copper layer surface.

[0044] The air pressure value is read from the environmental boundary of the initial parameter set P. If the air pressure is lower than 1000 Pa (i.e., one-thousandth of a standard atmosphere), it is determined that the mean free path of the gas molecules in the gap unit is much greater than the gap width, and the gas heat conduction mechanism fails.

[0045] Therefore, in the subsequent thermal field calculations, the heat conduction contribution of the gas in all gap elements is set to zero, that is, the heat transfer by air or rarefied gas through conduction is not considered.

[0046] Calculate the near-field radiative heat transfer thermal conductivity within the gap element. Obtain the surface emissivity of the copper layer material, which is derived from the physical properties of the copper layer material in the stacking parameters of the initial parameter set P. For standard electrolytic copper foil, its surface emissivity ranges from 0.03 to 0.10 at room temperature, with the specific value selected within the range of 0.05 to 0.30 depending on the actual surface treatment process (e.g., whether gold plating or passivation is used). Calculate the near-field radiative heat transfer theoretically based on the actual width of the gap element. Multiply Planck's constant by the speed of light and divide by twice pi to obtain the first intermediate value; multiply Boltzmann's constant by the absolute temperature and divide by Planck's constant to obtain the second intermediate value; multiply the gap width by the second intermediate value and divide by the speed of light to obtain a dimensionless gap parameter. When this dimensionless gap parameter is less than 0.1, the near-field radiative heat transfer thermal conductivity is calculated as follows: multiply the Stefan-Boltzmann constant by the cube of the absolute temperature, then multiply by the surface emissivity of the copper layer, and then multiply by a correction factor inversely proportional to the gap width. The value of the correction factor is equal to the second intermediate quantity multiplied by the imaginary part of the dielectric function of the copper layer material at the characteristic frequency, and then divided by the gap width.

[0047] For different gap widths, repeat the above calculation to obtain the corresponding near-field radiation heat transfer thermal conductivity value, and store this value in association with the actual width of the gap element.

[0048] Near-field radiative heat transfer thermal conductivity is used as an additional interfacial thermal resistance boundary condition at the interface between the copper layer and the dielectric layer, and is applied to the thermal field equation of the electro-thermal-fluid-mechanical multiphysics coupled simulation to form an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W. The specific implementation process is as follows.

[0049] The thermal conductivity of each gap element is converted into an interfacial heat transfer coefficient. For each gap element, the calculated near-field radiative heat transfer thermal conductivity is directly used as the interfacial heat transfer coefficient for that gap element. This coefficient represents the heat flux density transferred per unit area of ​​the gap between the copper layer and the dielectric layer under a unit temperature difference. When the ambient pressure is below 1000 Pa, the interfacial heat transfer coefficient corresponding to the empirical contact thermal resistance at normal pressure is completely discarded and no longer included in subsequent calculations.

[0050] In practice, in the material properties or boundary condition definition area of ​​the simulation software, locate the contact thermal resistance setting between the copper layer surface and the dielectric layer surface, and replace the original fixed heat transfer coefficient (e.g., 20,000 W / m² / Kelvin at atmospheric pressure) with the near-field radiative heat transfer thermal conductivity value calculated above. If different gap elements have different gap widths, each gap element should be set with its own independent interface heat transfer coefficient.

[0051] The thermal field equations are applied as heat flux boundary conditions. The thermal field equations are based on the transient Fourier heat conduction equation. At each gap element after finite element discretization, a heat flux boundary condition is applied between the copper layer surface nodes and the dielectric layer surface nodes. The thermal field equations are calculated as follows: the heat flux density equals the interfacial heat transfer coefficient multiplied by the absolute value of the difference between the copper layer surface temperature and the dielectric layer surface temperature, with the direction pointing from the high-temperature side to the low-temperature side. This heat flux boundary condition is manifested in the thermal field equations as follows: for the copper layer side, this heat flux appears as a negative heat load (i.e., a heat outflow term) on the right-hand side of the equation; for the dielectric layer side, this heat flux appears as a positive heat load (i.e., a heat inflow term) on the right-hand side of the equation.

[0052] In practice, the temperature values ​​of the nodes on both sides of the current gap element are read in each iteration of each time step through the user subroutine or custom boundary condition interface of the simulation software. The temperature difference is calculated, multiplied by the interface heat transfer coefficient corresponding to the gap element, and the current heat flux density is obtained. Then, the heat flux density is distributed to the heat load vector of the nodes on both sides.

[0053] After setting the heat flux boundary conditions, the thermal field equation is coupled with the flow field equation and the force field equation for solution. The flow field equation adopts the Navier-Stokes equations, applying slip boundary conditions at the wall for low-pressure environments; the force field equation adopts linear elastic or elastoplastic constitutive equations, considering thermal strain caused by temperature changes and mechanical strain generated by installation preload. The three physical fields are coupled in the following way: the internal heat source term in the thermal field equation is provided by the device power consumption timing; the temperature field T obtained from the thermal field solution is used as a thermal load input into the force field equation to calculate thermal stress and thermal deformation; the normal contact pressure in the deformation stress field M obtained from the force field solution is fed back to the gap element to update the actual contact gap width, thereby updating the near-field radiative heat transfer thermal conductivity, and then fed back to the thermal field equation; the flow pressure drop and velocity distribution obtained from the flow field solution affect the convective heat transfer coefficient, which is applied as a thermal boundary condition to the thermal field equation. The aforementioned coupling relationship is achieved by performing sequential iterations at each time step until the residuals of all physical fields are less than the convergence tolerance (e.g., temperature residual less than 0.01 degrees Celsius, velocity residual less than 0.001 meters per second, displacement residual less than 0.001 micrometers). After completing all settings, an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W is obtained, capable of simultaneously solving the temperature field T, the flow field V, and the deformation stress field M.

[0054] By using the above method, the near-field radiation heat transfer mechanism is accurately embedded into the multiphysics simulation framework, so that model W can truly reflect the near-field radiation heat transfer contribution generated by the submicron-level virtual contact gap between the copper layer and the dielectric layer under vacuum or low pressure environment, thereby overcoming the fundamental defects of the existing atmospheric pressure empirical contact thermal resistance model in this type of environment.

[0055] Based on the power consumption timing and environmental boundaries in the initial parameter set P, a twin simulation is performed using a multiphysics coupled twin simulation model W to obtain the corresponding temperature field T, flow field V, and deformation stress field M.

[0056] Read the device power consumption timing data from the initial parameter set P. This timing data records the power consumption values ​​of each heat-generating device at different times, in watts. Divide the entire operating time interval into several consecutive time steps. The choice of time step size depends on the drasticness of power consumption changes: for pulsed power consumption, the time step size is one-tenth to one-twentieth of the pulse width; for example, if the pulse width is 2 milliseconds, the time step size is 0.1 milliseconds to 0.2 milliseconds. For slowly changing power consumption, the time step size can be 0.5 seconds to 1 second.

[0057] At the start of each time step (i.e., the start of the nth time step), the currently solved deformation stress field M in the twin simulation model W is invoked (initially, the deformation stress field M is obtained from the static mechanical solution under the initial installation constraints). The normal contact pressure at the nodes located within each step corner region is extracted from the deformation stress field M. The normal contact pressure refers to the pressure exerted by the interaction between the copper layer surface and the dielectric layer surface along the normal direction at the contact interface.

[0058] For each gap cell within the stepped corner region, the contact pressure value calculated using the mechanical equilibrium equations is read from both the copper and dielectric nodes. For each stepped corner, the arithmetic mean of the normal contact pressures on all gap cells within that region is taken as the representative normal contact pressure for that corner at the current time step. This pressure value will be used to subsequently update the actual contact gap width.

[0059] Before proceeding with the contact gap update calculation at each time step, morphological parameters for the copper and dielectric surfaces are extracted from the stacking parameters of the initial parameter set P. These parameters include two items: roughness amplitude (root mean square roughness) in nanometers and micro-protrusion radius of curvature in micrometers. The roughness amplitude reflects the degree of vertical deviation of the surface from the average plane, while the micro-protrusion radius of curvature reflects the radius of curvature at the tips of microscopic protrusions on the surface. For the copper surface, if not directly provided in the initial parameter set, default values ​​are used based on the copper foil manufacturing process: for rolled copper foil, the roughness amplitude is 0.3 to 0.5 micrometers, and the micro-protrusion radius of curvature is 5 to 10 micrometers; for electrolytic copper foil, the roughness amplitude is 0.8 to 1.5 micrometers, and the micro-protrusion radius of curvature is 2 to 5 micrometers. For the dielectric surface (such as FR4 or polyimide), the roughness amplitude is typically 0.1 to 0.3 micrometers, and the micro-protrusion radius of curvature is 10 to 20 micrometers.

[0060] In this embodiment, the extracted surface roughness amplitude of the copper layer is 0.4 micrometers, and the radius of curvature of the micro-protrusion is 8 micrometers; the surface roughness amplitude of the dielectric layer is 0.2 micrometers, and the radius of curvature of the micro-protrusion is 15 micrometers. To simplify the calculation, the combined roughness amplitude of the two surfaces is defined as the square root of the sum of the squares of the roughness amplitudes of the copper layer and the dielectric layer, and the combined micro-protrusion radius of curvature is defined as the reciprocal of the sum of the reciprocals of the radii of curvature of the copper layer and the dielectric layer. The calculated combined roughness amplitude is approximately 0.447 micrometers, and the combined micro-protrusion radius of curvature is approximately 5.22 micrometers.

[0061] Substituting the normal contact pressure into the Greenwood-Williamson contact model, the actual contact area ratio and the average flattening amount of the micro-protrusions were calculated.

[0062] The Greenwood-Williamson contact model is a classical mechanical model used to describe the elastic deformation of micro-protrusions under contact pressure between two rough surfaces. This step uses this model to calculate two key parameters: the actual contact area ratio and the average flattening of the micro-protrusions. The actual contact area ratio is defined as the proportion of the area of ​​the micro-protrusions that actually engages in solid-state contact within the nominal contact area; it is dimensionless. The average flattening of the micro-protrusions is defined as the average compressive deformation of all contacting micro-protrusions at their contact point; it is measured in nanometers.

[0063] First, a dimensionless parameter called the Plass exponent is calculated based on the combined roughness magnitude (denoted as sigma) and the combined micro-protrusion radius of curvature (denoted as beta). The Plass exponent is calculated by multiplying the combined roughness magnitude by the square root of the combined micro-protrusion radius of curvature, and then multiplying by a constant (approximately 0.0004), but the specific value depends on the material's elastic modulus and Poisson's ratio. A more precise implementation involves calculating the combined elastic modulus using the elastic modulus of the copper layer and the dielectric layer (read from the initial parameter set P; the elastic modulus of copper is 110 GPa, and that of the FR4 dielectric layer in the plate direction is 18 GPa) and Poisson's ratio (0.34 for copper and 0.15 for FR4). The formula is: ;in, The combined elastic modulus is used to describe the combined elastic deformation characteristics of two contacting surfaces. This represents the elastic modulus of the copper layer material, with a typical value of 110 GPa. is the Poisson's ratio of the copper layer material, dimensionless, with a typical value of 0.34. This represents the elastic modulus of the dielectric layer material, with a typical value of 18 GPa in the plate direction for FR4 material. is the Poisson's ratio of the dielectric layer material, dimensionless, with a typical value of 0.15.

[0064] Multiply the normal contact pressure (denoted as Pc, in Pascals) by the radius of curvature of the combined micro-convexity, and then divide by the combined elastic modulus to obtain a dimensionless pressure parameter. Multiply this dimensionless pressure parameter by an empirical coefficient (usually 2.5), and then obtain the actual contact area ratio by looking up a table or using a polynomial approximation function.

[0065] For example, when the normal contact pressure is 1 MPa, the calculated dimensionless pressure parameter is approximately (1 x 10⁶ Pa x 5.22 x 10⁻⁶ m) divided by (17.6 x 10⁹ Pa), which is approximately 2.97 x 10⁻¹⁰. Multiplying this by 2.5 gives 7.42 x 10⁻¹⁰. The actual contact area ratio corresponding to this value is very small, approximately between 0.001 and 0.01. For this embodiment, the classic approximation formula of the Greenwood-Williamson model is used: the actual contact area ratio equals the normal contact pressure multiplied by the radius of curvature of the combined micro-protrusions, divided by (the combined elastic modulus multiplied by the combined roughness amplitude), and then multiplied by a constant (approximately 0.4). The calculated actual contact area ratio is: 1 x 10^6 Pascals x 5.22 x 10^-6 meters divided by (17.6 x 10^9 Pascals x 0.447 x 10^-6 meters) multiplied by 0.4, approximately equal to 0.00265. The formula for calculating the average flattening amount of the micro-protrusion is: multiply the actual contact area ratio by the combined roughness amplitude, then multiply by a correction factor (usually between 0.8 and 1.2). In this embodiment, the correction factor is 1.0, resulting in an average flattening amount of 0.00265 multiplied by 0.447 micrometers, approximately equal to 0.00118 micrometers, or 1.18 nanometers.

[0066] Subtract the average flattening amount of the micro-protrusions from the nominal gap width, and then multiply by the proportion of the non-contact area to obtain the updated actual contact gap width. Based on the updated gap width, the near-field radiative heat transfer thermal conductivity is recalculated.

[0067] The nominal gap width refers to the vertical distance between the average height line of the copper layer surface and the average height line of the dielectric layer surface under no contact pressure. This value is obtained by optical profile measurement and is usually tens to hundreds of nanometers.

[0068] In this embodiment, the nominal gap width at a certain step corner is 80 nanometers. The average flattening amount of the micro-protrusion is 1.18 nanometers. Due to the contact pressure causing elastic compression of the micro-protrusion, the actual distance between the two surfaces decreases. Therefore, subtracting the average flattening amount of the micro-protrusion from the nominal gap width yields an intermediate gap width of 78.82 nanometers. The non-contact area ratio is equal to 1 minus the actual contact area ratio. The actual contact gap width is defined as the true distance between the two surfaces only in the area where no solid contact occurs. Considering that under contact pressure, the gap width in the non-contact area is approximately equal to the intermediate gap width (because only the micro-protrusion is flattened, while the height of the recessed area remains unchanged), this approach is recommended.

[0069] Therefore, the actual contact gap width is equal to the intermediate gap width multiplied by the proportion of the non-contact area. This updated actual contact gap width replaces the initial optical measurement width for thermal conduction calculations at the current time step.

[0070] The near-field radiative heat transfer thermal conductivity is recalculated based on the updated actual contact gap width. It is determined whether the gas pressure in the environmental boundary is lower than 1000 Pa or 100 Pa. If so, the gas thermal conductivity contribution is set to zero. Then, the thermal conductivity is calculated according to the near-field radiative heat transfer theory based on the gap width (78.61 nm), the surface emissivity of the copper layer (0.05), and the imaginary part of the dielectric function (2.5).

[0071] Taking the absolute temperature as 300 Kelvin, the second intermediate quantity (Boltzmann constant 1.38 x 10⁻²³ joules per Kelvin x 300 Kelvin x Planck constant 6.626 x 10⁻³⁴ joules per second) is approximately 6.25 x 10¹² per meter. The dimensionless gap parameter is equal to the gap width 78.61 nanometers x the second intermediate quantity x the speed of light 3 x 10⁸ meters per second, approximately 1.64 x 10⁻³, much less than 0.1, therefore the near-field radiation formula is used. The correction factor is calculated as the second intermediate quantity x the imaginary part of the dielectric function (2.5) x the gap width (7.861 x 10⁻⁸ meters), approximately 1.99 x 10²⁰ per meter. The near-field radiative heat transfer thermal conductivity is equal to the Stefan-Boltzmann constant 5.67 multiplied by 10 to the power of -8 (W / m² / Kelvin), multiplied by 300 Kelvin cubed (27 million Kelvin cubed), multiplied by the emissivity 0.05, and multiplied by a correction factor, yielding approximately 0.124 W / m² / Kelvin. This value is slightly higher than the 0.12 W / m² / Kelvin calculated based on a 50 nm gap, because the increased gap width slightly weakens the near-field radiation effect, but the calculated value does not change significantly.

[0072] The recalculated thermal conductivity is substituted into the thermal field equation to solve for the temperature field T and flow field V at the current time step, and then the deformation stress field M at the next time step is updated based on the temperature field T.

[0073] The new near-field radiation heat transfer thermal conductivity value of each gap element calculated in the previous step is used as the interface heat transfer coefficient and loaded into the boundary conditions of the thermal field equation, replacing the original empirical contact thermal resistance at normal pressure or the thermal conductivity value of the previous time step.

[0074] Within the current time step, based on the interface heat transfer coefficient, combined with the device power consumption value at the current time step in the power consumption time series (as a volume heat source or surface heat source), and the convective heat transfer coefficient and far-field radiation parameters in the environmental boundary, the thermal field equation is solved. The thermal field equation adopts the transient Fourier heat conduction equation, which is discretized and solved iteratively using the finite element method to obtain the temperature of each node on the circuit board at the end of the current time step, i.e., the temperature field T. Simultaneously, based on the temperature distribution in the temperature field T, the flow field V is calculated: if there is forced convection in the environmental boundary, the Navier-Stokes equations are solved according to the fan characteristic curve and the flow channel geometry to obtain the velocity field and pressure field; if it is natural convection, the buoyancy-driven flow is solved according to the density change caused by the temperature difference. In this step, the flow field V and the temperature field T are solved simultaneously through coupling, that is, they are updated alternately in each iteration step.

[0075] After obtaining the temperature field T at the current time step, it is applied as a thermal load to the mechanical analysis. The coefficients of thermal expansion of each material on the circuit board (copper, dielectric, sandwich layer, interface material, etc.) are read from the initial parameter set P. The thermal strain at each node is calculated as the coefficient of thermal expansion multiplied by the temperature change (current temperature minus reference temperature, typically 22 degrees Celsius). The thermal strain, as the initial strain, works together with the mechanical loads (preload, displacement constraints) applied by the mounting constraints to solve the linear elastic or elastoplastic constitutive equations, obtaining the deformation field and stress field at the end of the current time step, i.e., the updated deformation and stress field M.

[0076] Finally, the temperature field T sequence, flow field V sequence, and deformation stress field M sequence are output throughout the entire time history.

[0077] Please see Figure 2 As shown, based on the temperature field T, flow field V, and deformation stress field M, the local failure risk zone R at the step corner is determined, and the through-hole density, heat spreader size, interface layer thickness, and heat dissipation fin size corresponding to the risk zone R are extracted as variables to be optimized X.

[0078] The magnitude of the temperature gradient reflects the degree of drastic temperature change in space and is a key indicator for assessing local thermal stress and the risk of heat accumulation. This step uses an anisotropic three-point difference method to extract the temperature gradient magnitude at each step corner.

[0079] At each step corner, using the center node of that corner as the reference point, two adjacent nodes are selected along three specific directions. These three directions are: the copper layer extension direction (i.e., the main direction of the copper foil trace within the circuit board plane), the dielectric layer normal direction (i.e., the direction perpendicular to the circuit board plane and pointing inwards into the dielectric layer), and the direction of the angle bisector between the copper layer extension direction and the dielectric layer normal direction (i.e., the direction forming a 45-degree angle with both directions). The selected adjacent nodes should be 1 to 2 times the grid size from the reference point, typically 0.2 to 0.5 mm. In each direction, the temperature values ​​at the two nodes are recorded (read from the temperature field T, in degrees Celsius or Kelvin), and the temperature difference between the two nodes is calculated (absolute value). Simultaneously, the straight-line distance between these two nodes is measured or extracted from the grid, in millimeters. For each direction, a set of data is obtained: temperature difference and straight-line distance.

[0080] Divide the temperature difference in each direction by the corresponding straight-line distance to obtain the component gradient value in that direction, in degrees Celsius per millimeter. For the copper layer extension direction, the calculated component gradient value is denoted as component gradient one; for the dielectric layer normal direction, the calculated component gradient value is denoted as component gradient two; and for the angle bisector direction, the calculated component gradient value is denoted as component gradient three. Then, square each of these three component gradient values, and add the three squared values ​​together to obtain an intermediate sum.

[0081] The initial temperature gradient amplitude is obtained by taking the square root of the above intermediate sum, in degrees Celsius per millimeter. This initial value does not yet consider the effects of anisotropic heat transfer. Since the thermal conductivity of the copper layer along the plane is much greater than that perpendicular to the plane, the same temperature gradient produces different thermal stress effects in different directions. Therefore, the initial temperature gradient amplitude is corrected: First, the ratio of the gradient component in the interface normal direction (i.e., component gradient two) to the gradient component in the tangential direction (taking component gradient one in the copper layer extension direction) is calculated. This ratio is then taken as the logarithm to base 10 (a commonly used logarithm). Finally, the initial temperature gradient amplitude is multiplied by this logarithm to obtain the corrected temperature gradient amplitude.

[0082] For example, if component gradient two is 5 degrees Celsius per millimeter and component gradient one is 20 degrees Celsius per millimeter, the ratio is 0.25, and the commonly used logarithm is -0.602. The corrected temperature gradient amplitude is the initial value multiplied by -0.602. Because the logarithm is negative, the actual correction result will reduce the amplitude, reflecting the weakening effect of thermal stress when the normal gradient is relatively small. If the normal gradient is greater than the tangential gradient, the ratio is greater than 1, and the logarithm is positive, then the corrected amplitude will increase. Through the above calculations, the final temperature gradient amplitude at each step corner is obtained.

[0083] Equivalent stress amplitude refers to the equivalent stress value converted from multiaxial stress state to uniaxial tensile stress, and is used to assess the yield risk of materials under multiaxial stress.

[0084] This step extracts six stress components from the deformation stress field M at the center node of each step corner, including: the normal stress perpendicular to the circuit board plane (denoted as stress component one), the tangential normal stress along the copper layer extension direction (denoted as stress component two), the third normal stress along the copper layer width direction (perpendicular to the extension direction and parallel to the board surface) (denoted as stress component three), and three shear stress components (corresponding to the normal and tangential shear stresses, and the shear stress between the two tangential directions, respectively). The equivalent stress is calculated using the von Mises yield criterion.

[0085] The calculation process is as follows: First, calculate the square of any pairwise difference among the three normal stress components, i.e., the square of the difference between stress component one and stress component two, the square of the difference between stress component two and stress component three, and the square of the difference between stress component three and stress component one. Add these three squared values ​​to obtain an intermediate sum. Then, add the squares of the three shear stress components and multiply by 3 to obtain another intermediate sum. Add the above two intermediate sums, divide by 2, and finally take the square root to obtain the equivalent stress amplitude, in Pascals or Megapascals.

[0086] For each node at the staircase corner, the above calculation is performed to obtain the equivalent stress amplitude of that node. The arithmetic mean of the equivalent stress amplitudes of all nodes in this area is taken as the representative equivalent stress amplitude of that staircase corner.

[0087] The coupling risk index is obtained by multiplying the temperature gradient amplitude by the equivalent stress amplitude, and the rate of change of flow pressure drop near the step angle is extracted from the flow field V.

[0088] The coupling risk index is defined as a dimensionless product used to simultaneously reflect the combined effects of thermal and mechanical loads. It is calculated by multiplying the corrected temperature gradient amplitude (in degrees Celsius per millimeter) by the equivalent stress amplitude (in megapascals). The resulting product is the coupling risk index, expressed in degrees Celsius·MPa per millimeter. A higher index indicates that the location is simultaneously subjected to larger temperature changes and higher stress, resulting in a higher risk of failure.

[0089] Simultaneously, the rate of change of flow pressure drop near each step angle is extracted from the obtained flow field V. The rate of change of flow pressure drop is defined as the rate of change of pressure drop per unit distance along the flow direction of the cooling medium.

[0090] Near the step corner, select three consecutive pressure monitoring points along the main flow direction (determined by the inlet direction in the environmental boundary), with an interval of 1 to 2 mm between adjacent monitoring points. Record the pressure values ​​at each point (in Pascals). Calculate the pressure difference between the first and second monitoring points, and the pressure difference between the second and third monitoring points. Divide each pressure difference by the corresponding interval to obtain two local pressure drop gradient values. Then, calculate the difference between the second and first pressure drop gradients, and divide it by the average interval between the two monitoring points to obtain the flow pressure drop rate, in Pascals per square millimeter.

[0091] When the rate of change is positive, it indicates that the flow resistance increases along the flow path, which may indicate local blockage or a sudden reduction in the flow cross-sectional area; when it is negative, it indicates that the flow accelerates.

[0092] When the coupling risk index exceeds the preset threshold and the rate of change of flow pressure drop is positive, the step corner is marked as a local failure risk zone R.

[0093] The preset threshold is determined as follows: First, for all stepped corners on the target sandwich stepped circuit board, calculate their respective coupling risk indices using the method described above. Sort these indices by value from smallest to largest, and take the 75th percentile (i.e., the upper quartile) as the baseline value. Multiply this baseline value by 1.2 to obtain the preset threshold. If the circuit board has historical failure data or accelerated life test data, the threshold is set based on the minimum coupling risk index at the corresponding location when the failure occurred, and 0.8 times this minimum value is taken as the safety threshold. In this embodiment, if there is no historical data, a statistical method is used: calculate the average value of the coupling risk indices of all stepped corners, and then multiply this average value by 1.5 to obtain the preset threshold.

[0094] Each stepped corner is determined to simultaneously meet two conditions: first, the coupling risk index of the corner is greater than or equal to a preset threshold; second, the rate of change of flow pressure drop near the corner is positive (i.e., greater than 0 Pascals per square millimeter). When both conditions are met, the stepped corner is marked as a local failure risk zone R. If only one condition is met, it is not marked. The marked risk zones R are output in the form of a coordinate list or geometric region, which includes the spatial location of each risk zone, the corresponding coupling risk index value, and the flow pressure drop rate value.

[0095] The variable X to be optimized is mapped to a multiphysics coupled twin simulation model W for iterative twin simulation, and the objective function Y consisting of peak temperature, local temperature difference, warpage and flow pressure drop is obtained.

[0096] The density of through holes, the size of the heat spreader, the thickness of the interface layer, and the size of the heat dissipation fins are treated as independent design variables, and the initial variable combinations are generated by Latin hypercube sampling within their respective allowable ranges.

[0097] First, determine the physical meaning and allowable value range of the four independent design variables. Via density refers to the number of thermally conductive vias per square millimeter within the local failure risk zone R, measured in vias per square millimeter. Its allowable range is typically 2 to 8 vias per square millimeter, depending on manufacturing process limitations. Vapor pack size refers to the side length of the metal vapor chamber embedded within the circuit board on the circuit board plane, measured in millimeters. For a square vapor chamber, the allowable side length ranges from 5 to 15 millimeters. Interface layer thickness refers to the compressed thickness of the thermal interface material between the chip and the heatsink, measured in micrometers, with an allowable range of 25 to 150 micrometers. Heatsink fin size refers to the fin height of the heatsink mounted on the circuit board surface, measured in millimeters, with an allowable range of 10 to 30 millimeters.

[0098] The initial variable combinations are generated using the Latin hypercube sampling method. Latin hypercube sampling is a stratified sampling technique that can uniformly cover the entire range of values ​​in a multidimensional variable space with a smaller number of sample points. The specific implementation steps are as follows: The initial number of sample points is set to 30 (i.e., generating 30 different combinations of variables). For each design variable, its allowable range is divided into 30 non-overlapping intervals. For example, the aperture density range is 2 to 8 apertures per square millimeter, with an interval length of 0.2 apertures per square millimeter. The first interval is 2.00 to 2.20, the second interval is 2.20 to 2.40, and so on until the 30th interval is 7.80 to 8.00. A value is randomly selected from each interval as a candidate value for that variable in the corresponding sample point. The above stratified random sampling is performed independently for each of the four variables, resulting in 30 candidate value sequences of length 4. Then, each variable is randomly permuted among its 30 candidate values, using a random permutation sequence from 1 to 30. The random permutations of the four variables are combined in the same order, meaning the k-th sample point is composed of the k-th value from each of the four variable permutations.

[0099] The resulting 30 sets of variable combinations are evenly distributed in four-dimensional space, avoiding clustering issues that may occur with simple random sampling. These 30 sets of variable combinations are denoted as the initial variable combination set, with each combination represented as a vector containing the through-hole density value, the heat spreader side length value, the interface layer thickness value, and the heat sink fin height value.

[0100] Each set of variables is written into the corresponding position of the twin simulation model W by parameter substitution, and transient simulation is executed sequentially to obtain peak temperature, local temperature difference, warpage and flow pressure drop.

[0101] For each combination of variables in the initial set of variable combinations, perform the following operations: First, open the input file of the twin simulation model W (e.g., a simulation script or parametric geometry file saved in text format). Locate the parameter definitions corresponding to the four design variables in this file. Specifically: the via density parameter is located in the definition area of ​​the thermally conductive via array, typically stored as a quantity per square millimeter; the vapor chamber size parameter is located in the edge length field of the vapor chamber in the geometry modeling area; the interface layer thickness parameter is located in the thickness attribute field of the interface material; and the heat sink fin size parameter is located in the fin height field of the heat sink geometry. Replace the current values ​​of these four fields with the corresponding values ​​from the variable combinations. After the replacements are complete, save the modified input file and call the simulation solver to perform a transient simulation.

[0102] After the simulation, four target quantities are extracted from the output: peak temperature (the maximum temperature across all nodes of the circuit board, in degrees Celsius), local temperature difference (the difference between the highest and lowest temperatures within the local failure risk zone R, in degrees Celsius), warpage (the maximum displacement of the circuit board in the thickness direction, i.e., the maximum positive displacement perpendicular to the board surface in the deformation stress field M minus the minimum negative displacement, in micrometers), and flow pressure drop (the pressure loss of the cooling medium from the inlet to the outlet, i.e., the total inlet pressure minus the total outlet pressure, in Pascals). These four values ​​are used as the objective function value Y corresponding to this set of variable combinations, denoted as a four-dimensional vector. The above parameter substitution and simulation solution are performed one by one for each of the 30 initial variable combinations to obtain 30 sets of corresponding data pairs of input variables and output target values.

[0103] With the goal of minimizing peak temperature and ensuring that local temperature difference, warpage, and flow pressure drop do not exceed their respective upper limits, radial basis function interpolation is used to construct the response relationship between variables and the target. After optimization on the response surface, the variable combination is updated for the next iteration.

[0104] Radial basis function interpolation is a numerical method for fitting high-dimensional nonlinear functional relationships. This step uses radial basis functions to construct approximate functions from the four-dimensional input variables (orifice density, heat spreader size, interface layer thickness, and heat sink fin size) to each target quantity (peak temperature, local temperature difference, warpage, and flow pressure drop), called the response surface. The specific construction process is as follows: Given 30 sample points, each sample point contains an input vector (denoted as vector x, containing 4 components) and a corresponding target value (a scalar value, with response surfaces constructed independently for each of the four target values). A Gaussian function is chosen as the radial basis function, which takes the square of the Euclidean distance between any two sample point input vectors as the independent variable, and calculates the exponential function value by dividing the negative square by twice the square of the width parameter. The width parameter is 0.5 times the maximum value of the Euclidean distance between all sample points. When constructing the response surface, the target value is expressed as a weighted sum of the 30 radial basis functions, plus a linear polynomial term (including a constant term and a linear term for the four variables). The weight coefficients of each radial basis function and the coefficients of the linear polynomial are determined by solving a system of linear equations, ensuring that the predicted value of the response surface at each sample point is exactly equal to the actual target value at that sample point. This system of linear equations is solved using Gaussian elimination, and consists of 31 equations (30 radial basis functions plus 1 constant constraint) corresponding to 31 unknown coefficients.

[0105] After constructing the four response surfaces, optimization is performed in the four-dimensional variable space. The optimization objective is to minimize the peak temperature response surface value while constraining the local temperature difference response surface value to not exceed an upper limit (determined by design specifications, e.g., not exceeding 20 degrees Celsius), the warpage response surface value to not exceed an upper limit (e.g., not exceeding 50 micrometers), and the flow pressure drop response surface value to not exceed an upper limit (e.g., not exceeding 500 Pascals). A sequential quadratic programming algorithm is used to search for the combination of input variables that minimizes the peak temperature under the constraints in the four-dimensional space. The search starts with the sample point with the smallest peak temperature among the current 30 sample points, and iterative calculations stop when the variable change is less than 1% or the objective function change is less than 0.1 degrees Celsius. The new variable combination obtained from the search is used as the first set of optimization candidate combinations.

[0106] Then, the new variable combination is added to the sample point set, and the radial basis function response surface is reconstructed (the number of sample points becomes 31). The optimization is then performed again on the updated response surface to obtain the second set of optimization candidate combinations.

[0107] Repeat the above process until the difference between two consecutive optimal peak temperature predictions is less than 0.5 degrees Celsius, or the number of iterations reaches 20. The final optimal variable combination is denoted as the output variable combination. In each iteration, to verify the accuracy of the response surface prediction, after adding 5 new sample points, the best-predicted combination is selected for actual simulation verification (i.e., a complete transient simulation is performed after parameter replacement). The simulation results are compared with the response surface prediction values. If the relative error exceeds 5%, the actual simulation value is added to the sample point set to correct the response surface.

[0108] Based on the objective function Y, the optimal variable to be optimized is selected under preset constraints. The output is the same as the optimal variable to be optimized. Corresponding circuit board heat dissipation structure optimization scheme.

[0109] The preset constraints are given by the design specifications or user input, and specifically include: the upper limit of local temperature difference, the upper limit of warpage, and the upper limit of flow pressure drop.

[0110] In this embodiment, the upper limit of local temperature difference is set at 15 degrees Celsius, the upper limit of warpage is set at 40 micrometers, and the upper limit of flow pressure drop is set at 300 Pascals. For a set of candidate variable combinations, the corresponding objective function Y includes the measured local temperature difference value (denoted as measured temperature difference), the measured warpage value (denoted as measured warpage), and the measured flow pressure drop value (denoted as measured pressure drop). First, the excess ratio of local temperature difference is calculated. The measured temperature difference is subtracted from the upper limit of local temperature difference (15 degrees Celsius). If the difference is positive, it indicates that the measured temperature difference exceeds the allowable range. In this case, the difference is divided by the upper limit value (15 degrees Celsius), and the quotient is taken as the excess ratio of local temperature difference. If the difference is negative or zero, it indicates that the measured temperature difference does not exceed the upper limit. In this case, the excess ratio of local temperature difference is set to zero. Similarly, the excess ratio of warpage is calculated: the measured warpage is subtracted from the upper limit of warpage (40 micrometers). If the difference is positive, divide the difference by the upper limit (40 micrometers) to obtain the warpage excess ratio; if the difference is negative or zero, the warpage excess ratio is zero. Calculate the flow pressure drop excess ratio: Subtract the upper limit of the flow pressure drop (300 Pascals) from the measured pressure drop. If the difference is positive, divide the difference by the upper limit (300 Pascals) to obtain the flow pressure drop excess ratio; if the difference is negative or zero, the flow pressure drop excess ratio is zero.

[0111] For example, suppose a candidate variable combination corresponds to a measured temperature difference of 18 degrees Celsius, a measured warpage of 35 micrometers, and a measured pressure drop of 320 Pascals. Then, the local temperature difference excess ratio is (18 - 15) divided by 15, which equals 0.2; the warpage excess ratio is (35 - 40), which is negative, so we take 0; and the flow pressure drop excess ratio is (320 - 300) divided by 300, which is approximately 0.0667. Through the above calculations, the degree of exceeding the three constraint indicators is quantified into a dimensionless value between 0 and 1 (or greater), where zero indicates that the constraint is met, and a positive number indicates the degree of exceeding the constraint.

[0112] The weighted peak temperature is a comprehensive index derived from the measured peak temperature, after penalizing behaviors that exceed constraints based on the three excess ratios mentioned above. It is used to uniformly measure the merits of variable combinations. The specific calculation process is as follows: First, a penalty factor is constructed. The penalty factor equals 1 plus the sum of the following three terms: the first term is the excess ratio of local temperature difference multiplied by a weighting coefficient of 0.3; the second term is the excess ratio of warpage multiplied by a weighting coefficient of 0.4; and the third term is the excess ratio of flow pressure drop multiplied by a weighting coefficient of 0.3. The sum of these three weighting coefficients is 1, reflecting the degree of importance placed on local temperature difference, warpage, and flow pressure drop exceeding the constraints. In this embodiment, since excessive warpage may lead to interface delamination and stress cracking, a higher weight of 0.4 is assigned; local temperature difference and flow pressure drop are each assigned a weight of 0.3. Then, the measured peak temperature (denoted as the measured peak temperature) is multiplied by this penalty factor to obtain the weighted peak temperature. The unit of the weighted peak temperature is the same as that of the measured peak temperature, in degrees Celsius.

[0113] Using the values ​​from the previous example: assuming the measured peak temperature is 85 degrees Celsius, the penalty factor equals 1 plus (0.2 x 0.3) plus (0 x 0.4) plus (0.0667 x 0.3), which equals 1.08 (1 + 0.06 + 0 + 0.02). The weighted peak temperature equals 85 degrees Celsius multiplied by 1.08, which equals 91.8 degrees Celsius. If the measured peak temperature of another variable combination is slightly higher (e.g., 86 degrees Celsius), but all excess ratios are zero (i.e., the constraints are fully satisfied), then the penalty factor is 1, and the weighted peak temperature is 86 degrees Celsius, lower than 91.8 degrees Celsius. Therefore, this combination is better. Through this weighting method, even with a slight increase in peak temperature, as long as the constraints are avoided, it may still be selected as the optimal solution.

[0114] The combination of variables with the smallest weighted peak temperature is selected as the optimal variable to be optimized. And output an optimized solution for the heat dissipation structure of the circuit board.

[0115] For all candidate variable combinations generated iteratively (including the initial 30 groups from the Latin hypercube sampling and several new groups added in subsequent iterations, typically totaling 30 to 50 groups), calculate their respective weighted peak temperatures following the steps described above. Then, compare these weighted peak temperature values ​​and find the minimum value. The variable combination that produces this minimum value is determined as the optimal variable to be optimized, X. The optimal variable to be optimized, X, contains four specific values: optimal via density, optimal heat spreader size, optimal interface layer thickness, and optimal heat dissipation fin size.

[0116] Determine the optimal variable to be optimized Then, output the corresponding optimized heat dissipation structure scheme for the circuit board. The output format includes, but is not limited to, the following: a text-formatted list of optimized parameters, explicitly listing the optimal via density (in units per square millimeter, accurate to 0.1 units per square millimeter), the optimal heat spreader size (in millimeters, accurate to 0.5 millimeters), the optimal interface layer thickness (in micrometers, accurate to 1 micrometer), and the optimal heat sink fin size (in millimeters, accurate to 1 millimeter); and a parameter configuration file that can be directly used for production or simulation, containing adjustment values ​​for all stack-up parameters.

[0117] In addition, the predicted value of the objective function Y corresponding to the optimal solution can be output, including the values ​​of peak temperature, local temperature difference, warpage and flow pressure drop obtained from the simulation under this combination of variables, for designers to verify.

[0118] If multiple combinations of variables have the same weighted peak temperature (difference less than 0.1 degrees Celsius) during the calculation, the measured peak temperatures of each combination are further compared, and the combination with the lower measured peak temperature is selected. If the measured peak temperatures are still equal, the measured values ​​of local temperature difference, warpage, and flow pressure drop are compared sequentially, and the combination with the smaller value for each index is selected. If all comparisons fail to distinguish between the two, one combination is randomly selected as the optimal solution, and the existence of an equally optimal solution is noted in the output scheme.

[0119] Through the above steps, and taking into full account the constraints of local temperature difference, warpage, and flow pressure drop, the optimal variable to be optimized is selected from numerous candidate variable combinations based on the criterion of minimizing the weighted peak temperature. This process generates an engineering-featured optimized heat dissipation structure for circuit boards. This optimized solution can be directly used to guide the manufacturing or improvement of circuit boards.

[0120] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A circuit board heat dissipation structure optimization design method based on multi-physics coupling, characterized in that: include: Obtain the stack-up parameters, device power consumption timing, installation constraints, and environmental boundaries of the target sandwich stepped circuit board to obtain the initial parameter set P; Based on the initial parameter set P, establish an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W, which includes a copper layer, a dielectric layer, thermally conductive vias, interface materials, and heat dissipation components; Based on the power consumption timing and environmental boundary in the initial parameter set P, the twin simulation is solved using the multiphysics coupled twin simulation model W to obtain the corresponding temperature field T, flow field V and deformation stress field M. Based on the temperature field T, flow field V, and deformation stress field M, the local failure risk zone R at the step corner is determined, and the through hole density, heat spreader size, interface layer thickness, and heat dissipation fin size corresponding to the risk zone R are extracted as variables to be optimized X. The variable X to be optimized is mapped to a multiphysics coupled twin simulation model W for iterative twin simulation to obtain the objective function Y consisting of peak temperature, local temperature difference, warpage and flow pressure drop; Based on the objective function Y, the optimal variable to be optimized is selected under preset constraints. The output is the same as the optimal variable to be optimized. Corresponding circuit board heat dissipation structure optimization scheme.

2. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 1, characterized in that: Specific methods for establishing a coupled electro-thermal-fluid-mechanical twin simulation model W include: A gap grid is generated between the copper layer and the dielectric layer at the step junction, and the actual width of each gap cell is obtained by optical profile measurement. When the environmental boundary in the initial parameter set P contains a gas pressure below 1 kPa, the gas thermal conductivity contribution in the gap element is set to zero, and the near-field radiation heat transfer thermal conductivity is calculated based on the gap width and the emissivity of the copper layer surface. The near-field radiation heat transfer thermal conductivity is used as an additional interface thermal resistance boundary condition at the interface between the copper layer and the dielectric layer, and is loaded into the thermal field equation of the electro-thermal-fluid-mechanical multiphysics coupled simulation to form an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W.

3. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 2, characterized in that: Specific methods for solving twin simulations using the electro-thermal-fluid-force multiphysics coupled twin simulation model W include: The power consumption timing is divided into continuous time steps, and the normal contact pressure at the step corner in the deformation stress field M is read at the beginning of each time step. The actual contact gap width at the interface between the copper layer and the dielectric layer is updated based on the normal contact pressure, and the near-field radiative heat transfer thermal conductivity is recalculated based on the updated gap width. The recalculated thermal conductivity is substituted into the thermal field equation to solve for the temperature field T and flow field V at the current time step, and then the deformation stress field M at the next time step is updated based on the temperature field T.

4. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 3, characterized in that: The actual contact gap width at the interface between the copper layer and the dielectric layer is updated based on the normal contact pressure, including: Extract the roughness amplitude and micro-protrusion radius of curvature of the copper layer surface and the dielectric layer surface from the initial parameter set P; Substituting the normal contact pressure into the Greenwood-Williamson contact model, the actual contact area ratio and the average flattening amount of the micro-protrusions were calculated. Subtract the average flattening amount of the micro-protrusions from the nominal gap width, and then multiply by the proportion of the non-contact area to obtain the updated actual contact gap width.

5. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 1, characterized in that: Specific methods for determining the local failure risk zone R include: The temperature gradient amplitude of each step angle is extracted from the temperature field T, and the equivalent stress amplitude of each step angle is extracted from the deformation stress field M. The coupling risk index is obtained by multiplying the temperature gradient magnitude by the equivalent stress magnitude, and the rate of change of flow pressure drop near the step angle is extracted from the flow field V. When the coupling risk index exceeds the preset threshold and the rate of change of flow pressure drop is positive, the step corner is marked as a local failure risk zone R.

6. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 5, characterized in that: Extract the temperature gradient magnitudes at each step turn from the temperature field T, including: At each step corner, select two adjacent nodes along the copper layer extension direction, the dielectric layer normal direction, and the bisector of the angle between the two, and record the temperature difference and straight-line distance between the three sets of nodes. Divide each temperature difference by the corresponding straight-line distance to obtain the component gradient values ​​in three directions. Then, square the three component gradient values ​​respectively and sum them to obtain the temperature gradient amplitude.

7. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 1, characterized in that: The variable X to be optimized is mapped to an electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W for iterative twin simulation, including: The density of through holes, the size of the heat spreader, the thickness of the interface layer and the size of the heat dissipation fins are treated as independent design variables, and the initial variable combination is generated by Latin hypercube sampling within their respective allowable ranges. Each set of variables is written into the corresponding position of the electro-thermal-fluid-mechanical multiphysics coupled twin simulation model W by parameter substitution, and transient simulation is performed sequentially to obtain peak temperature, local temperature difference, warpage and flow pressure drop; With the goal of minimizing peak temperature and ensuring that local temperature difference, warpage, and flow pressure drop do not exceed their respective upper limits, radial basis function interpolation is used to construct the response relationship between variables and the target. After optimization on the response surface, the variable combination is updated for the next iteration.

8. The circuit board heat dissipation structure optimization design method based on multi-physics coupling according to claim 1, characterized in that: Based on the objective function Y, the optimal variable to be optimized is selected under preset constraints. ,include: Subtract the preset upper limit of the local temperature difference from the measured value of the local temperature difference in the objective function Y. If the difference is positive, divide the difference by the preset upper limit to obtain the local temperature difference excess ratio. Similarly, calculate the warpage excess ratio and the flow pressure drop excess ratio. If the difference is negative or zero, the excess ratio is zero. The weighted peak temperature is obtained by weighting and summing the measured peak temperature, the local temperature difference excess ratio, the warpage excess ratio, and the flow pressure drop excess ratio. The combination of variables with the smallest weighted peak temperature is selected as the optimal variable to be optimized. .

Images