Simulation method for flowing phase change heat transfer of pulse high heat flow chip

Abstract

The invention belongs to the technical field of chip heat management simulation, and particularly relates to a simulation method for flowing phase change heat transfer of a pulse high heat flow chip. Comprising the steps of physical property parameter input, unsteady state time and coulomb number management, fluid-structure interaction grid and boundary condition establishment, surface tension and temperature association, micro-scale phase change model construction, latent heat of vaporization model, wall surface roughness model and pulse heat flow time relation setting, equation set solving, speed and pressure correction, energy equation additional source item and algorithm optimization. Judging an equation residual error and a convergence condition; and carrying out time step length self-adaption and the like. According to the simulation method, a simulation algorithm is specially designed for the flowing phase change heat transfer process of a micro-channel ultrahigh-power local heat source, factors such as the non-isothermal phase change process, the micro-channel wall surface roughness and the wall surface wettability are comprehensively considered, detailed simulation of a micro-level evaporation mechanism is focused on, the flowing heat transfer characteristics of the micro-channel can be predicted more accurately, and the flow heat transfer efficiency of the micro-channel is improved. And the physical adaptability and numerical stability of the numerical method are improved.

Description

Technical Field

[0001] The present application belongs to the technical field of chip thermal management simulation, and specifically relates to a simulation method for pulse high heat flux chip flow phase change heat transfer. Background Art

[0002] When processing different types of tasks, the workload of high-performance chips will change significantly, which will cause a sharp increase in the chip's instantaneous power consumption. Since the size of the transistors that generate heat in the chip is very small, the heating surface of the microchannel presents a point-like heat source distribution. The heating surface is very small, and it is very easy to have ultra-high heat flow areas during instantaneous high-power heating. As electronic devices develop towards higher performance and smaller size, the integration of chips is getting higher and higher, and the size of components is getting smaller and smaller. For example, under 2.5D and 3D chip packaging technology, more transistors and circuits are concentrated in a limited physical space. The pulsed high heat flow during chip operation has become a key challenge for chip thermal management.

[0003] The pulsed high heat flow of the chip will suddenly increase its internal resistance, causing energy to be dissipated in the form of heat in a short period of time, making it impossible to perform effective computing and processing, and the chip operating frequency will drop rapidly in a short period of time. The pulsed high heat flow density will accelerate the aging process of chip materials, including transistor degradation, metal migration, etc. In addition, the pulsed heat flow of high-power hot spots will cause the temperature of the hot spot area to be much higher than the average temperature of the chip in a short period of time, causing greater mechanical stress due to local thermal expansion differences, further affecting chip performance and reliability.

[0004] Based on the above background, it is very important to carry out research and optimization of chip thermal management. Due to the small size of chips, the high price of high-performance chips, the high bottleneck of experimental precise temperature measurement and microfluidic control technology, and the difficulty in visualizing fluid flow and two-phase flow patterns, most chip thermal management design and development adopt numerical simulation methods. A large number of studies have shown that the application of numerical simulation technology has enabled chip thermal management technology to achieve significant and effective development. Microchannel liquid cooling and phase change heat dissipation thermal management technologies can help chips with the same thermal design power consumption (TDP) maintain higher and more stable operating frequencies when running high-complexity computing tasks for a long time, and continuously release stronger computing performance.

[0005] Traditional research considers the flow and heat transfer characteristics of the chip to achieve thermal stability during long-term operation. It is usually assumed that the wall boundary has a constant heat flux density (the first type of heat transfer boundary condition) or contains a constant internal heat source. However, the chip's pulsed high heat flux is more technically difficult, and the existing simulation technology cannot meet the needs of thermal management simulation of pulsed high heat flux chips. The main reasons are as follows:

[0006] 1) Existing common numerical simulation techniques mostly use the assumptions of isothermal phase change and constant thermal boundary conditions. However, there is a significant temperature gradient between the ultra-high power hotspot and other calculation domains. The phase change cannot be completed at a single temperature. The phase change of the fluid is a non-isothermal process. The traditional macro-based phase change equation has low prediction accuracy and cannot accurately predict the supercooled boiling phenomenon under micro-scale point-like high heat flux density.

[0007] 2) Existing common numerical simulation techniques usually take the surface tension of multiphase flow fluid as a constant, which sometimes has negligible error effects in ordinary multiphase flow heat transfer at conventional scales. However, the flow simulation of pulsed high heat flux density point heat sources needs to consider the thermal capillary pumping effect, and the constant surface tension cannot simulate the Marangoni effect caused by temperature gradient under microscale conditions.

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

[0009] 4) The fluid phase change of existing common numerical simulation technology is usually based on macroscopic physical models, and only uses the physical properties, temperature and pressure of the fluid as the criteria for judging whether a phase change occurs. However, the fluid dynamic conditions near the evaporation interface of a pulsed high heat flux density point heat source are more complex, and the evaporation rate is limited by molecular motion. It is necessary to consider the molecular dynamics process between the gas phase and the liquid phase at the microscale, and focus on the evaporation mechanism at the microscopic level.

[0010] Therefore, the existing common numerical simulation technology cannot meet the flow simulation needs of pulsed high heat flux density point heat sources at the microscale, and it is urgent to develop a simulation technology for pulsed high heat flux chip phase change thermal management to overcome the above problems. This simulation technology is expected to promote the design of advanced cooling solutions to meet the strict requirements of future high-performance electronic devices for thermal management, and has immeasurable value in improving energy efficiency and reliability, and extending the service life of electronic equipment. Summary of the invention

[0011] In order to solve at least one technical problem existing in the prior art, the present application provides a simulation method for pulse high heat flux chip flow phase change heat transfer.

[0012] The present application discloses a simulation method for pulse high heat flux chip flow phase change heat transfer, comprising the following steps:

[0013] Step S1, inputting the physical property parameters of all solids and fluids in the calculation domain for flow heat transfer simulation calculation;

[0014] Step S2, inputting the non-steady-state time, the initial time step and the maximum coulomb number;

[0015] Step S3, inputting the solid calculation domain model and the fluid calculation domain model with the boundary condition type marked;

[0016] Step S4, establishing a fluid-solid coupling calculation grid and completing the initialization of the velocity field, pressure field and temperature field;

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

[0018] Step S6, solving a group of fluid dynamics equations, wherein the surface tension, microscale phase change parameters, latent heat of vaporization, near-wall turbulent viscosity, and heat flux density involved in the group of fluid dynamics equations are obtained by traversing each grid unit in step S4 by each calculation model in step S5, and in the process of solving the group of fluid dynamics equations, the density and the inverse of the drag coefficient of each grid surface are interpolated, and then the inverse of the drag coefficient is used to achieve velocity correction;

[0019] Step S7, using the time derivative to correct the additional volume flux, the additional volume flux is the additional volume flow that appears in the numerical calculation and is introduced due to the movement of the free surface or the change of the interface shape, which is determined by the rate of phase change;

[0020] Step S8, correcting the flow rate by solving the pressure correction equation, updating the pressure field using the pressure correction amount, correcting the velocity, and calculating the turbulent kinetic energy, where the turbulent kinetic energy is the average value of the velocity fluctuation energy per unit mass of the fluid;

[0021] Step S9, introducing an additional source term into the energy equation of the fluid dynamics equations;

[0022] Step S10, performing non-orthogonal correction on the grid established in step S4;

[0023] Step S11, numerically solving the fluid dynamics equations to obtain calculation results of the heat transfer coefficient over simulation time, wherein the fluid dynamics equations at least include a mass conservation equation, a momentum conservation equation, and an energy conservation equation.

[0024] Optionally, 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.

[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 σ is the gradient of surface tension coefficient changing with temperature; T0 represents the reference temperature.

[0028] Optionally, in step S5, the microscale phase change parameter calculation model includes a mass flux correlation function and a phase change 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 represents the phase transition temperature; u v represents the gas phase velocity at the phase interface; R represents the universal gas constant; represents the gas phase volume fraction;

[0031] The correlation function of the phase change rate is as follows:

[0032]

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

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

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

[0036] Among them, L a represents 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] Among them, ν tw represents the near-wall turbulent viscosity; ν tlim =max(ν tw ,ν w ); νw represents 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 + represents the sand roughness height of the wall grid, and the expression is:

[0042]

[0043] Among them, R s Indicates the wall roughness.

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

[0045]

[0046] Among them, q max 10kW / cm 2 , t is 0~30ms;

[0047] Alternatively, the correlation function of the heat flux density calculation model is a continuous high-frequency pulse fluctuation function over 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 of the additional volume flux is as follows:

[0050]

[0051] in, represents the phase change rate; ρ v represents the gas phase density; ρ l Represents the density of the liquid phase.

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

[0053]

[0054] Among them, u', v', and w' represent the fluctuation components of the 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, and the correction source term expression is as follows:

[0056]

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

[0058] Optionally, in step S8, the pressure field is updated using the 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 represents the relaxation factor;

[0061] The updated speeds are as follows:

[0062]

[0063] in, Indicates the updated speed; represents the speed before the update; and is the relaxation factor; Δu represents the correction term of 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 the density; C p represents specific heat; represents flow rate; T represents temperature; t represents time term; represents partial derivative.

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

[0068] 1) The simulation method for pulsed high heat flux chip flow phase change heat transfer of the present application describes the phase change process by considering the molecular collision and evaporation or condensation rate between the gas phase and the liquid phase, focusing on the phase change mechanism at the microscopic level. Compared with the existing technology, it can provide a more accurate prediction of the non-isothermal phase change phenomenon under micro-scale point-like high heat flux density;

[0069] 2) The simulation method for pulsed high heat flux chip flow phase change heat transfer of the present application establishes a functional relationship between surface tension and temperature, calculates and traverses each grid unit of the fluid domain to accurately correct the surface tension of each area at the microscale. Compared with the existing technology, it can better capture the details of the thermal capillary 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 pulsed high heat flux chip flow phase change heat transfer in this application establishes the near-wall turbulent viscosity ν at the microscale tw and wall roughness R s The correlation model takes into account the complexity of turbulence, the influence of wall roughness and wall wettability on the flow characteristics in the turbulent boundary layer. Compared with traditional technologies, it can better simulate the fluid dynamic characteristics at the microscale, the growth and evolution of bubbles or droplets, and more accurately calculate the resistance pressure drop along the way;

[0071] 4) The simulation method for pulse high heat flux chip flow phase change heat transfer in this application introduces the algorithm correction term and additional source term S for the microscale flow phase change of pulse point ultra-high heat flux density heat source ad , with 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 of pulsed high heat flux chips in this application takes into account the non-isothermal phase change process of ultra-high power point pulse heat sources, the roughness of the microchannel wall, the wettability of the wall, the microscale thermal capillary pumping effect and other processes, and focuses on the calculation method or simulation technology of flow phase change heat transfer of the evaporation mechanism at the microscopic level. This technology is expected to promote the design of advanced cooling solutions, thereby meeting the strict requirements of future high-performance electronic devices for thermal management, and has immeasurable value in improving energy efficiency and reliability, and extending the service life of electronic equipment. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 It is a flowchart of the simulation technology for phase change thermal management of pulsed high heat flux chips of the present application;

[0074] Figure 2 It is a dot-line diagram of the calculation results of the correlation model of latent heat of vaporization with temperature and pressure in a preferred embodiment of the present application;

[0075] Figure 3A is a curve diagram of pulse step heat flux density changing with simulation time in a preferred embodiment of the present application;

[0076] Figure 3BIt is a curve diagram of a continuous high-frequency pulse heat flux density changing with simulation time in a preferred embodiment of the present application;

[0077] Figure 4 It is a comparison chart of the calculated results of the change of the lower microchannel wall heat transfer coefficient with the simulation time in a preferred embodiment of the present application. DETAILED DESCRIPTION

[0078] In order to make the purpose, technical solutions and advantages of the implementation of this application clearer, the technical solutions in the embodiments of this application will be described in more detail below in conjunction with the drawings in the embodiments of this application.

[0079] This application discloses a simulation method for pulse high heat flux chip flow phase change heat transfer, 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 calculation domain 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, and 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 stage, the liquid phase density of the fluid in the calculation domain is defined as 958.35 kg / m 3 , gas phase density is 0.59817kg / m 3 The gas and liquid phases of the fluid are incompatible with each other. The compressibility of the fluid is ignored under low-speed flow. 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 Pa·s, and the dynamic viscosity of the gas phase is 1.223×10 -5 Pa·s.

[0083] Step S2, input the non-steady-state time, initial time step and maximum coulomb number.

[0084] Among them, the non-steady-state time is used to determine whether the numerical calculation is terminated, the initial time step is the starting parameter of the calculation, and 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 residual and the set maximum Coulomb number during the calculation process to ensure that the global Coulomb number in the calculation domain is less than or equal to the set maximum Coulomb number, thereby avoiding cumulative errors and ensuring the stability of the calculation.

[0085] Furthermore, in this embodiment, the input unsteady-state simulation time is 40 ms and the initial time step is 1×10 -6 s, and the maximum Coulomb number is 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 topological 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 unnecessary point, line and surface features, completing the geometric topology, repairing model errors, extracting the fluid domain, and naming and marking the solid domain and the fluid domain respectively.

[0088] Step S4: Establish a fluid-solid coupling calculation grid and complete the initialization of the velocity field, pressure field and temperature field.

[0089] In this embodiment, this step specifically includes matching and coupling grids at the interface where the fluid contacts the solid, ensuring that the fluid-solid interaction can be accurately transmitted, ensuring the transmission of computational data, drawing high-quality grids and saving the msh file output in ASCII encoding, wherein all computational grids are hexahedral structured grids, the number of boundary layers is 6, and the maximum grid distortion is 0.31; finally, the velocity field, pressure field, temperature field, etc. are initialized based on the physical scene and actual boundary conditions.

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

[0091] The above calculation models are further explained below:

[0092] 1) Surface tension calculation model

[0093] This embodiment is based on experimental test data or querying the NISTREFPROP 10 standard library physical properties to obtain the corresponding data of surface tension and temperature, so as to establish the correlation function σ=f(k σ ,T,T ref ), specifically, the correlation function is as follows:

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

[0095] Wherein, σ0 represents the surface tension coefficient at the reference temperature, which is 0.0589 N / m in this embodiment; k σ is the gradient of surface tension coefficient with temperature, which is 2.014×10-4 ; T0 represents the reference temperature, which is 373.15K in this embodiment.

[0096] The above correlation function is integrated into the algorithm through code. After the calculation starts, each grid unit is traversed. Each grid node determines the size of the surface tension according to the temperature value. Microfluidics can realize the Marangoni effect caused by temperature gradient. The curvature radius r of the gas-liquid two-phase flow interface is calculated by r=f(σ,ΔP σ ) is established, where ΔP σ Representing the pressure difference due to surface tension, the algorithm can capture flow details of microscale capillary effects and thermocapillary pumping effects.

[0097] 2) Microscale phase transition parameter calculation model

[0098] This embodiment provides a microscale phase change model based on the kinetic gas theory, in which the mass flux F is expressed by considering molecular dynamics collisions, thereby capturing the fluid density ρ, diffusion coefficient D, concentration gradient ▽C, fluid molar mass M, fluid temperature T, phase change temperature T sat The gas velocity u at the interface v The combined effect on mass flux F.

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

[0100]

[0101] Where: m mo represents the mass of the molecule, k B represents the Boltzmann constant, and the molecular mass m mo Substituting for the molar mass M and combining with the ideal gas constant R, we get the average velocity u of the gas molecules mo The expression is:

[0102]

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

[0104]

[0105] Where: v represents the gas phase density; ρ l represents the liquid phase density, and R represents the universal gas constant; represents the volume fraction of the gas phase; further, the first term on the right side of the equation represents the diffusion part of the mass flux driven by the concentration gradient due to the local temperature change of the microchannel fluid under the pulsed ultra-high heat flux density, the second term on the right side of the equation takes into account the thermal motion of molecules at the microscale, the third term on the right side of the equation 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 side of the equation represent the corrections of the density term and the temperature term respectively, and the last term takes into account other influencing factors.

[0106] Furthermore, the expression of mass flux F is dimensionalized, and the dimensions on the left and right sides of the equal sign are both [kg·m -2 ·s -1 ], it is verified that the expression for the mass flux F has the correct dimension.

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

[0108]

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

[0110] 3) Calculation model of latent heat of vaporization

[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 represents latent heat of vaporization, P represents pressure, and T represents temperature;

[0114] Then, the correlation function is integrated into the algorithm through code to calculate and update the latent heat of vaporization of each grid in the calculation domain as the temperature and pressure change (that is, the latent heat of vaporization acts on each grid unit as the temperature and pressure change, and changes in real time on the spatial scale and time scale), thereby realizing accurate calculation under complex thermal boundaries and flow boundaries. After calculation, the relationship curve between the latent heat of vaporization and temperature and pressure output by the model of this embodiment is shown in the attached figure. Figure 2 As shown, the deviation from the NIST standard library data is less than 0.5%.

[0115] 4) Calculation model of near-wall turbulent viscosity

[0116] In this embodiment, the near-wall turbulent viscosity ν is established at the microscale tw and wall roughness R s The correlation model ν tw =f(R s ,k,y,ν), specifically, the correlation function of the near-wall turbulent viscosity calculation model is as follows:

[0117]

[0118] Among them, ν tw represents the near-wall turbulent viscosity; ν tlim =max(ν tw ,ν w ); ν w represents the viscosity of the fluid near the wall; dimensionless distance y represents the height of the first grid 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, represents the sand roughness height of the wall grid, and the expression is:

[0121]

[0122] Among them, R s Indicates the wall roughness.

[0123] 5) Heat flux calculation model

[0124] In this embodiment, according to the chip power and working conditions under study, a step-wise function of the pulse heat flux density q over time t is established, such as Figure 3A As shown, q max 10kW / cm 2 , t is 0~30ms, the expression is as follows:

[0125]

[0126] As an optional embodiment, a continuous high-frequency pulse fluctuation function of the heat flux q over time t can also be realized, and the expression is as follows:

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

[0128] like Figure 3BAs shown in the figure, it is the function of the heat flux density expressed by the above formula changing with time. It can be seen from the figure that the initial value of the heat flux density is 8kW / cm 2 , it reaches its maximum value during the first pulse fluctuation and pulsates continuously at a higher frequency. The overall fluctuation amplitude decays over time and gradually approaches the initial value.

[0129] The model parameters of heat flux density changing with time are adjusted according to the specific situation, acting on the specified thermal boundary condition area and integrated into the algorithm through the code.

[0130] Step S6, solving 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, momentum conservation equation and energy conservation equation. The surface tension, microscale phase change parameters, latent heat of vaporization, near-wall turbulent viscosity and heat flux density involved in the equations are obtained by traversing each grid unit in step S4 by each calculation model in step S5.

[0132] In addition, in the process of solving the fluid dynamics equations, the density and the inverse of the resistance coefficient of each grid surface are interpolated, and the inverse of the resistance coefficient is used to realize velocity correction, which 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 realizes accurate simulation of local fluid property differences (viscosity changes) caused by pulsed ultra-high heat flow.

[0133] Step S7: Correct the additional volume flux Φ using the time derivative ex , the additional volume flux is the additional volume flow introduced in the numerical calculation due to the movement of the free surface or the change of the interface shape, which is determined by the rate of phase change, and the correlation function is as follows:

[0134]

[0135] in, represents the phase change rate; ρ v represents the gas phase density; ρ l Represents the density of the liquid phase.

[0136] In this embodiment, the extra volume flux is corrected by using the time derivative based on the VOF multiphase flow model method; as time goes on, the position of the interface will change with the movement of the fluid. In order to maintain volume conservation, the change in fluid volume is corrected, so an extra volume flux is introduced in the calculation, that is, the virtual volume flux generated between units due to interface movement or deformation. This flux is not a real physical flux, but is introduced to maintain volume consistency in the calculation. The time derivative correction is to calculate this extra volume flux more accurately, reflect the movement speed of the interface through the time derivative, and ensure the correction measure of volume conservation. This process can help reduce non-physical volume changes caused by numerical errors. Ultra-high power hot spot pulses will bring huge changes to the temperature field of the calculation domain in a very short time. The correction of the time derivative improves the prediction accuracy of the rapid dynamic changes of the model.

[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, which is the average value of the velocity fluctuation energy per unit mass of the fluid, and its mathematical expression is as follows:

[0138]

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

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

[0141]

[0142] Among them, S k represents an additional source term; C μ represents a constant related to turbulent kinetic energy, which is 0.09 in this embodiment; ζ represents a constant related to the turbulent kinetic energy dissipation rate, which is 0.1 in this embodiment. k The first term takes into account the increased turbulence generation due to the presence of the phase interface when the bubble changes dynamically near the phase interface, and the second term represents the dissipation of turbulent energy caused by the gas-liquid phase change.

[0143] Further, the 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 represents the relaxation factor;

[0146] Furthermore, the velocity field is corrected again according to the above updated pressure field results, which can better meet the mass conservation, improve numerical stability and convergence speed, and the updated velocity is as follows:

[0147]

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

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

[0150] In this step, an additional source term S is introduced ad This is to prevent the numerical instability caused by convection or time derivatives that may occur during the numerical solution of pulse ultra-high heat flux point heat sources. The additional source term expression is Specifically:

[0151]

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

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

[0154] Among them, dealing with grid non-orthogonality can affect the accuracy of numerical calculations. Specifically, in actual complex geometric structures, grids are usually not completely orthogonal, especially when dealing with unstructured grids or grids near surfaces. The VOF method involves tracking and reconstruction of interfaces, and relies on solving discrete equations to capture the interface of the fluid. During the discretization process, the gradient and flux at the fluid interface need to be calculated. Non-orthogonal grids will cause errors in the calculated gradients and fluxes, thereby affecting the accuracy of interface tracking and, in turn, affecting the simulation results of the flow field. By correcting the non-orthogonality of the grid, the gradient calculation error caused by grid non-orthogonality can be corrected, the calculation accuracy of the VOF method can be improved, the numerical error caused by grid non-orthogonality can be reduced, and possible numerical instability can be avoided.

[0155] In this embodiment, the PIMPLE algorithm is used to process non-orthogonal correction, and multiple iterations are performed in each time step. A larger time step is used while ensuring the stability and accuracy of the numerical calculation, thereby improving the calculation efficiency.

[0156] In summary, this embodiment takes the non-orthogonality of the grid into consideration, and introduces non-orthogonal correction terms to correct the pressure field and the velocity field, thereby gradually reducing the error introduced by the non-orthogonality of the grid.

[0157] Step S11, numerically solving the fluid dynamics equations to obtain calculation results of the heat transfer coefficient over simulation time.

[0158] In this embodiment, the finite volume method is used to divide the computational domain into control volumes. Conservation laws are applied to the boundaries of these control volumes to transform the partial differential equations into a set of algebraic equations. The Gaussian linear interpolation method is used to calculate the flow on the control volume boundary by linear reconstruction of the surface. For a certain physical quantity φ on the control volume boundary, its value φ on the surface f f , through the value φ of the adjacent cell center i-1 and φ i+1 Linear interpolation calculation:

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

[0160] Furthermore, the upwind scheme only uses information from the upstream direction of the flow when calculating the boundary flux. It adopts the Gaussian linear interpolation algorithm to solve the smooth flow, and adopts the upwind difference format to solve the convection-dominated flow.

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

[0162] Furthermore, the simulation method for pulse high heat flux chip flow phase change heat transfer of the present application may also include the following steps:

[0163] Step S12: determine whether the equation residual and monitoring parameters meet the convergence condition. During the iterative solution process, the residual continues to decrease with the increase of the number of iterations until it reaches a preset threshold, which includes the absolute residual (less than or equal to 10 -5 ) and relative residual (the decrease ratio relative to the initial residual), and at the same time monitor the changes of key physical quantities in the flow field with the number of iterations, and execute the next step after the convergence condition is met; when the convergence condition is not met, judge whether the residual continues to decrease and monitor whether the key value tends to be stable, if the judgment result is yes, return to step S6, if the judgment result is otherwise, continue to judge whether the calculation result diverges, if the judgment result is otherwise, adjust the sub-relaxation factor and the gas-liquid interface permeability factor, return 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 judgment result is yes, return to step S4 to re-establish the calculation model;

[0164] Step S13, judging whether the numerical solution meets the output result condition, if the judgment result is yes, outputting the calculation result of the corresponding time step, if the judgment result is no, directly executing the next step operation;

[0165] Step S14, determine whether the simulation time reaches the preset time. If the judgment result is no, correct the time step according to the maximum coulomb number, return to step S6 to continue iterative calculation, and if the judgment result is yes, output the final result and end the algorithm program.

[0166] Finally, if Figure 4 As shown, a calculation result curve of the heat transfer coefficient versus simulation time in the above embodiment of the present application is obtained, and the result is compared with the true value of the heat transfer coefficient used for verification and the heat transfer coefficient value calculated by the traditional technology.

[0167] The results show that the technical solution proposed in this application is in good agreement with the verification value, indicating that the technical solution of this application can more accurately evaluate the growth movement of bubbles and the influence of ultra-high heat flux pulse changes on the microchannel wall on heat transfer than the existing technology. In addition, under the optimization of the algorithm, the simulation method of this application has strong convergence and numerical stability.

[0168] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application shall be based on the protection scope of the claims.

Claims

1. A simulation method for pulse high heat flux chip flow phase change heat transfer, characterized in that: The following steps are involved: Step S1, inputting the physical property parameters of all solids and fluids in the calculation domain for flow heat transfer simulation calculation; Step S2, inputting the non-steady-state time, the initial time step and the maximum coulomb number; Step S3, inputting the solid calculation domain model and the fluid calculation domain model with the boundary condition type marked; Step S4, establishing a fluid-solid coupling calculation grid and completing the initialization of the velocity field, pressure field and temperature field; Step S5, respectively establishing a surface tension calculation model, a microscale phase change parameter calculation model, a vaporization latent heat calculation model, a near-wall turbulent viscosity calculation model, and a heat flux density calculation model; Step S6, solving a group of fluid dynamics equations, wherein the surface tension, microscale phase change parameters, latent heat of vaporization, near-wall turbulent viscosity, and heat flux density involved in the group of fluid dynamics equations are obtained by traversing each grid unit in step S4 by each calculation model in step S5, and in the process of solving the group of fluid dynamics equations, the density and the inverse of the drag coefficient of each grid surface are interpolated, and then the inverse of the drag coefficient is used to achieve velocity correction; Step S7, using the time derivative to correct the additional volume flux, the additional volume flux is the additional volume flow introduced in the calculation of step S6 due to the movement of the free surface or the change of the interface shape, which is determined by the rate of phase change; Step S8, correcting the flow rate by solving the pressure correction equation, updating the pressure field using the pressure correction amount, correcting the velocity, and calculating the turbulent kinetic energy, where the turbulent kinetic energy is the average value of the velocity fluctuation energy per unit mass of the fluid; Step S9, introducing an additional source term into the energy equation of the fluid dynamics equations; Step S10, performing non-orthogonal correction on the grid established in step S4; Step S11, numerically solving the fluid dynamics equations to obtain calculation results of the heat transfer coefficient over simulation time, wherein the fluid dynamics equations at least include a mass conservation equation, a momentum conservation equation, and an energy conservation equation.

2. The simulation method according to claim 1, characterized in that: 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 according to claim 1, characterized in that: In step S5, the correlation function of the surface tension calculation model is as follows: σ=σ0-k σ (T-T0)(1); Where σ0 represents the surface tension coefficient at the reference temperature; k σ is the gradient of surface tension coefficient changing with temperature; T0 represents the reference temperature.

4. The simulation method according to claim 1, characterized in that: In step S5, the microscale phase change parameter calculation model includes a mass flux correlation function and a phase change rate correlation function, wherein the mass flux correlation function is as follows: Where F represents mass flux; D represents diffusion coefficient; represents the concentration gradient; M represents the molar mass of the fluid; T represents the fluid temperature; T sat represents the phase transition temperature; u v represents the gas phase velocity at the phase interface; R represents the universal gas constant; represents the gas phase volume fraction; The correlation function of the phase change rate is as follows: in, represents the phase change rate; C represents the gas-liquid interface permeability factor; A represents the area of ​​the gas-liquid interface; P v Represents the saturated vapor pressure at the phase interface; P ∞ represents the vapor pressure away from the interface; T int Represents the actual fluid temperature at the phase interface.

5. The simulation method according to claim 1, characterized in that: In step S5, the correlation function of the vaporization latent heat calculation model is as follows: L a =f(P,T)(4); Among them, L a represents latent heat of vaporization, P represents pressure, and T represents temperature.

6. The simulation method according to claim 1, characterized in that: In step S5, the correlation function of the near-wall turbulent viscosity calculation model is as follows: Among them, ν tw represents the near-wall turbulent viscosity; ν tlim =max(ν tw ,ν w ); ν w represents 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: in, represents the sand roughness height of the wall grid, and the expression is: Among them, R s Indicates the wall roughness.

7. The simulation method according to claim 1, characterized in that: In step S5, the correlation function of the heat flux calculation model is a step piecewise function with time t, and its expression is as follows: Among them, q max 10kW / cm 2 , t is 0~30ms; Alternatively, the correlation function of the heat flux density calculation model is a continuous high-frequency pulse fluctuation function over time t, and its expression is as follows: q=8e -0.5t ·sin(31.42t-4.71)+8(9)。 8. The simulation method according to claim 1, characterized in that: When calculating the turbulent kinetic energy in step S8, a correction source term is also introduced, and the correction source term expression is as follows: Among them, S k represents an additional source term; C μ represents a constant related to turbulent kinetic energy; ζ represents a constant related to the dissipation rate of turbulent kinetic energy.

9. The simulation method according to claim 1, characterized in that: In step S8, the pressure field is updated using the pressure correction value, and the updated pressure is as follows: P new =P+α P P'(11); Among them, P new Indicates the updated pressure; α P represents the relaxation factor; The updated speeds are as follows: in, Indicates the updated speed; represents the speed before the update; and is the relaxation factor; Δu represents the correction term of the velocity field.

10. The simulation method according to claim 1, characterized in that: The additional source term expression introduced in step S9 is: Among them, S ad represents the additional source term; ρ represents the density; C p represents specific heat; represents flow rate; T represents temperature; t represents time term; represents partial derivative.

Images