Lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling

Abstract

The invention relates to the technical field of multiphase flow micro-channel radiators, in particular to a lattice Boltzmann numerical simulation method for spontaneous nucleation during boiling of three-phase flow droplets. The method can be used for numerical simulation of spontaneous nucleation boiling of three-phase flow droplets in the micro-channel radiator to obtain the accurate three-phase distribution and bubble motion evolution process in the micro-channel radiator, combines an improved phase change model and a spontaneous nucleation mechanism, can effectively simulate the three-phase flow boiling process under a large density ratio, and improves the stability of the three-phase flow in the micro-channel radiator. The numerical instability of a traditional model under complex conditions is overcome; meanwhile, by introducing a new source item and a new control equation, interaction and bubble dynamic behaviors among three-phase fluids can be accurately captured.

Description

Technical Field

[0001] The invention relates to the technical field of multiphase flow microchannel radiators, in particular to a numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling based on a lattice Boltzmann flux solver. Background Art

[0002] Two-phase (steam-water) and three-component (steam-water-oil / petroleum-natural gas-water) multiphase boiling is widely used in petrochemical, energy and power, low-temperature refrigeration, electronic equipment heat dissipation and other related fields. The significant phase change ability and rapid detachment of bubbles during multiphase boiling can effectively improve the heat dissipation performance and promote heat transfer and diffusion. This is the basic principle for improving the cooling capacity of multiphase microchannel heat sinks. Dielectric fluids such as fluorocarbon fluids (FC-72, FC-87 and PF-5060) and hydrofluoroethers (HFE-7000, HFE-7100, HFE-7300) are widely used in the heat dissipation of wide-bandgap semiconductor chips due to their excellent heat dissipation performance. Such dielectric fluids usually exhibit a three-phase state during the phase change process of microchannel boiling, and their boiling nucleation mechanism and bubble evolution process are different from those of traditional multiphase flow. The study of the spontaneous nucleation mechanism during three-phase boiling and the precise directional characterization of bubble dynamics are crucial for successful applications in these fields.

[0003] A lot of results have been achieved in both experimental and numerical simulations in the research on the nucleation mechanism and behavior of three-phase flow boiling. However, the existing visualization experiments are limited to the study of macroscopic boiling phenomena and cannot characterize the processes at the microscopic scale. Traditional multiphase flow numerical simulation methods, such as the level set method and the volume of fluid (VOF) method, cannot complete the simulation of the spontaneous nucleation process of bubbles. The lattice Boltzmann method (LBM) has developed rapidly in the study of multiphase flow problems due to its characteristics of automatically capturing phase interfaces and spontaneously realizing boiling bubble nucleation. However, the existing multiphase LBM model has low numerical stability and cannot simulate the nucleation phase change problem of three-phase flow with large density ratios.

[0004] The Lattice Boltzmann Flux Solver (LBFS) has become an effective solution for numerical simulation of large density ratio three-phase flow. LBFS can effectively simulate incompressible multiphase flow with low density ratio and high density ratio. Extending LBFS to three-component fluid can effectively simulate isothermal three-component fluid flow, but the numerical stability is very poor. As a result, existing research is limited to isothermal two-phase flow or three-phase flow problems, without considering heat transfer problems and boiling phase change.

[0005] In summary, existing research mainly focuses on two-phase or three-phase flow problems under isothermal conditions. The simulation of three-phase flow boiling process with high density ratio or high viscosity ratio has the problem of poor numerical stability. Summary of the invention

[0006] The purpose of the present invention is to overcome the above-mentioned defects of the prior art, and proposes a lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling, which can be used for numerical simulation of three-phase flow droplet flow boiling in a microchannel radiator, and obtain accurate three-phase distribution and bubble motion evolution process in the microchannel radiator. The method described in the present application combines an improved phase change model and the mechanism of spontaneous nucleation, and can effectively simulate the three-phase flow boiling process under a large density ratio, overcoming the numerical instability of traditional models under complex conditions; at the same time, by introducing new source terms and control equations, it can accurately capture the interaction between the three-phase fluids and the bubble dynamics behavior.

[0007] The technical solution of the present invention is: a lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling, which includes the following steps: S1. Establish a physical model and determine the computational domain; initialize the physical parameters at the center point of each grid cell at the boundary of the computational domain, and perform initialization calculations on the physical parameters on the cell interface of the entire computational domain through linear interpolation; S2. Calculate density distribution function , temperature distribution function , phase distribution function ; S3. Calculate the chemical potential function and obtain the potential function value of the three-phase fluid , , ; S4. Calculate the interface macro flux and the newly added source term using LBFS-TPFB; S5. Use the fourth-order Runge-Kutta method to obtain the total flux , and calculate the physical parameters at the kth grid unit, including , , , , T ; Use the equation to calculate and obtain the physical parameters , , , so far the iterative calculation of the physical parameters on the interface of the computational domain unit is completed; S6. Repeat the iterative calculation process of steps S2 to S5 until a cycle of the evolution process from droplet nucleation to bubble detachment and ascent is completed.

[0008] In the present invention, in step S1, a physical model is established according to the actual physical characteristics of the chip microchannel heat sink; A cross section within the specified area of ​​the microchannel is selected, the upper boundary is an open boundary condition, the left and right boundaries are periodic boundaries, the center of the lower boundary is the heating core boundary, and the rest of the lower boundary is set as an adiabatic wall boundary.

[0009] In step S2, the density distribution function of the current node , phase distribution function and temperature distribution function is the equilibrium distribution function based on the current moment , , and the distribution function of the previous moment , , get.

[0010] Density equilibrium distribution function at the current moment for: ; in, Indicates the current position. Indicates the current moment, represents the particle direction at the current moment, represents the weight factor, Indicates pressure, represents the mass density of the fluid at a specific location and time, represents the lattice sound speed, represents the discrete velocity of the lattice, It expresses the velocity of a fluid at a specific location and time; Temperature equilibrium distribution function at the current moment for: ; Where T represents the macroscopic thermodynamic temperature at the current position; Phase balance distribution function at the current moment for: ; in, Indicates j The volume fraction of the phase fluid, Indicates j The potential function value of the phase fluid, Indicates j Thermal diffusivity of the phase fluid, Indicates j The relaxation time related to the diffusivity of the phase fluid, represents a constant; The density distribution function at the previous moment for: ; in, represents the time step, Indicates the previous moment; Temperature distribution function at the previous moment for: ; Phase distribution function at the previous moment for: ; Density distribution function The calculation formula is: ; in, Indicates j relaxation time related to the kinematic viscosity of the phase fluid; Temperature distribution function The calculation formula is: ; in, Indicates j relaxation time related to the thermal diffusivity of the phase fluid; Phase distribution function The calculation formula is: .

[0011] In step S3, the specific calculation process of the potential function value of the three-phase fluid is as follows: The simplified body free energy formula is: ; in, Indicates j Volume fraction of phase fluid; Indicates j Surface tension of the phase fluid; ; Where W represents the interface thickness; It represents the weighted average value of the surface tension of the three-phase fluid, and its calculation formula is: ; Indicates j The volume fraction of the phase fluid satisfies the functional relationship, and its calculation formula is: ; Surface tension between two phases The calculation formula is: ; Body force in the entire computational domain The calculation formula is: ; in, Indicates high density, Indicates low density, Represents the acceleration due to gravity.

[0012] In step S4, the evolution equation of LBFS-TPFB is: ; in, represents the total flux; represents the flux parallel to the x-axis; represents the flux parallel to the y-axis; Indicates the volume of a control element whose size can be adjusted; represents the enclosed area of ​​the kth grid cell; represents the source term; , represents a unit vector associated with a spatial direction; Total flux The calculation formula is: ; in, represents density, represents the vector velocity in the x direction, represents the vector velocity in the y direction, represents the volume fraction of the first phase in the mixture, represents the volume fraction of the second phase in the mixture; Flux parallel to the x-axis The calculation formula is: ; in, The discrete velocity of the lattice is x Component in the axial direction; The discrete velocity of the lattice is y Component in the axial direction; Flux parallel to the y-axis The calculation formula is: ; The calculation formula is: ; in, represents the component of surface tension in the x-axis direction, represents the component of surface tension in the y-axis direction, represents the component of the body force in the x-axis direction, represents the component of body force in the y-axis direction; Phase change source The calculation formula is:

[0013] in, represents the thermal conductivity, Indicates the ambient temperature, is the thermal diffusion coefficient, represents the phase volume fraction; The phase field equation is modified as follows: ; Where M represents thermal diffusivity; The temperature equation is modified to: ; in, and Both represent constants determined by latent heat.

[0014] In step S5, the fourth-order Runge-Kutta method is used to calculate the equation ; Get total flux ; pressure The calculation formula is: ; Fluid velocity The calculation formula is: ; Through the j Volume fraction of phase fluid The volume fraction of the first phase is obtained by the calculation formula , volume fraction of the second phase : ; temperature The calculation formula is: .

[0015] Using the calculated first phase volume fraction and the second phase volume fraction , the volume fraction of the third phase is determined by the following formula : ; density The calculation formula is: ; Viscosity The calculation formula is: ; At this point, an iterative process is completed.

[0016] A cycle of bubble evolution process includes a cycle from spontaneous nucleation boiling in droplet form, detachment and rise until leaving the calculation domain of the entire physical model; The iterative calculation process of the physical parameters on the unit interface of the physical model calculation domain is completed once within one time step. The above iterative process is repeated continuously. When a cycle of the bubble evolution process is completed, the iteration ends, and the simulation target is obtained.

[0017] The beneficial effects of the present invention are: This application is based on the lattice Boltzmann flux solver, i.e., LBFS-TPFB, to complete the spontaneous nucleation phenomenon of droplet boiling in a three-phase flow microchannel radiator. The numerical simulation method can greatly reduce the cost of manual intervention and the influence of uncontrollable factors such as experiments, thereby improving the reliability and repeatability of the results. Compared with traditional multiphase flow numerical simulation methods such as the level set method and the volume of fluid (VOF) method, the LBFS-TPFB method can effectively complete the spontaneous nucleation process of boiling bubbles, and can simulate the phase change evolution process of three-phase flow boiling in nano-micrometer cross-scale microchannel radiators, and can provide a better understanding of the mechanism of the spontaneous nucleation process of three-phase flow boiling from the mesoscopic level; The advantage of the LBFS-TPFB method over the general multiphase flow LBM model is that it can simulate the phase change boiling nucleation process of three-phase flow with a large density ratio, and has high numerical stability. A newly proposed spontaneous nucleation phase change model can easily solve the newly added source term without modifying the entire evolution equation. The present application can better simulate and present the distribution of three-phase flow in a three-phase flow microchannel radiator and reproduce the process of three-phase interface evolution.

[0018] In summary, the method of the present invention can provide a theoretical basis for the design and optimization of microchannel heat exchangers, especially in improving heat transfer efficiency and suppressing boiling instability in microchannels, and has broad application potential. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 It is a schematic diagram of the process of the method of the present invention; Figure 2 is a schematic diagram of a physical model provided by the present invention; Figure 3 Schematic diagram of the three-phase density distribution of a single three-phase droplet boiling spontaneous nucleation evolution process in a three-phase flow microchannel of the present invention; Figure 4 It is a schematic diagram of the bubble shape of LBFS-TPFB provided in this application for different Eotvos numbers; Figure 5 This is a schematic diagram of the bubble shape for different Eotvos numbers according to the LBM model proposed by Yi; Figure 6 This is a schematic diagram of the bubble shape for different Eotvos numbers in the LBM model proposed by Takada; Figure 7 This application is for The verification diagram of spontaneous nucleation of three-phase single droplet boiling under the condition of Figure 8 This application is for The verification diagram of spontaneous nucleation of three-phase single droplet boiling under the condition of Fig. 9 This application is for Verification diagram of spontaneous nucleation of three-phase single droplet boiling under the condition of . DETAILED DESCRIPTION

[0020] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0021] In the following description, specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in a variety of other ways than those described herein, and those skilled in the art can make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0022] Combine the following Figure 1 , the specific steps of a numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling based on a lattice Boltzmann flux solver described in this application are introduced in detail.

[0023] The first step is to establish a physical model based on the actual physical characteristics of the chip microchannel heat sink and determine the calculation domain.

[0024] The calculation domain of the physical model selected in this application is as follows Figure 2As shown, a cross section within a certain area of ​​the microchannel is selected, the upper boundary is an open boundary condition, and the left and right boundaries are periodic boundaries. The physical model calculation domain in this embodiment selects 150×150 grid units. The center of the lower boundary is the heating core boundary, which has a certain degree of superheat, and the other areas of the lower boundary are set as adiabatic wall boundaries. In this embodiment, there are 3 grid units of heating core boundaries at the center of the lower boundary.

[0025] The initial fluoride droplets are the second phase and are placed at the center of the lower boundary of the entire calculation domain. The initial droplet diameter is set to 60 grid units. The other areas are all liquid water phases, and the liquid water phase is the first phase. The fluoride droplet vapor bubbles generated by the initial fluoride droplets being heated and boiled are recorded as the third phase.

[0026] The second step is to set the initialization conditions at the center of the unit and the physical parameters at the boundary according to the actual situation. In other words, the physical parameters such as density, velocity, temperature, etc. are initialized directly at the center point of each grid unit at the boundary. These parameters also need to be specially processed at the boundary conditions to ensure that they conform to the actual physical situation, such as the boundary conditions of the flow inlet, outlet, solid wall, etc.

[0027] In the present application, the initial temperature of the first phase, i.e., the liquid water phase, is the saturation temperature of water, which is usually 100°C at room temperature and pressure. The temperature of the second phase, i.e., the fluoride droplets, is the saturation temperature of fluoride, which is usually 61°C at room temperature and pressure. The temperature at the heating core is between the saturation temperature of fluoride and the saturation temperature of liquid water. Therefore, the fluoride droplets boil but the liquid water does not boil. After heating, a three-phase region is generated in the entire physical model, i.e., the first phase of liquid water, the second phase of initial fluoride droplets, and the third phase of fluoride vapor generated by boiling due to heating, as shown in the attached figure. Figure 5 As shown in the figure, the green area represents the first phase, that is, the liquid water phase; the red area represents the second phase, that is, the initial fluoride droplets; and the blue area represents the third phase, that is, the fluoride vapor produced by boiling due to heat.

[0028] The third step is to initialize the calculation of the physical parameters on the unit interface of the entire physical model calculation domain through linear interpolation.

[0029] The physical parameters are dynamically calculated on the interface of the grid cells in the entire physical model calculation domain by linear interpolation, so as to determine the initial physical parameters on the grid cell interface, that is, at the junction between grid cells. The purpose is to transition the physical quantities more smoothly during the numerical calculation process, thereby improving the accuracy and stability of the simulation.

[0030] Initialization at the cell center and setting of boundary conditions are usually important steps at the beginning of simulation to ensure that the simulation starts from a reasonable initial state. The linear interpolation of physical parameters at the cell interface is an ongoing process that is involved at each time step to ensure the rationality and accuracy of the transfer of physical quantities between different cells when calculating the flow field.

[0031] Step 4: Calculate the density distribution function , temperature distribution function , phase distribution function .

[0032] Density distribution function of the current node , phase distribution function and temperature distribution function is the equilibrium distribution function based on the current moment , , and the distribution function of the previous moment , , To approximate.

[0033] and It mainly describes the density distribution and temperature distribution of the fluid, which are derived based on the fluid velocity u and corresponding physical quantities such as pressure, temperature, etc. It describes the distribution of phase field, chemical potential, etc., and is used to capture the three-phase interface and obtain the three-phase distribution.

[0034] Density equilibrium distribution function at the current moment for: ; The above function represents the value of the distribution function of particles along direction i at position x and time t when local thermodynamic equilibrium is reached. Indicates the current position. Indicates the current moment, represents the weight factor, Indicates pressure, represents the mass density of the fluid at a specific location and time, represents the lattice sound speed, represents the discrete velocity of the lattice, Represents the velocity of a fluid at a specific location and time.

[0035] Temperature equilibrium distribution function at the current moment for: ; Where T represents the macroscopic thermodynamic temperature at the current position.

[0036] Phase balance distribution function at the current moment for: ; in, Indicates j The volume fraction of the phase fluid, Indicates j The potential function value of the phase fluid, Indicates j Thermal diffusivity of the phase fluid, Indicates j The relaxation time related to the diffusivity of the phase fluid, is a constant, which is j A factor related to the phase fluid used to correct or adjust the physical effects of the distribution function.

[0037] Density distribution function at the previous moment for: ; in, represents the time step, Indicates the previous moment.

[0038] Temperature distribution function at the previous moment for: ; Phase distribution function at the previous moment for: ; Therefore, the density distribution function The calculation formula is: ; in, Indicates j The relaxation time is related to the kinematic viscosity of the phase fluid.

[0039] Temperature distribution function The calculation formula is: ; in, Indicates j The relaxation time is related to the thermal diffusivity of the phase fluid.

[0040] Phase distribution function The calculation formula is: .

[0041] Step 5: Calculate the chemical potential function and obtain the potential function value .

[0042] In multiphase flow simulation, in order to accurately describe the interface information of two-phase flow, in addition to the flow equation, it is also necessary to introduce the interface capture equation. In the diffusion interface method, the Cahn-Hilliard (CH) equation is an important and widely used tool for calculating the capture phase interface.

[0043] The body free energy formula is simplified to the following form: ; in, Indicates j Volume fraction of phase fluid; Indicates j Surface tension of the phase fluid.

[0044] ; Where W represents the interface thickness, It represents the weighted average value of the surface tension of the three-phase fluid, and its calculation formula is: ; in, Indicates the mixture j The volume fraction of the phase fluid satisfies the functional relationship, and its calculation formula is: .

[0045] Surface tension between two phases The calculation formula is: ; Body force in the entire computational domain The calculation formula is: ; in, Indicates high density, Indicates low density, Represents the acceleration due to gravity.

[0046] When one of the phases disappears and its order parameter When it is 0, = 0, the three-phase phase field equation becomes the classic two-phase CH equation, satisfying .

[0047] In the sixth step, the improved LBFS-TPFB is used to calculate the interface macro flux and the newly added source terms.

[0048] In this step, the lattice Boltzmann flux solver (abbreviated as LBFS-TPFB, the same below) is used for the LBM flux solver of three-phase boiling to solve the governing equations of three-phase fluid boiling. Compared with LBM, LBFS can easily solve the newly added source terms, such as and , without modifying the entire evolution equation. The evolution equation of LBFS-TPFB in this application can be written as: ; in, represents the total flux; represents the flux parallel to the x-axis; represents the flux parallel to the y-axis; Indicates the volume of a control element whose size can be adjusted; represents the enclosed area of ​​the kth grid cell; represents the source term; , Represents a unit vector associated with a direction in space.

[0049] Total flux The calculation formula is: ; in, represents density, represents the vector velocity in the x direction, represents the vector velocity in the y direction, represents the volume fraction of the first phase in the mixture, Represents the volume fraction of the second phase in the mixture.

[0050] Flux parallel to the x-axis The calculation formula is: ; in, The discrete velocity of the lattice is x Component in the axial direction; The discrete velocity of the lattice is y Component in the direction of the axis.

[0051] Flux parallel to the y-axis The calculation formula is: ; Source Term The calculation formula is: ; in, represents the component of surface tension in the x-axis direction, represents the component of surface tension in the y-axis direction, represents the component of the body force in the x-axis direction, represents the component of the body force in the y-axis direction.

[0052] The core idea of ​​LBFS-TPFB is to apply LBM to the magnetic flux on the interface between two control units.

[0053] In many previous studies, the phase field LBM model has been successfully applied to the study of various phase change processes. Its governing equation is obtained by adding a phase change source term to the CH equation. However, the existing original phase change model cannot realize spontaneous nucleation calculation and can only simulate the phase change in the interface region. Based on this problem, the present invention proposes a new phase change source term , as shown below: ; in, represents the thermal conductivity, Indicates the ambient temperature; is the thermal diffusion coefficient, represents the phase volume fraction.

[0054] In addition, the phase field equation can be modified as follows: ; Here, M represents the thermal diffusivity.

[0055] The temperature equation can be modified as follows: ; in, and Both represent constants determined by latent heat.

[0056] Step 7: Calculate the total flux using the fourth-order Runge-Kutta method , get the physical parameters at the kth grid unit , , , , T .

[0057] The fourth-order Runge-Kutta method is used to calculate the equation , get the total flux , and use the following equations to define the flux physical parameters, including: pressure , fluid velocity , the volume fraction of the first phase , volume fraction of the second phase , and temperature .

[0058] pressure The calculation formula is: ; Fluid velocity The calculation formula is: ; No. j Volume fraction of phase fluid The calculation formula is: ; temperature The calculation formula is: .

[0059] Step 8: Use the equation to continue calculating and obtain the physical parameters , , .

[0060] Using the physical parameters calculated in step 7 and , the volume fraction of the third phase is determined by the following formula : ; density The calculation formula is: ; Viscosity The calculation formula is: ; At this point, an iterative process is completed.

[0061] In the ninth step, repeat steps 3 to 8. Each time the processing of steps 3 to 8 is completed, the iteration of the physical parameters on the unit interface of the physical model calculation domain is completed. The above iteration process is repeated until a cycle of the evolution process from droplet nucleation to bubble detachment and rise is completed, and the iteration ends, thereby obtaining the simulation target.

[0062] In the present application, a cycle of a bubble evolution process includes a cycle from spontaneous nucleation boiling in droplet form, detachment and rise until leaving the entire physical model calculation domain.

[0063] like Figure 3As shown in the figure, the evolution process can be divided into three stages, namely, the first stage is the initial stage of droplet boiling when the liquid film has not fallen off, the second stage is the stage when the liquid film falls off and the bubble just detaches, and the third stage is the process of the bubble rising from detachment. The red area represents the initial droplet phase, the blue area represents the steam bubble phase after boiling nucleation, the green area represents the third phase, and the yellow area represents the liquid film that has not fallen off during the boiling phase change. The parameter t in the figure represents the number of iterations.

[0064] Figures 4 to 6 The figure shows the comparison of bubble shapes for different Eotvos numbers using the LBFS-TPFB of this application and the LBM model of other researchers. The top group of figures are the results of the research of this application. Figure 5 and Figure 6 The results of other researchers can be observed. As the Eotvos number increases, the bubble shape changes from round to ellipsoidal and finally to a mushroom shape. At a small Eotvos number, the degree of bubble deformation is similar to that of previous researchers, which can basically verify the accuracy of the LBFS-TPFB described in this application for the spontaneous nucleation evolution process of three-phase flow boiling.

[0065] Figures 7 to 9 The figure shows a schematic diagram of the three-phase flow droplet boiling nucleation process under different density ratios. For the simulation of the three-phase flow droplet boiling nucleation process under different density ratios, the density ratio that can be simulated by the present application reaches 1000, and the maximum density ratio can reach 2000. Figures 7 to 9 The evolution results shown can illustrate the feasibility of the method described in this application in simulating the spontaneous nucleation of droplet boiling and the evolution of the three-phase interface in a large density ratio three-phase flow.

[0066] The above is a detailed introduction to the lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling provided by the present invention. Specific examples are used in this article to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method and core ideas of the present invention. It should be pointed out that for ordinary technicians in this technical field, without departing from the principles of the present invention, the present invention can also be improved and modified in a number of ways, and these improvements and modifications also fall within the scope of protection of the claims of the present invention. The above description of the disclosed embodiments enables professional and technical personnel in this field to implement or use the present invention. Various modifications to these embodiments will be obvious to professional and technical personnel in this field, and the general principles defined in this article can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown in this article, but will conform to the widest range consistent with the principles and novel features disclosed herein.

Claims

1. A lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling, characterized in that: The following steps are involved: S1. Establish a physical model and determine the computational domain; initialize the physical parameters at the center point of each grid cell at the boundary of the computational domain, and perform initialization calculations on the physical parameters on the cell interface of the entire computational domain through linear interpolation; S2. Calculate density distribution function , temperature distribution function , phase distribution function ; S3. Calculate the chemical potential function and obtain the potential function value of the three-phase fluid , , ; S4. Calculate the interface macro flux and the newly added source term using LBFS-TPFB; S5. Use the fourth-order Runge-Kutta method to obtain the total flux , and calculate the physical parameters at the kth grid unit, including , , , , T ; Use the equation to calculate and obtain the physical parameters , , , so far the iterative calculation of the physical parameters on the interface of the computational domain unit is completed; S6. Repeat the iterative calculation process of steps S2 to S5 until a cycle of the evolution process from droplet nucleation to bubble detachment and ascent is completed.

2. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 1 is characterized in that: In step S1, a physical model is established according to the actual physical characteristics of the chip microchannel heat sink; A cross section within the specified area of ​​the microchannel is selected, the upper boundary is an open boundary condition, the left and right boundaries are periodic boundaries, the center of the lower boundary is the heating core boundary, and the rest of the lower boundary is set as an adiabatic wall boundary.

3. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 1 is characterized in that: In step S2, the density distribution function of the current node , phase distribution function and temperature distribution function is the equilibrium distribution function based on the current moment , , and the distribution function of the previous moment , , get.

4. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 3 is characterized in that: Density equilibrium distribution function at the current moment for: ； in, Indicates the current position. Indicates the current moment, represents the grid direction of the particle at the current moment, represents the weight factor, Indicates pressure, represents the mass density of the fluid at a specific location and time, represents the lattice sound speed, represents the discrete velocity of the lattice, It expresses the velocity of a fluid at a specific location and time; Temperature equilibrium distribution function at the current moment for: ； Where T represents the macroscopic thermodynamic temperature at the current position; Phase balance distribution function at the current moment for: ； in, Indicates j The volume fraction of the phase fluid, Indicates j The potential function value of the phase fluid, Indicates j Thermal diffusivity of the phase fluid, Indicates j The relaxation time related to the diffusivity of the phase fluid, represents a constant; The density distribution function at the previous moment for: ； in, represents the time step, Indicates the previous moment; Temperature distribution function at the previous moment for: ； Phase distribution function at the previous moment for: ； Density distribution function The calculation formula is: ； in, Indicates j relaxation time related to the kinematic viscosity of the phase fluid; Temperature distribution function The calculation formula is: ; in, Indicates j relaxation time related to the thermal diffusivity of the phase fluid; Phase distribution function The calculation formula is: 。 5. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 1, characterized in that: In step S3, the specific calculation process of the potential function value of the three-phase fluid is as follows: The simplified body free energy formula is: ； in, Indicates j Volume fraction of phase fluid; Indicates j Surface tension of the phase fluid; ； Where W represents the interface thickness; It represents the weighted average value of the surface tension of the three-phase fluid, and its calculation formula is: ； Indicates j The volume fraction of the phase fluid satisfies the functional relationship, and its calculation formula is: ； Surface tension between two phases The calculation formula is: ； Body force in the entire computational domain The calculation formula is: ； in, Indicates high density, Indicates low density, Represents the acceleration due to gravity.

6. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 1, characterized in that: In step S4, the evolution equation of LBFS-TPFB is: ； in, represents the total flux; represents the flux parallel to the x-axis; represents the flux parallel to the y-axis; Indicates the volume of a control element whose size can be adjusted; represents the enclosed area of ​​the kth grid cell; represents the source term; , represents a unit vector associated with a spatial direction; Total flux The calculation formula is: ； in, represents density, represents the vector velocity in the x direction, represents the vector velocity in the y direction, represents the volume fraction of the first phase in the mixture, represents the volume fraction of the second phase in the mixture; Flux parallel to the x-axis The calculation formula is: ； in, The discrete velocity of the lattice is x Component in the axial direction; The discrete velocity of the lattice is y Component in the axial direction; Flux parallel to the y-axis The calculation formula is: ； The calculation formula is: ； in, represents the component of surface tension in the x-axis direction, represents the component of surface tension in the y-axis direction, represents the component of the body force in the x-axis direction, represents the component of body force in the y-axis direction; Phase change source The calculation formula is: ； in, represents the thermal conductivity, Indicates the ambient temperature, is the thermal diffusion coefficient, represents the phase volume fraction; The phase field equation is modified as follows: ； Where M represents the thermal diffusivity; The temperature equation is modified to: ； in, and Both represent constants determined by latent heat.

7. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 1, characterized in that: In step S5, the fourth-order Runge-Kutta method is used to calculate the equation: ； Get total flux ; pressure The calculation formula is: ； Fluid velocity The calculation formula is: ； Through the j Volume fraction of phase fluid The volume fraction of the first phase is obtained by the calculation formula , volume fraction of the second phase : ； temperature The calculation formula is: 。 8. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 7, characterized in that: Using the calculated first phase volume fraction and the second phase volume fraction , the volume fraction of the third phase is determined by the following formula : ； density The calculation formula is: ； Viscosity The calculation formula is: ； At this point, an iterative process is completed.

9. The lattice Boltzmann numerical simulation method for spontaneous nucleation of three-phase flow droplet boiling according to claim 1, characterized in that: A cycle of bubble evolution process includes a cycle from spontaneous nucleation boiling in droplet form, detachment and rise until leaving the calculation domain of the entire physical model; The iterative calculation process of the physical parameters on the unit interface of the physical model calculation domain is completed once within one time step. The above iterative process is repeated continuously. When a cycle of the bubble evolution process is completed, the iteration ends, and the simulation target is obtained.

Images