A simulation method for phase change heat transfer in pulsed high heat flux chips.

Abstract

This application belongs to the field of chip thermal management simulation technology, specifically involving a simulation method for phase change heat transfer in pulsed high heat flux chips. The method includes steps such as inputting physical property parameters, managing unsteady-state time and Coulomb number, establishing fluid-structure interaction meshes and boundary conditions, correlating surface tension with temperature, constructing a microscale phase change model, establishing a latent heat of vaporization model, a wall roughness model, setting the pulsed heat flux time relationship, solving the equations, correcting velocity and pressure, adding source terms to the energy equation, optimizing the algorithm, judging equation residuals and convergence conditions, and adaptive time step. This simulation method is specifically designed for the phase change heat transfer process of ultra-high power local heat sources in microchannels. It comprehensively considers factors such as non-isothermal phase change processes, microchannel wall roughness, and wall wettability, focusing on detailed simulation of the microscopic evaporation mechanism. This allows for more accurate prediction of microchannel flow heat transfer characteristics and improves the physical adaptability and numerical stability of numerical methods.

Description

Technical Field

[0001] This application belongs to the field of chip thermal management simulation technology, specifically relating to a simulation method for phase change heat transfer in pulsed high heat flux chips. Background Technology

[0002] High-performance chips experience significant workload variations when handling different types of tasks, leading to a sharp increase in instantaneous power consumption. Due to the extremely small size of the heat-generating transistors within the chip, the microchannels exhibit a point-like heat source distribution. This extremely small surface area makes it prone to regions of ultra-high heat flux during instantaneous high-power heating. As electronic devices evolve towards higher performance and smaller sizes, chip integration is increasing, and component dimensions are shrinking. For example, with 2.5D and 3D chip packaging technologies, more transistors and circuits are concentrated within a limited physical space. The pulsed high heat flux during chip operation becomes a key challenge for chip thermal management.

[0003] High-frequency pulsed heat flux in a chip causes a sudden increase in its internal resistance, resulting in energy dissipation as heat within a short period, hindering efficient computation and causing a rapid decrease in the chip's operating frequency. High-frequency pulsed heat flux density accelerates the aging process of chip materials, including transistor degradation and metal migration. Furthermore, the pulsed heat flux from high-power hotspots causes the temperature in the hotspot area to be significantly higher than the chip's average temperature for a short period, creating substantial mechanical stress due to localized thermal expansion differences, further impacting chip performance and reliability.

[0004] Against this backdrop, research and optimization of chip thermal management are crucial. Due to the small size of chips, the high cost of high-performance chips, and the significant bottlenecks in precise experimental temperature measurement and microflow control technology, visualization of fluid flow and two-phase flow patterns is difficult. Therefore, most chip thermal management design and development employs numerical simulation methods. Extensive research has demonstrated that the application of numerical simulation technology has significantly advanced chip thermal management techniques. Microchannel liquid cooling and phase change thermal management technologies can help chips with the same thermal design power (TDP) maintain higher and more stable operating frequencies during prolonged high-complexity computational tasks, continuously delivering stronger computational performance.

[0005] Traditional research considers the flow and heat transfer characteristics of chips that have reached thermal stability during long-term operation. It usually assumes that the wall boundary has a constant heat flux density (a heat transfer boundary condition of the first kind) or contains a constant internal heat source. However, the technical difficulty of simulating pulsed high heat flux in chips is much greater. Existing simulation techniques cannot meet the needs of simulating the thermal management of chips with pulsed high heat flux. The main reasons are as follows:

[0006] 1) Existing conventional numerical simulation techniques mostly adopt the assumptions of isothermal phase transition and constant thermal boundary conditions. However, there is a significant temperature gradient between ultra-high power hotspots and other computational domains. Phase transition cannot be completed at a single temperature. The phase transition of fluids is a non-isothermal process. Traditional macroscopic phase transition equations have low prediction accuracy and cannot accurately predict the supercooled boiling phenomenon under microscale point-like high heat flux density.

[0007] 2) Existing conventional numerical simulation techniques usually take the surface tension of multiphase flow as constant. The error effect of this in conventional multiphase flow heat transfer at the conventional scale is sometimes negligible. However, the flow simulation of pulsed high heat flux density point heat sources needs to consider the thermocapillary pumping effect. Constant surface tension cannot simulate the Marangoni effect caused by temperature gradient under microscale conditions.

[0008] 3) Existing conventional numerical simulation techniques usually treat fluid channels as smooth hydraulic pipes, while the flow simulation of microchannel pulsed high heat flux density point heat sources needs to consider microscale wall effects. The wall roughness and wall wettability of microchannels are important factors affecting fluid dynamics characteristics, bubble growth and evolution mechanisms, and channel wall heat conduction.

[0009] 4) Existing conventional numerical simulation techniques for fluid phase transitions are usually based on macroscopic physical models, using only fluid properties, temperature, and pressure as criteria for judging whether a phase transition has occurred. However, the fluid dynamics near the evaporation interface of a pulsed high heat flux density point heat source are more complex. The evaporation rate is limited by molecular motion, and the molecular dynamics between the gas and liquid phases at the microscale must be considered, focusing on the evaporation mechanism at the microscopic level.

[0010] Therefore, existing conventional numerical simulation techniques cannot meet the flow simulation requirements of pulsed high heat flux density point heat sources at the microscale. There is an urgent need to develop a simulation technique for phase change thermal management of pulsed high heat flux chips to overcome these problems. This simulation technique is expected to drive the design of advanced cooling solutions, thereby meeting the stringent thermal management requirements of future high-performance electronic devices. It will have immeasurable value in improving energy efficiency and reliability, and extending the lifespan of electronic devices. Summary of the Invention

[0011] To address at least one technical problem in the prior art, this application provides a simulation method for phase change heat transfer in pulsed high heat flux chips.

[0012] This application discloses a simulation method for phase change heat transfer in pulsed high-heat-flux chips, comprising the following steps:

[0013] Step S1: Input the physical property parameters of all solids and fluids in the computational domain used for flow heat transfer simulation calculation;

[0014] Step S2: Input the unsteady-state time, initial time step, and maximum Coulomb number;

[0015] Step S3: Input the solid computational domain model and the fluid computational domain model with the boundary condition types marked;

[0016] Step S4: Establish the fluid-structure interaction computational mesh and complete the initialization of the velocity field, pressure field, and temperature field;

[0017] Step S5: Establish the surface tension calculation model, microscale phase change parameter calculation model, latent heat of vaporization calculation model, near-wall turbulent viscosity calculation model, and heat flux density calculation model respectively;

[0018] Step S6: Solve the fluid dynamics equations. The surface tension, microscale phase change parameters, latent heat of vaporization, near-wall turbulent viscosity, and heat flux density involved in the fluid dynamics equations are obtained by each calculation model in step S5 traversing each grid cell in step S4. In the process of solving the fluid dynamics equations, the density and the reciprocal of the drag coefficient of each grid surface are interpolated, and the velocity correction is achieved by using the reciprocal of the drag coefficient.

[0019] Step S7: Correct the additional volumetric flux using the time derivative. The additional volumetric flux is the additional volumetric flow that appears in the numerical calculation due to the motion of the free surface or the change of the interface shape, and it is determined by the rate of phase transition.

[0020] Step S8: Correct the flow rate by solving the pressure correction equation, update the pressure field using the pressure correction amount, correct the velocity, and calculate the turbulent kinetic energy, where the turbulent kinetic energy is the average value of the velocity fluctuation energy per unit mass of fluid.

[0021] Step S9: Introduce an additional source term into the energy equation of the fluid dynamics equation set;

[0022] Step S10: Perform non-orthogonal correction on the mesh established in step S4;

[0023] Step S11: Numerically solve the fluid dynamics equations to obtain the calculated results of the heat transfer coefficient as a function of simulation time. The fluid dynamics equations include at least the mass conservation equation, the momentum conservation equation, and the energy conservation equation.

[0024] Optionally, in step S1, the physical properties of the solid include density, specific heat capacity, and thermal conductivity, while the physical properties of the fluid include density, specific heat capacity, thermal conductivity, viscosity, and molar mass.

[0025] Optionally, in step S5, the correlation function of the surface tension calculation model is as follows:

[0026] σ=σ0-kσ (T-T0)(1);

[0027] Where σ0 represents the surface tension coefficient at the reference temperature; k σ T0 represents the gradient of surface tension coefficient with temperature; T0 represents the reference temperature.

[0028] Optionally, in step S5, the microscale phase transition parameter calculation model includes a mass flux correlation function and a phase transition rate correlation function, wherein the mass flux correlation function is as follows:

[0029]

[0030] Where F represents mass flux; D represents diffusion coefficient; ∠C represents concentration gradient; M represents fluid molar mass; T represents fluid temperature; T sat Indicates the phase transition temperature; u v R represents the gas phase velocity at the phase interface; R represents the universal gas constant. Indicates the volume fraction of the gas phase;

[0031] The correlation function for the phase transition rate is as follows:

[0032]

[0033] in, C represents the phase change rate; A represents the gas-liquid interface permeability factor; P represents the area of ​​the gas-liquid interface. v P represents the saturated vapor pressure at the phase interface. ∞ T represents the vapor pressure away from the interface. int This indicates the actual fluid temperature at the phase interface.

[0034] Optionally, in step S5, the correlation function of the latent heat of vaporization calculation model is as follows:

[0035] L a =f(P,T)(4);

[0036] Among them, L a P represents the latent heat of vaporization, P represents pressure, and T represents temperature.

[0037] Optionally, in step S5, the correlation function of the near-wall turbulent viscosity calculation model is as follows:

[0038]

[0039] Where, ν tw Indicates near-wall turbulent viscosity; ν tlim =max(ν tw ,ν w );νw Indicates the viscosity of the fluid near the wall; dimensionless distance The dimensionless velocity u* is proportional to the square root of the turbulent kinetic energy k; E is the wall roughness parameter, expressed as:

[0040]

[0041] Among them, R s + The expression representing the sand grain roughness height of the wall mesh is:

[0042]

[0043] Among them, R s This indicates the surface roughness.

[0044] Optionally, in step S5, the correlation function of the heat flux density calculation model is a step piecewise function with respect to time t, and its expression is as follows:

[0045]

[0046] Where, q max 10kW / cm 2 t is 0 to 30 ms;

[0047] Alternatively, the correlation function of the heat flux density calculation model is a continuous high-frequency pulse fluctuation function with time t, and its expression is as follows:

[0048] q = 8e -0.5t ·sin(31.42t-4.71)+8(9).

[0049] Optionally, in step S7, the correlation function for the additional volumetric flux is as follows:

[0050]

[0051] in, ρ represents the phase transition rate; v ρ represents the gas phase density. l This indicates the density of the liquid phase.

[0052] Optionally, in step S8, the mathematical expression for turbulent kinetic energy k is as follows:

[0053]

[0054] Where u', v', and w' represent the fluctuating components of fluid velocity in the x, y, and z directions, respectively.

[0055] Optionally, when calculating the turbulent kinetic energy in step S8, a correction source term is also introduced, the expression of which is as follows:

[0056]

[0057] Among them, S k Indicates additional source terms; C μ ζ represents a constant related to turbulent kinetic energy; ζ represents a constant related to the turbulent kinetic energy dissipation rate.

[0058] Optionally, in step S8, the pressure field is updated using a pressure correction amount, and the updated pressure is as follows:

[0059] P new =P+α P P'(11);

[0060] Among them, P new Indicates the updated pressure; α P Indicates the relaxation factor;

[0061] The updated speed is as follows:

[0062]

[0063] in, Indicates the updated speed; Indicates the speed before the update; and It is the relaxation factor; Δu represents the correction term for the velocity field.

[0064] Optionally, the additional source term expression introduced in step S9 is:

[0065]

[0066] Among them, S ad Represents the additional source term; ρ represents density; C p Indicates specific heat; The flow rate is represented by T; temperature is represented by t; and time is represented by t. This represents the partial derivative.

[0067] This application has at least the following beneficial technical effects:

[0068] 1) The simulation method for phase change heat transfer in pulsed high heat flux chip flow of this application describes the phase change process by considering molecular collisions and evaporation or condensation rates between the gas and liquid phases. It focuses on the phase change mechanism at the microscopic level and can provide more accurate prediction of non-isothermal phase change phenomena under microscale point high heat flux density compared with the prior art.

[0069] 2) The simulation method for phase change heat transfer in pulsed high heat flux chip flow in this application establishes the functional relationship between surface tension and temperature, calculates and accurately corrects the surface tension of each region at the microscale by traversing each grid cell of the fluid domain. Compared with the prior art, it can better capture the details of the thermocapillary pumping effect in the flow simulation of pulsed high heat flux density point heat source. Under specific flow and heat transfer boundary conditions, the Marangoni effect caused by temperature gradient under microscale conditions can be observed.

[0070] 3) The simulation method for phase change heat transfer in pulsed high heat flux chip flow proposed in this application establishes the near-wall turbulent viscosity ν at the microscale. tw With wall roughness R s The correlation model considers the influence of turbulence complexity, wall roughness and wall wettability on the flow characteristics within the turbulent boundary layer. Compared with traditional techniques, it can better simulate the fluid dynamics characteristics at the microscale, the growth and evolution of bubbles or droplets, and more accurately calculate the friction drag pressure drop.

[0071] 4) The simulation method for heat transfer of pulsed high heat flux chip flow phase change in this application introduces an algorithm correction term and an additional source term S for microscale flow phase change of pulsed point-like ultra-high heat flux density heat source. ad It exhibits smaller numerical oscillations and lower divergence probability, effectively improving the stability and accuracy of numerical calculations under complex boundary conditions;

[0072] 5) The simulation method for flow phase change heat transfer in pulsed high heat flux chips of this application considers the non-isothermal phase change process of ultra-high power point pulsed heat sources, the microchannel wall roughness, wall wettability, microscale thermocapillary pumping effect, and other processes. It focuses on the calculation method or simulation technology of flow phase change heat transfer based on the microscopic evaporation mechanism. This technology is expected to promote the design of advanced cooling solutions, thereby meeting the stringent requirements of future high-performance electronic devices for thermal management. It has immeasurable value for improving energy efficiency and reliability and extending the service life of electronic devices. Attached Figure Description

[0073] Figure 1 This is a flowchart of the analog technology for phase change thermal management of pulsed high heat flux chips according to this application;

[0074] Figure 2 This is a line graph showing the calculation results of the correlation model between latent heat of vaporization and temperature and pressure in a preferred embodiment of this application.

[0075] Figure 3A This is a graph showing the variation of pulsed step heat flux density with simulation time in a preferred embodiment of this application;

[0076] Figure 3BThis is a graph showing the change of continuous high-frequency pulsed heat flux density over simulation time in a preferred embodiment of this application.

[0077] Figure 4 This is a comparison chart of the calculated results of the change of the heat transfer coefficient of the lower microchannel wall with simulation time in a preferred embodiment of this application. Detailed Implementation

[0078] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions in the embodiments of this application will be described in more detail below with reference to the accompanying drawings.

[0079] This application discloses a simulation method for phase change heat transfer in pulsed high-heat-fluidity chips, such as... Figure 1 As shown, the simulation method includes the following steps:

[0080] Step S1: Input the physical property parameters of all solids and fluids in the computational domain used for flow heat transfer simulation calculation.

[0081] Furthermore, the physical properties of the solid in this step include density, specific heat capacity, and thermal conductivity, while the physical properties of the fluid include density, specific heat capacity, thermal conductivity, viscosity, and molar mass.

[0082] In a preferred embodiment, during this calculation preparation phase, the liquid phase density of the fluid within the computational domain is defined as 958.35 kg / m³. 3 The gas phase density is 0.59817 kg / m³. 3 The liquid and gas phases of the fluid are immiscible. Under low-speed flow, the compressibility of the fluid is neglected. The specific heat capacity of the liquid phase is 4215.7 J / (kg·K), the specific heat capacity of the gas phase is 2080.0 J / (kg·K), the thermal conductivity of the liquid phase is 0.6772 W / (m·K), the thermal conductivity of the gas phase is 0.0246 W / (m·K), the molar mass is 18.02 g / mol, and the dynamic viscosity of the liquid phase is 2.816 × 10⁻⁶. -4 The dynamic viscosity of the gas phase is 1.223 × 10 Pa·s. -5 Pa·s.

[0083] Step S2: Input the unsteady-state time, initial time step, and maximum Coulomb number.

[0084] The unsteady-state time is used to determine whether the numerical calculation has terminated. The initial time step is the starting parameter of the calculation. The Coulomb number represents the ratio of the distance traveled by the fluid (based on its characteristic velocity) to the size of the spatial grid within a time step. In this application, the time step is automatically adjusted according to the equation residuals and the set maximum Coulomb number during the calculation process to ensure that the global Coulomb number in the computational domain is less than or equal to the set maximum Coulomb number, thereby avoiding accumulated errors and ensuring the stability of the calculation.

[0085] Furthermore, in this embodiment, the input unsteady-state simulation time is 40ms, and the initial time step is 1×10⁻⁶. -6 s, with a maximum Coulomb number of 0.1.

[0086] Step S3: Input the solid computational domain model and the fluid computational domain model with the boundary condition types marked, complete the solid domain modeling and fluid domain extraction, perform topology repair on the model geometry, and mark and name the main boundaries.

[0087] In this embodiment, this step specifically includes checking the quality of the geometric model, deleting redundant points, lines, and surfaces, completing the geometric topology, fixing model errors, extracting the fluid domain, and naming and labeling the solid domain and the fluid domain respectively.

[0088] Step S4: Establish the fluid-structure interaction computational mesh and complete the initialization of the velocity field, pressure field, and temperature field.

[0089] In this embodiment, this step specifically includes matching and coupling the mesh at the interface where the fluid and solid come into contact, ensuring that the fluid-structure interaction can be accurately transmitted, guaranteeing the data transmission of the computation, drawing a high-quality mesh and saving it as an ASCII encoded msh file. The computational mesh is entirely a hexahedral structured mesh with 6 boundary layers and a maximum mesh distortion of 0.31. Finally, the initialization of the velocity field, pressure field, and temperature field is completed based on the physical scene and actual boundary conditions.

[0090] Step S5: Establish the surface tension calculation model, microscale phase change parameter calculation model, latent heat of vaporization calculation model, near-wall turbulent viscosity calculation model, and heat flux density calculation model respectively.

[0091] The above calculation models will be further explained below:

[0092] 1) Surface tension calculation model

[0093] This embodiment establishes a correlation function σ = f(k) between surface tension σ and temperature change based on experimental test data or by querying the NISTREFPROP 10 standard library. σ ,T,T ref Specifically, the association function is as follows:

[0094] σ=σ0-k σ (T-T0)(1);

[0095] Where σ0 represents the surface tension coefficient at the reference temperature, which is taken as 0.0589 N / m in this embodiment; k σ The surface tension coefficient varies with temperature; in this embodiment, it is taken as 2.014 × 10⁻⁶.-4 T0 represents the reference temperature, which is 373.15K in this embodiment.

[0096] The aforementioned correlation function is integrated into the algorithm through code. After the calculation begins, each grid cell is traversed. Each grid node determines the surface tension based on the temperature value. Microfluidics can realize the Marangoni effect caused by the temperature gradient. The radius of curvature r of the gas-liquid two-phase flow interface is obtained through r = f(σ, ΔP). σ ) is established, where ΔP σ This algorithm captures flow details caused by microscale capillary effects and thermocapillary pumping effects, representing the pressure difference due to surface tension.

[0097] 2) Calculation model for microscale phase transition parameters

[0098] This embodiment presents a microscale phase transition model based on kinetic gas theory, where the mass flux F is expressed by considering molecular dynamics collisions, thereby capturing fluid density ρ, diffusion coefficient D, concentration gradient ∇C, fluid molar mass M, fluid temperature T, and phase transition temperature T. sat Gas phase velocity u at the phase interface v The combined impact on mass flux F.

[0099] Based on the kinetic theory of gases, average kinetic energy It is directly related to temperature, as shown by the following formula:

[0100]

[0101] Where: m mo k represents the mass of a molecule. B Representing the Boltzmann constant, the molecular mass m mo Replacing it with the molar mass M and combining it with the ideal gas constant R, we obtain the average velocity u of the gas molecules. mo The expression:

[0102]

[0103] According to the Maxwell-Boltzmann distribution, the velocity distribution of molecules is a continuous probability distribution, describing the relative probability of molecules at different velocities. The expression for the mass flux F in this embodiment is as follows:

[0104]

[0105] In the formula: ρ v ρ represents the gas phase density. l R represents the liquid phase density, and R represents the universal gas constant. The gas phase volume fraction is represented by the first term on the right-hand side of the equation, which represents the diffusion portion of the mass flux driven by the concentration gradient due to the sudden change in local temperature in the microchannel fluid under pulsed ultra-high heat flux density. The second term on the right-hand side takes into account the thermal motion of molecules at the microscale. The third term on the right-hand side takes into account the influence of the gas phase velocity at the phase interface on the mass flux. The fourth and fifth terms on the right-hand side represent the corrections for the density and temperature terms, respectively. The last term takes into account other influencing factors.

[0106] Furthermore, a dimensional analysis is performed on the expression for mass flux F, and the dimensions of both the left and right sides of the equation are [kg·m]. -2 ·s -1 It has been verified that the expression for mass flux F has the correct dimensions.

[0107] At the same time, the phase transition rate at the gas-liquid interface was established. The equation (correlation function) Specifically, the association function is as follows:

[0108]

[0109] Where C represents the gas-liquid interface permeability factor; A represents the area of ​​the gas-liquid interface; P v P represents the saturated vapor pressure at the phase interface. ∞ T represents the vapor pressure away from the interface. int This indicates the actual fluid temperature at the phase interface.

[0110] 3) Latent heat of vaporization calculation model

[0111] In this embodiment, the following correlation model between latent heat of vaporization and temperature and pressure is first established:

[0112] L a =f(P,T)(4);

[0113] Among them, L a P represents the latent heat of vaporization, and T represents the pressure and temperature, respectively.

[0114] The correlation function is then integrated into the algorithm via code to realize the calculation and updating of the latent heat of vaporization of each grid in the computational domain as temperature and pressure change (i.e., the latent heat of vaporization changes with temperature and pressure in real time in both spatial and temporal scales). This enables accurate calculation under complex thermal and flow boundaries. The relationship curves between latent heat of vaporization and temperature and pressure output by the model in this embodiment are shown in the attached figure. Figure 2 As shown, the deviation from the NIST standard library data is less than 0.5%.

[0115] 4) Near-wall turbulent viscosity calculation model

[0116] In this embodiment, the near-wall turbulent viscosity ν at the microscale is established. tw With wall roughness R s Association model ν tw =f(R) s Specifically, the correlation function for the near-wall turbulent viscosity calculation model is as follows:

[0117]

[0118] Where, ν tw Indicates near-wall turbulent viscosity; ν tlim =max(ν tw ,ν w );ν w Indicates near-wall fluid viscosity; dimensionless distance y represents the height of the first mesh layer; ν represents the kinematic viscosity of the fluid near the wall, ν is a function of temperature ν = f(T); the dimensionless velocity u* is proportional to the square root of the turbulent kinetic energy k; E is the wall roughness parameter, expressed as:

[0119]

[0120] in, The expression representing the sand grain roughness height of the wall mesh is:

[0121]

[0122] Among them, R s This indicates the surface roughness.

[0123] 5) Heat flux density calculation model

[0124] In this embodiment, based on the chip power and operating conditions under study, a step piecewise function of the pulse heat flux density q versus time t is established, as follows: Figure 3A As shown, where q max 10kW / cm 2 t is 0 to 30 ms, and the expression is as follows:

[0125]

[0126] As an optional embodiment, a continuous high-frequency pulse oscillation function of heat flux density q over time t can also be realized, as expressed below:

[0127] q = 8e -0.5t ·sin(31.42t-4.71)+8(9);

[0128] like Figure 3BThe graph shows the change of heat flux density over time, as expressed by the above equation. It can be seen from the graph that the initial value of the heat flux density is 8 kW / cm². 2 It reaches its maximum value during the first pulse fluctuation and pulsates continuously at a high frequency. The overall fluctuation amplitude decays over time and gradually approaches the initial value.

[0129] Adjust the model parameters of heat flux density over time according to specific circumstances, apply them to the specified thermal boundary condition region, and integrate them into the algorithm through code.

[0130] Step S6: Solve the fluid dynamics equations.

[0131] Among them, the fluid dynamics equations are the basic equations of fluid motion, which include, but are not limited to, the mass conservation equation, the momentum conservation equation, and the energy conservation equation. The surface tension, microscale phase transition parameters, latent heat of vaporization, near-wall turbulent viscosity, and heat flux density involved in the equations are obtained by each calculation model in step S5 traversing each grid cell in step S4.

[0132] Furthermore, during the solution of the fluid dynamics equations, the density and the reciprocal of the drag coefficient of each grid surface are interpolated, and the velocity correction is achieved by using the reciprocal of the drag coefficient. This reflects the influence of microchannel flow resistance or friction loss on fluid velocity, improves the microscale physical adaptability and numerical stability of the model algorithm, and achieves accurate simulation of local fluid characteristic differences (viscosity changes) caused by pulsed ultra-high heat flux.

[0133] Step S7: Correct the additional volumetric flux Φ using the time derivative. ex The additional volumetric flux refers to the additional volumetric flow introduced in numerical calculations due to the motion of free surfaces or changes in interface shape, which is determined by the rate of phase transition, and the correlation function is as follows:

[0134]

[0135] in, ρ represents the phase transition rate; v ρ represents the gas phase density. l This indicates the density of the liquid phase.

[0136] In this embodiment, the VOF multiphase flow model method is used, and the time derivative is used to correct for additional volumetric flux. As time progresses, the position of the interface changes with fluid movement. To maintain volume conservation, the change in fluid volume is corrected. Therefore, an additional volumetric flux is introduced into the calculation, which is a virtual volumetric flux generated between cells due to interface movement or deformation. This flux is not a true physical flux, but is introduced to maintain volume consistency in the calculation. The time derivative correction is used to calculate this additional volumetric flux more accurately, reflecting the interface movement speed through the time derivative, ensuring volume conservation. This process helps reduce non-physical volume changes caused by numerical errors. Ultra-high power hotspot pulses can cause huge changes in the temperature field of the computational domain in a very short time. The time derivative correction improves the prediction accuracy of the model for rapid dynamic changes.

[0137] Step S8: Correct the flow rate by solving the pressure correction equation, update the pressure field using the pressure correction amount, correct the velocity, and calculate the turbulent kinetic energy k. The turbulent kinetic energy is the average value of the velocity fluctuation energy per unit mass of fluid, and its mathematical expression is as follows:

[0138]

[0139] Where u', v', and w' represent the fluctuating components of fluid velocity in the x, y, and z directions, respectively.

[0140] In calculating turbulent kinetic energy, this embodiment introduces a correction source term for turbulent kinetic energy to improve the model's ability to predict turbulent energy dissipation caused by dynamic changes of bubbles near the phase interface (i.e., transient flow behavior) and gas-liquid phase transitions. The expression for the correction source term is as follows:

[0141]

[0142] Among them, S k Indicates additional source terms; C μ ζ represents a constant related to turbulent kinetic energy, which is taken as 0.09 in this embodiment; ζ represents a constant related to the turbulent kinetic energy dissipation rate, which is taken as 0.1 in this embodiment. Correction source term S k The first term considers the increased turbulence generated due to the presence of the phase interface when the bubbles undergo dynamic changes near the phase interface, while the second term represents the dissipation of turbulent energy caused by the gas-liquid phase transition.

[0143] Furthermore, speed correction Related to the pressure gradient correction P', it can be expressed as The flow rate is corrected by solving the pressure correction equation, and the pressure field is updated using the pressure correction. The updated pressure is as follows:

[0144] P new =P+α PP'(11);

[0145] Among them, P new Indicates the updated pressure; α P Indicates the relaxation factor;

[0146] Furthermore, the velocity field is further corrected based on the updated pressure field results to better satisfy mass conservation, improve numerical stability and convergence speed. The updated velocity is as follows:

[0147]

[0148] in, Indicates the updated speed; Indicates the speed before the update; and It is the relaxation factor; Δu represents the correction term for the velocity field.

[0149] Step S9: Introduce an additional source term into the energy equation of the fluid dynamics equation set.

[0150] In this step, an additional source term S is introduced. ad To prevent numerical instability caused by convection or time derivative terms that may occur during the numerical solution of pulsed ultra-high heat flux density point heat sources, the additional source term expression is provided. Specifically:

[0151]

[0152] Among them, S ad The negative sign indicates an additional source term used to offset numerical instability terms in the energy equation caused by convection or time derivative terms; ρ represents density; C p Indicates specific heat; The flow rate is represented by T; temperature is represented by t; and time is represented by t. This represents the partial derivative.

[0153] Step S10: Perform non-orthogonal correction on the mesh established in step S4.

[0154] Addressing mesh non-orthogonality can significantly impact numerical computation accuracy. Specifically, in complex real-world geometries, meshes are typically not perfectly orthogonal, especially when dealing with unstructured meshes or meshes near curved surfaces. The VOF method involves interface tracking and reconstruction, relying on solving discrete equations to capture the fluid interface. During discretization, the gradient and flux at the fluid interface need to be calculated. Non-orthogonal meshes can lead to errors in the calculated gradients and fluxes, affecting the accuracy of interface tracking and consequently the simulation results. By correcting the mesh non-orthogonality, gradient calculation errors caused by mesh non-orthogonality can be corrected, improving the computational accuracy of the VOF method, reducing numerical errors caused by mesh non-orthogonality, and avoiding potential numerical instability.

[0155] In this embodiment, the PIMPLE algorithm is used to handle non-orthogonal correction. Multiple iterations are performed within each time step, which ensures the stability and accuracy of numerical calculation while using a larger time step, thereby improving computational efficiency.

[0156] In summary, this embodiment takes into account the non-orthogonality of the mesh and gradually reduces the error introduced by the non-orthogonality of the mesh by introducing a non-orthogonal correction term to correct the pressure field and velocity field.

[0157] Step S11: Solve the fluid dynamics equations numerically to obtain the calculated results of the heat transfer coefficient as a function of simulation time.

[0158] In this embodiment, the computational domain is divided into control volumes using the finite volume method. Conservation laws are applied to the boundaries of these control volumes, transforming the partial differential equations into a system of algebraic equations. The Gaussian linear interpolation method calculates the flow rate at the control volume boundaries through linear reconstruction of the surfaces. For a certain physical quantity φ at the control volume boundaries, its value φ on surface f is... f The value φ through the center of the adjacent unit i-1 and φ i+1 Linear interpolation calculation:

[0159] φ f =βφ i-1 +(1-β)φ i+1 ;

[0160] Furthermore, the upwind scheme uses only information from the upstream direction of the flow when calculating boundary fluxes, employs Gaussian linear interpolation in gentle flows, and uses upwind difference schemes when dealing with convection-dominated flows.

[0161] Furthermore, in this embodiment, the absolute residual threshold is set to 10. -5The relative residual threshold is set to 0.001. At the same time, the amount of bubble generation and the wall heat transfer rate in the flow field are monitored, and the calculation results are output every 1ms until the calculation is completed.

[0162] Furthermore, the simulation method for phase change heat transfer in pulsed high heat flux chips according to this application may also include the following steps:

[0163] Step S12: Determine whether the equation residuals and monitoring parameters meet the convergence conditions. During the iterative solution process, the residuals continuously decrease with the number of iterations until a preset threshold is reached. This threshold includes the absolute residuals (less than or equal to 10). -5 The algorithm calculates the relative residual (the rate of decrease relative to the initial residual) and monitors the changes in key physical quantities in the flow field with the number of iterations. If the convergence condition is met, the algorithm proceeds to the next step. If the convergence condition is not met, the algorithm determines whether the residual continues to decrease and whether the key values ​​tend to stabilize. If the determination result is yes, the algorithm returns to step S6. If the determination result is no, the algorithm continues to determine whether the calculation result diverges. If the determination result is no, the algorithm adjusts the sub-relaxation factor and the gas-liquid interface permeability factor and returns to step S6 to continue iterative calculation. The gas-liquid interface permeability factor is the algorithm adjustment parameter of the microscale evaporation model. If the determination result is yes, the algorithm returns to step S4 to rebuild the calculation model.

[0164] Step S13: Determine whether the numerical solution meets the output result conditions. If the determination result is yes, output the calculation result of the corresponding time step. Otherwise, proceed directly to the next step.

[0165] Step S14: Determine whether the simulation time has reached the preset time. If the result is yes, otherwise adjust the time step according to the maximum coulomb number, return to step S6 to continue iterative calculation, and if the result is yes, output the final result and end the algorithm program.

[0166] Ultimately, as Figure 4 As shown, the calculated heat transfer coefficient versus simulation time curves of the above embodiments of this application are obtained, and the results are compared with the true value of the heat transfer coefficient used for verification and the heat transfer coefficient value calculated by conventional technology.

[0167] The results show that the technical solution proposed in this application agrees well with the verification values, indicating that the technical solution of this application can more accurately evaluate the growth motion of bubbles and the impact of ultra-high heat flux pulse changes on heat transfer at the microchannel wall compared with the prior art. Furthermore, with algorithm optimization, the simulation method of this application exhibits strong convergence and numerical stability.

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

Claims

1. A simulation method for pulsed high heat flux chip flow phase change heat transfer, characterized by, The method comprises the following steps: Step S1, inputting physical property parameters of all solids and fluids in a calculation domain for flow heat transfer simulation calculation; Step S2, inputting unsteady state time, initial time step and maximum Coulomb number; Step S3, inputting solid calculation domain model and fluid calculation domain model marked with boundary condition types; Step S4, establishing fluid-structure coupling calculation grid, and initializing velocity field, pressure field and temperature field; Step S5, respectively establishing surface tension calculation model, microscale phase change parameter calculation model, vaporization latent heat calculation model, near-wall turbulent viscosity calculation model and heat flux density calculation model; Step S6, solving fluid dynamics equation set, wherein surface tension, microscale phase change parameter, vaporization latent heat, near-wall turbulent viscosity and heat flux density involved in the fluid dynamics equation set are respectively obtained by traversing each grid cell in step S4 by each calculation model in step S5, and during solving the fluid dynamics equation set, the reciprocal of the density and drag coefficient of each grid surface is interpolated, and then the reciprocal of the drag coefficient is used to correct the velocity; Step S7, correcting extra volume flux by using time derivative, wherein the extra volume flux is extra volume flow introduced due to movement of a free surface or change of an interface shape, and is determined by the rate of phase change; Step S8, correcting flow by solving a pressure correction equation, updating the pressure field by using a pressure correction amount, correcting the velocity, and calculating turbulent kinetic energy, wherein the turbulent kinetic energy is the average value of velocity fluctuation energy per unit mass of fluid; Step S9, introducing an additional source term in an energy equation of the fluid dynamics equation set; Step S10, correcting the grid established in step S4 in a non-orthogonal manner; Step S11, numerically solving the fluid dynamics equation set to obtain calculation results of heat transfer coefficients with simulation time, wherein the fluid dynamics equation set at least comprises a mass conservation equation, a momentum conservation equation and an energy conservation equation; In step S5, the correlation function of the near-wall turbulent viscosity calculation model is as follows: （5）； wherein, represents the near-wall turbulent viscosity; ; v w represents the near-wall fluid viscosity; dimensionless distance ; dimensionless velocity u is proportional to the square root of the turbulent kinetic energy; k is a wall roughness parameter, expressed as:​ （6）； wherein, represents the sand roughness height of the wall surface grid, and the expression is: （7）； wherein R s represents the wall roughness.

2. The simulation method of claim 1, wherein, In step S1, the physical property parameters of the solid include density, specific heat capacity and thermal conductivity, and the physical property parameters of the fluid include density, specific heat capacity, thermal conductivity, viscosity and molar mass.

3. The simulation method of claim 1, wherein, In step S5, the correlation function of the surface tension calculation model is as follows: （1）； wherein denotes the surface tension coefficient at the reference temperature; is the gradient of the surface tension coefficient with respect to temperature; denotes the reference temperature.

4. The simulation method of claim 1, wherein, In step S5, the microscale phase change parameter calculation model comprises a mass flux correlation function and a phase change rate correlation function, wherein the mass flux correlation function is as follows: （2）； wherein, F represents mass flux; D represents diffusion coefficient; represents concentration gradient; M represents fluid molar mass; T represents fluid temperature; T sat represents phase transition temperature; u v represents gas phase velocity at phase interface; R represents universal gas constant; The phase change rate correlation function is as follows: represents gas phase volume fraction; In step S5, the correlation function of the vaporization latent heat calculation model is as follows: （3）； wherein, represents the phase change rate; C represents the gas-liquid phase interface permeation factor; A represents the area of the gas-liquid phase interface; represents the saturated vapor pressure at the phase interface; represents the vapor pressure far from the interface; represents the actual fluid temperature at the phase interface.

5. The simulation method of claim 1, wherein, In step S8, when calculating the turbulent kinetic energy, a correction source term is also introduced, and the expression of the correction source term is as follows: （4）； wherein, represents the latent heat of vaporization, P represents the pressure, T represents the temperature.

6. The simulation method of claim 1, wherein, In said step S5, the correlation function of the heat flow density calculation model is a stepwise function in time t whose expression is as follows: （8）； wherein, is 10 kW / cm 2 , t is 0~30 ms; Alternatively, the correlation function of the heat flow density calculation model is a continuous high-frequency pulsating function over time t whose expression is as follows: （9）。 7. The simulation method of claim 1, wherein, In step S8, the pressure field is updated by using the pressure correction amount, and the updated pressure is as follows: （10）； wherein S k represents an additional source term; represents a constant related to the turbulent kinetic energy; represents a constant related to the dissipation rate of the turbulent kinetic energy.

8. The simulation method of claim 1, wherein, The updated velocity is as follows: （11）； wherein, represents the updated pressure; represents the relaxation factor; The expression of the additional source term introduced in step S9 is as follows: （12）； wherein denotes the updated velocity; denotes the velocity before the update; and is a relaxation factor; u denotes a correction term for the velocity field.

9. The simulation method of claim 1, wherein, ​ （13）； wherein, represents an additional source term; represents a density; represents a specific heat; represents a flow rate; represents a temperature; represents a time term; represents a partial derivative.

Images