Three-dimensional numerical simulation method of rayleigh-benard convection based on coupled lattice boltzmann method and finite difference method

Abstract

The application discloses a Leidenfrost effect three-dimensional numerical simulation method based on a double distribution function lattice Boltzmann method coupled with a finite difference method, and steps are as follows: 1. constructing a three-dimensional calculation domain and defining a grid point type; 2. initializing the grid point type, a fluid field density field, a velocity field, a temperature field and a distribution function; 3. calculating the fluid field density, pressure and thermal diffusion coefficient according to the distribution function; 4. calculating the resultant force of the molecular interaction force and the volume force between different phases based on a pseudo potential model; 5. calculating the fluid field collision process; 6. calculating the migration process; 7. processing the wall boundary condition; 8. updating the fluid field temperature field by using the temperature distribution function, and updating the solid field temperature field by using the finite difference method; and 9. iteratively executing steps 3 to 8 until a preset simulation time is reached. The application can stably and accurately simulate the whole process of the Leidenfrost phenomenon caused by droplet impact on a high-temperature wall surface, and provides an accurate and reliable numerical simulation tool for the design of high-efficiency heat dissipation schemes such as spray cooling.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation and calculation technology of fluid-structure conjugate heat transfer involving phase change, specifically relating to a three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method. Background Technology

[0002] The Leidenfrost effect, discovered by John Gottlob Leidenfrost in 1756, refers to the phenomenon where a vapor layer forms at the bottom of a liquid upon contact with a high-temperature surface, causing suspension and delaying evaporation. In the field of energy and power, the atomization and combustion efficiency of fuel droplets in internal combustion engines are directly affected by their interaction with the high-temperature cylinder wall. In advanced manufacturing, the spray quenching process for metallic materials relies on the rapid cooling of high-temperature metal surfaces by droplet impact to regulate their microstructure and properties. In thermal management, spray cooling for high-heat-flux-density electronic devices is a promising and highly efficient heat dissipation technology. Furthermore, in-depth research and accurate prediction of droplet impacts on high-temperature walls are crucial in scenarios such as nuclear reactor safety, inkjet printing, chip cooling, and spray fire suppression.

[0003] When the wall temperature is high, the impact of a droplet on the wall is a process accompanied by intense heat and mass transfer. A vapor layer forms between the bottom of the droplet and the high-temperature wall, resulting in the Leidenfrost phenomenon.

[0004] The impact of droplets on high-temperature walls is characterized by extremely short time scales (milliseconds or even microseconds), extremely small spatial scales (micrometers), and drastic changes in the topological structure of multiphase interfaces, posing significant challenges to studying this phenomenon using traditional experimental methods. Experimental observations are not only costly and demanding, but also struggle to capture transient details of the velocity, pressure, and temperature fields within the droplet, as well as the vaporization process near the phase interface, thus limiting a deeper understanding of the physical mechanisms underlying this phenomenon.

[0005] Therefore, numerical simulation has become an indispensable and powerful tool for studying such problems. Currently, numerical simulations for this problem mainly employ traditional numerical simulation methods based on macroscopic continuum mechanics, such as the fluid volume method or the level set method. These methods describe fluid motion by solving the Navier-Stokes equations on a fixed Eulerian grid and combine interface tracking or capture algorithms to handle the gas-liquid two-phase interface. However, existing macroscopic numerical simulation methods have the following inherent limitations when dealing with the problem of droplet impact on high-temperature walls:

[0006] Interface processing complexity: Fluid volume methods introduce numerical dissipation during interface reconstruction, leading to interface blurring; while level set methods can maintain sharp interfaces, they typically do not guarantee mass conservation. For the intense interface deformation, fragmentation, and fusion during droplet impact, the computational accuracy and robustness of these methods face severe challenges.

[0007] Phase transition model dependence: When simulating phase transition processes using macroscopic methods, it is usually necessary to introduce empirical or semi-empirical heat and mass transfer models at the interface, such as the Lee model. The phase transition coefficients in these models are often difficult to determine accurately and have poor universality, causing the simulation results to heavily depend on the selection of empirical parameters and making it difficult to accurately reflect the real physical process.

[0008] In recent years, the lattice Boltzmann method, as a mesoscale numerical method bridging macroscopic fluid dynamics and microscopic molecular dynamics, has demonstrated unique advantages in simulating complex multiphase flow problems. Based on molecular kinetic theory, the lattice Boltzmann method simulates the macroscopic behavior of fluids by solving discrete Boltzmann equations. Its unique advantages lie in its simple algorithm structure, convenient boundary condition handling, and inherently excellent parallel computing performance. For multiphase flow problems, the lattice Boltzmann method can automatically capture the generation, evolution, and even topological changes of the gas-liquid interface by introducing an interphase interaction force model, eliminating the need for complex interface tracking. The dual-distribution-function lattice Boltzmann method uses two or more independent distribution functions to solve for different macroscopic physical quantities, resulting in a more stable program.

[0009] Currently, some scholars have used the lattice Boltzmann method to study the Leidenfrost effect when droplets collide with high-temperature walls. However, two-dimensional numerical simulations cannot accurately capture the heat and mass transfer phenomena in this process. Furthermore, current analyses only cover the heating effect of the high-temperature wall on the droplets and cannot simulate the cooling effect of the droplets on the wall during this process. Summary of the Invention

[0010] To fill the research gaps in the existing technologies, the present invention aims to provide a three-dimensional numerical simulation method for the Leidenfrost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method. This method can accurately calculate the formation of thin vapor layers and droplet bouncing in the Leidenfrost phenomenon. The calculation is stable, does not easily diverge, and is highly efficient. It provides a powerful numerical simulation tool for studying the thermal physics under extreme conditions and is of great significance for guiding and optimizing the design of efficient heat dissipation schemes for spray cooling.

[0011] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0012] A three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method includes the following steps:

[0013] Step 1: Set basic parameters:

[0014] Step 1-1: Set initial physical property parameters: liquid phase density is set to... The gas phase density is set to The specific heat at constant volume is set to The liquid phase diffusion coefficient is and gas phase thermal diffusivity This is used for subsequent calculations of the fluid's thermal diffusivity; the thickness Z of the solid wall is set. solid Thermal diffusivity of solid walls Density of solids relaxation time of solids ;

[0015] Steps 1-2: Set simulation parameters: Use the D3Q19 model as the discrete model for velocity and temperature. This model contains 19 discrete velocity directions. The flow field uses a multi-relaxation time collision operator with higher numerical stability, and the temperature field uses a single-relaxation time collision operator. Construct a three-dimensional computational domain of size Nx×Ny×Nz, i.e., Nx grid points in the x-direction, Ny grid points in the y-direction, and Nz grid points in the z-direction, for a total of Nx×Ny×Nz grid points. Set the total number of simulation steps N. step ;

[0016] Step 2: Initialize the grid type, fluid domain density field, velocity field, temperature field, and distribution function:

[0017] Step 3: Calculate the fluid domain density, pressure, and thermal diffusivity:

[0018] Step 4: Calculate the resultant force of intermolecular interaction forces and volume forces between different phases:

[0019] Step 5: Computation of the fluid domain collision process:

[0020] Step 6: Calculate the migration process:

[0021] Step 7: Handle wall boundary conditions:

[0022] Step 7-1: The bottom boundary of the fluid domain, i.e., the fluid-solid conjugate heat transfer boundary, is handled using the standard half-step bounce scheme to achieve a no-slip boundary condition on the solid wall:

[0023]

[0024] In the formula:

[0025] --exist At time, located at the boundary fluid node At that point, along with discrete velocity i The input distribution function in the opposite direction, where, t For the current moment, For time step;

[0026] --exist At time, located at the boundary fluid node At that point, along the discrete velocity i The distribution function of the directional output;

[0027] Step 7-2: Top boundary of the fluid domain. Pressure boundary conditions are achieved using a non-equilibrium extrapolation scheme. This method involves setting a target density on the virtual boundary. And assuming that the non-equilibrium distribution function at the boundary is the same as that of the adjacent fluid internal grid points, the unknown distribution function flowing into the computational domain is determined. :

[0028]

[0029] In the formula:

[0030] ——Grid Input distribution function at grid point It is a boundary lattice point at k = Nz-1;

[0031] --use and The calculated equilibrium distribution function uses the pressure equilibrium boundary in this method, assuming the velocity gradient at the upper boundary is 0, and uses the velocities of adjacent internal fluid grid points. To approximate, It is the density at lattice point k=Nz-1;

[0032] —Inner fluid grid points immediately adjacent to the upper boundary Input distribution function at grid point It is a boundary lattice point at k = Nz-2;

[0033] --use and The calculated equilibrium distribution function, The lattice point represents the density at k = Nz - 2;

[0034] Step 8: Solve for the temperature field:

[0035] Calculate the temperature distribution function :

[0036] Step 8-1: Calculate the relaxation time of the temperature distribution function:

[0037] Relaxation time of the fluid domain , The fluid thermal diffusivity;

[0038] Step 8-2: Calculate the temperature distribution function using the source term-corrected lattice Boltzmann equation. The collision process:

[0039]

[0040] In the formula:

[0041] —Temperature source item;

[0042] Discrete direction i Temperature field output distribution function;

[0043] Discrete direction i Temperature field input distribution function;

[0044] —Temperature is T When, discrete direction i The temperature field equilibrium distribution function;

[0045] Step 8-3: Calculate the migration process:

[0046] The migration process of the temperature distribution function will involve the output distribution function after collision relaxation. Along discrete velocity Direction, migrate to its adjacent grid points At this point, it becomes the input distribution function for the next time step. For the x and y directions, due to the use of periodic boundary conditions, the migration is cyclical; for the z direction, when migrating to the boundary or fluid-structure interface, its direction is unknown. The value will be handled in the boundary conditions step;

[0047] The migration equation is as follows:

[0048]

[0049] In the formula:

[0050] —— time, Location, discrete direction i Temperature field input distribution function;

[0051] —— t time, Location, discrete direction i Temperature field output distribution function;

[0052] Step 8-4: Update macroscopic temperature parameters:

[0053] Update the temperature of the entire fluid domain :

[0054]

[0055] In the formula:

[0056] Discrete direction i Temperature field distribution function;

[0057] Step 8-5: Calculate the fluid-structure conjugate heat transfer boundary conditions:

[0058] On the solid side of the fluid-solid interface k s =Z solid -1 grid points fluid side of the fluid-solid interface k f =Z solid Grid points A non-equilibrium extrapolation method is used to achieve fluid-solid conjugate heat transfer, ensuring the continuity of temperature and heat flux density at the interface:

[0059] Fluid-structure interaction interface temperature ,in, It is a grid. temperature, It is a grid. Temperature;

[0060] Fluid-side temperature field input distribution function ,in It is the internal lattice point k=Z of the fluid. solid +1, For grid points Discrete direction i The corresponding temperature field input distribution function, Discrete directions at the fluid-structure interaction interface i The corresponding temperature field equilibrium distribution function, For grid points Discrete direction i The corresponding temperature field input distribution function, For grid points Discrete direction i The corresponding temperature field equilibrium distribution function;

[0061] Solid side: ,in The internal lattice point of the solid is k=Z solid -2, For grid points Discrete direction i The corresponding temperature field input distribution function, For grid points Discrete direction i The corresponding temperature field input distribution function, For grid points Discrete direction i The corresponding temperature field equilibrium distribution function;

[0062] External boundary: bottom k=0, isothermal T wb The bounce boundary scheme is adopted; the top k=Nz-1 isothermal. T s All of these are implemented using an unbalanced extrapolation scheme.

[0063] Step 8-6: Solve the heat conduction equation in the solid domain using the finite difference method:

[0064] Heat conduction equation:

[0065]

[0066] —The thermal conductivity of a solid;

[0067] T —Grid point temperature;

[0068] Discretizing the heat conduction equation using the explicit finite difference method:

[0069]

[0070] In the formula:

[0071] ——Current moment ( i , j , k Temperature of the grid points;

[0072] —The previous moment ( i , j , k The temperature of the grid points, ;

[0073] in, = =1, i , j , k These represent the grid point labels in the x, y, and z directions, respectively.

[0074] Step 9: Repeat steps 3 through 8 until the set total number of simulation steps is reached.

[0075] Preferably, step 2 is as follows:

[0076] Step 2-1: Initialize grid type:

[0077] The grid points constructed in steps 1-2 are classified using the grid identifier pointsflag, dividing the entire computational domain into internal fluid domain grid points, fluid domain boundary grid points, solid domain grid points, and virtual fluid domain grid points:

[0078] Step 2-2: Initialize the fluid domain density field:

[0079] The droplet radius is set to R0, and the droplet center position is set to (x0, y0, z0);

[0080] The gas-liquid interface width is set to ;

[0081] Calculate the distance between the grid point (x, y, z) and the droplet center position (x0, y0, z0). d ;

[0082] The initial grid density of the fluid domain is calculated using the following formula. :

[0083]

[0084] The densities of solid virtual grid points and boundary grid points are not included in the calculation, so they do not need to be set.

[0085] Step 2-3: Initialize the velocity field:

[0086] Determine the distance between the grid point (x, y, z) and the center of the droplet (x0, y0, z0). d Is the initial velocity of the lattice point (x, y, z) less than or equal to the droplet radius R0? If so, then... Set to (0, 0, -0.2), otherwise set to 0;

[0087] Step 2-4: Initialize the temperature field:

[0088] Set the temperature of the bottom solid domain grid point to T b Initialize the grid point temperature within the fluid domain to T s =0.86 T c ,in T c It is the critical temperature;

[0089] Step 2-5: Initialize the distribution function:

[0090] Initialize the temperature field equilibrium distribution function across the entire computational domain. Temperature field input distribution function and temperature field output distribution function ;

[0091] Initial velocity field equilibrium distribution function in the fluid domain Velocity field input distribution function and velocity field output distribution function ;

[0092] Initial Moment Space Equilibrium Distribution Function in Fluid Domain Moment space input distribution function Sum-moment space output distribution function .

[0093] Preferably, step 3 is as follows:

[0094] Step 3-1: Calculate the density of the fluid domain :

[0095]

[0096] Step 3-2: Calculate pressure using the equation of state The following formula gives the calculation formula for the PR state equation, or it can be replaced with other state equations for calculation:

[0097]

[0098] In the formula;

[0099] R —Gas constant;

[0100] T --temperature;

[0101] —Constant, take ;

[0102] b --constant, ;

[0103] —Critical pressure;

[0104] K — Scaling factor, used to adjust the interface thickness during simulation;

[0105] Step 3-3: Calculate the thermal diffusivity:

[0106] Within the fluid domain, the thermal diffusivity is modeled using a density-dependent linear interpolation model, and based on... and diffusion coefficient at time and Dynamic calculation:

[0107]

[0108] For fluid thermal diffusivity, in the solid domain, the solid thermal diffusivity is... Maintain its preset fixed value.

[0109] Preferably, step 4 is as follows:

[0110] Step 4-1: Calculate the potential function of the fluid domain :

[0111]

[0112] In the formula:

[0113] G —Inter-particle interaction coefficient, G =-1;

[0114] c —Speed ​​of sound;

[0115] —Grid speed of sound, ;

[0116] Step 4-2: Calculate the intermolecular interaction forces between different phases :

[0117]

[0118] In the formula:

[0119] —The i-th discrete velocity vector in the D3Q19 model;

[0120] —D3Q19 model The corresponding weighting coefficients;

[0121] —Time step;

[0122] —Current grid point position;

[0123] Step 4-3: Calculate volume forces :

[0124]

[0125] In the formula:

[0126] —The average density of the fluid domain;

[0127] g — acceleration due to gravity;

[0128] Step 4-4: Calculate the resultant force F:

[0129]

[0130] Steps 4-5: Calculation speed :

[0131]

[0132] .

[0133] Preferably, step 5 is as follows:

[0134] Step 5-1: Obtain the input distribution function of the velocity field from initialization step 2. By projecting the vector onto the moment space through an orthogonal transformation matrix M, the moment vector before the collision is obtained. :

[0135]

[0136] In the formula:

[0137] M — a 19×19 orthogonal transformation matrix;

[0138] Step 5-2:

[0139] Based on the density calculated in step 3 Calculate the corresponding equilibrium moment vector based on the velocity u calculated in step 4. These moment vectors are analytical expressions of macroscopic physical quantities in moment space, representing the local equilibrium state that the system tends to.

[0140] Step 5-3: Calculate the moment of force term:

[0141] To incorporate macroscopic forces into the evolution equations of the lattice Boltzmann method, it is necessary to transform the forces into moment space and calculate the moment of the force terms. ;

[0142] This force term moment originates from a higher-order analysis of the evolution equations of the lattice Boltzmann method, and aims to compensate for the errors introduced by discrete effects and force terms;

[0143] Step 5-4: Calculate the collision in the moment space:

[0144]

[0145] In the formula:

[0146] S – Diagonal relaxation matrix;

[0147] —A regulating term that enhances mechanical stability;

[0148] S – Diagonal relaxation matrix;

[0149] Step 5-5: Inverse transformation from moment space back to velocity space:

[0150] .

[0151] Compared with the prior art, the present invention has the following advantages:

[0152] 1. In steps 8-6 of this invention, the finite difference method is used to solve the temperature field of the solid domain. The finite difference method does not require iteration or solving large linear equations, so it is faster than using the lattice Boltzmann method.

[0153] 2. This invention uses the temperature distribution function g. i The solution avoids complex nonlinear terms in heat balance calculations, simplifies the calculation process, enhances the stability of the calculation, and can calculate fluid viscosity and thermal diffusivity separately. Furthermore, it couples the temperature fields of the solid domain and the fluid domain for solution, enabling the method to simulate conjugate heat transfer problems.

[0154] 3. This invention can clearly analyze the formation, thickness evolution, and internal flow and heat transfer details of the micron-scale vapor layer, which has scientific and engineering value for understanding and predicting the "boiling crisis" and critical heat flux density point of heat transfer efficiency. Attached Figure Description

[0155] Figure 1 The flowchart shows the three-dimensional numerical simulation method of the Leiden Frost effect based on the lattice Boltzmann method coupled with the finite difference method using the dual distribution function.

[0156] Figure 2 This is a schematic diagram of the calculation simulation;

[0157] Figure 3 Diagram showing grid type initialization;

[0158] Figure 4 This is a graph showing the change in droplet mass over time.

[0159] Figure 5 The graph shows the change of wall temperature over time at the location z=18, x=75.

[0160] Figure 6 This is a graph showing the change in droplet position over time. Detailed Implementation

[0161] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0162] like Figure 1 As shown, this invention provides a three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method, comprising the following steps:

[0163] Step 1: Set basic parameters:

[0164] Step 1-1: Set initial physical property parameters: liquid phase density is set to... The gas phase density is set to The specific heat at constant volume is set to The liquid phase diffusion coefficient is and gas phase thermal diffusivity This is used for subsequent calculations of the fluid's thermal diffusivity; the thickness Z of the solid wall is set. solid = 20, thermal diffusivity of solid wall Density of solids relaxation time of solids ;

[0165] Steps 1-2: Taking spray cooling of a high-temperature chip as an example, simulating the impact of droplets ejected from the nozzle on the high-temperature chip and the occurrence of the Leidenfrost effect: the chip and the area above the chip constitute the computational domain, such as... Figure 2 As shown,

[0166] Simulation parameters are set as follows: The D3Q19 model is used as the discrete model for velocity and temperature. This model contains 19 discrete velocity directions. The flow field uses a multi-relaxation time collision operator with higher numerical stability, while the temperature field uses a single-relaxation time collision operator. The size of the three-dimensional computational domain is Nx×Ny×Nz. The number of grid points in the x-direction is set to Nx=150, the number of grid points in the y-direction is set to Ny=150, and the number of grid points in the z-direction is set to Nz=300; that is, there are 150 grid points in the x-direction, 150 grid points in the y-direction, and 300 grid points in the z-direction, for a total of 6,750,000 grid points. The total number of simulation steps is set to Nstep=100,000.

[0167] Step 2: Initialization:

[0168] Step 2-1: Initialize the grid type, such as Figure 3 As shown:

[0169] The grid points constructed in steps 1-2 are classified using the grid identifier pointsflag, dividing the entire computational domain into internal fluid domain grid points, fluid domain boundary grid points, solid domain grid points, and virtual fluid domain grid points:

[0170] Solid domain lattice points (pointsflag=2): Solid region k ranges from 0 to Z solid The grid points are set to -2, and only the temperature field in this region is solved.

[0171] Grid points within the fluid domain (pointsflag=1): k from Z solid The flow field and temperature field are solved simultaneously at grid points from +1 to Nz-3.

[0172] Fluid domain boundary grid points (pointsflag=0): Includes fluid grid points at the top boundary k=Nz-2, used to set external fluid boundary conditions, and adjacent solid wall points k=Z. solid Fluid grid points are used to realize fluid-structure interaction boundary conditions and non-slip boundaries of the flow field;

[0173] Virtual grid points in the fluid domain (pointsflag=-2): k=Nz-1 and k=Z solid -1 grid point;

[0174] Step 2-2: Initialize the fluid domain density field:

[0175] The droplet radius is set to R0=20, and the droplet center position (x0, y0, z0) is set to x0=Nx / 2.0, y0=Ny / 2.0, z0= Nz / 2.0;

[0176] The gas-liquid interface width is set to w d =3;

[0177] Calculate the distance between the grid point (x, y, z) and the droplet center position (x0, y0, z0). d :

[0178]

[0179] The initial grid density of the fluid domain in step 2-2 is calculated using the following formula. :

[0180]

[0181] The lattice density of the solid domain is initialized to ;

[0182] Step 2-3: Initialize the velocity field:

[0183] Determine the distance between the grid point (x, y, z) and the center of the droplet (x0, y0, z0). d Is it less than or equal to the droplet radius? R 0, if so, then the initial velocity of the lattice point (x, y, z) is... Set to (0, 0, -0.2), otherwise set to 0;

[0184] Step 2-4: Initialize the temperature field:

[0185] Set the temperature of the bottom solid domain grid point to T b =1.16 T c Initialize the grid point temperature within the fluid domain to T s =0.86 T c ,in T c It is the critical temperature;

[0186] Steps 2-5: Initialize the distribution functions in the lattice Boltzmann method, including the temperature distribution function and the velocity distribution function:

[0187] Initialize the temperature field equilibrium distribution function Input distribution function and output distribution function :

[0188]

[0189]

[0190]

[0191] In the formula:

[0192] —D3Q19 model and discrete velocity vector The corresponding weighting coefficients;

[0193] —The i-th discrete velocity vector in the D3Q19 model;

[0194] T —Grid temperature, fluid domain grid T = Ts solid domain lattice T = Tb ;

[0195] Initial velocity field equilibrium distribution function in the fluid domain Velocity field input distribution function and velocity field output distribution function ;

[0196]

[0197] In the formula:

[0198] c s — Grid speed of sound;

[0199] —Vector velocity, ;

[0200]

[0201]

[0202] Initial Moment Space Equilibrium Distribution Function in Fluid Domain Moment space input distribution function Sum-moment space output distribution function ;

[0203]

[0204]

[0205]

[0206] Step 3: Calculate the fluid domain density, pressure, and thermal diffusivity:

[0207] Step 3-1: Calculate density :

[0208]

[0209] Step 3-2: Calculate pressure using the equation of state The following formula gives the calculation formula for the PR state equation, or it can be replaced with other state equations for calculation:

[0210]

[0211] In the formula;

[0212] R —Gas constant;

[0213] —Constant, take ;

[0214] b --constant, ;

[0215] —Critical pressure;

[0216] K — Scaling factor, used to adjust the interface thickness during simulation.

[0217] Step 3-3: Calculate the thermal diffusivity:

[0218] Within the fluid domain, the thermal diffusivity is modeled using a density-dependent linear interpolation model, and based on... and diffusion coefficient at time and Dynamic calculation:

[0219]

[0220] For fluid thermal diffusivity, in the solid domain, the solid thermal diffusivity is... Maintain its preset fixed value;

[0221] Step 4: Calculate the forces:

[0222] Step 4-1: Calculate the potential function of the fluid domain :

[0223]

[0224] In the formula:

[0225] G—Inter-particle interaction coefficient, G=-1;

[0226] c —Speed ​​of sound;

[0227] —Grid speed of sound, ;

[0228] Step 4-2: Calculate the intermolecular interaction forces between different phases :

[0229]

[0230] In the formula:

[0231] —Time step;

[0232] —Current grid point position;

[0233] Step 4-3: Calculate volume forces :

[0234]

[0235] In the formula:

[0236] —The average density of the fluid domain;

[0237] g — acceleration due to gravity;

[0238] Step 4-4: Calculate the resultant force F:

[0239]

[0240] Steps 4-5: Calculation speed :

[0241]

[0242]

[0243] Step 5: Computation of the fluid domain collision process:

[0244] Step 5-1: Obtain the input distribution function of the velocity field from initialization step 2. By projecting the vector onto the moment space through an orthogonal transformation matrix M, the moment vector before the collision is obtained. :

[0245]

[0246] In the formula:

[0247] M — a 19×19 orthogonal transformation matrix;

[0248] Step 5-2:

[0249] Based on the macroscopic density calculated in step 3 Calculate the corresponding equilibrium moment vector based on the velocity u calculated in step 4. These moment vectors are analytical expressions of macroscopic physical quantities in moment space, representing the local equilibrium state that the system tends towards.

[0250] Step 5-3: Calculate the moment of force term:

[0251] To incorporate macroscopic forces (including interphase forces, gravity, etc.) into the evolution equations of the lattice Boltzmann method, it is necessary to transform the forces into moment space and calculate the moment of the force terms. ;

[0252] This force term moment originates from a higher-order analysis of the evolution equations of the lattice Boltzmann method, and aims to compensate for the errors introduced by discrete effects and force terms;

[0253] Step 5-4: Calculate the collision in the moment space:

[0254]

[0255] In the formula:

[0256] —The moment vector after the collision;

[0257] —A regulating term that enhances mechanical stability;

[0258] S – Diagonal relaxation matrix;

[0259]

[0260] Step 5-5: Inverse transformation from moment space back to velocity space:

[0261]

[0262] Step 6: Calculate the migration process:

[0263]

[0264] In the formula:

[0265] t —The current moment;

[0266] Step 7: Handle boundary conditions:

[0267] Step 7-1: Bottom boundary of the fluid domain, i.e., fluid-structure interaction boundary:

[0268] Processed using standard half-step bounce format Boundary layer to achieve no-slip boundary conditions on solid walls:

[0269]

[0270] In the formula:

[0271] i —Index of discrete velocity directions;

[0272] opp( i — Index of the discrete velocity direction (this velocity direction is opposite to the direction indexed i);

[0273] Step 7-2: Top boundary of the fluid domain:

[0274] Pressure boundary conditions are achieved using a non-equilibrium extrapolation scheme: This method sets the target density on the virtual boundary. And assume that the non-equilibrium distribution function at the boundary is similar to that of the adjacent fluid internal grid points. That is, k = Nz-2 is the same, so the unknown distribution function flowing into the computational domain is determined. :

[0275]

[0276] In the formula:

[0277] —Boundary grid point at k = Nz-1;

[0278] —The density of the grid points at the top k = Nz-1;

[0279] —Boundary grid point at k = Nz-2;

[0280] —The velocity at grid point Nz-2;

[0281] —Density at Nz-2 lattice points;

[0282] Step 8: Solve for the temperature field:

[0283] Calculate the temperature distribution function:

[0284] Step 8-1: Calculate the relaxation time of the temperature distribution function:

[0285] Relaxation time of the fluid domain ;

[0286] relaxation time of solid domain ;

[0287] Step 8-2: Calculate the collision process of the temperature distribution function using the source-term modified lattice Boltzmann equation:

[0288]

[0289] In the formula:

[0290] —Temperature source item;

[0291] Step 8-3: Calculate the migration process:

[0292] The migration process of the temperature distribution function will involve the output distribution function after collision relaxation. Along discrete velocity Direction, migrate to its adjacent grid points At this point, it becomes the input distribution function for the next time step. For the x and y directions, due to the use of periodic boundary conditions, the migration is cyclical; for the z direction, when migrating to the boundary or fluid-structure interface, its direction is unknown. The value will be handled in the boundary conditions step.

[0293] The migration equation is as follows:

[0294]

[0295] Step 8-4: Update macroscopic temperature parameters:

[0296] Update the temperature of the entire fluid domain :

[0297]

[0298] Step 8-5: Calculate the fluid-structure conjugate heat transfer boundary conditions:

[0299] On the solid side of the fluid-solid interface k s =Z solid -1 grid points fluid side of the fluid-solid interface k f =Z solid Grid points A non-equilibrium extrapolation method is used to achieve fluid-solid conjugate heat transfer, ensuring the continuity of temperature and heat flux density at the interface:

[0300] Fluid-structure interaction interface temperature ;

[0301] Fluid-side temperature field input distribution function ,in It is the internal lattice point k=Z of the fluid. solid +1;

[0302] Solid side: ,in The internal lattice point of the solid is k=Z solid -2;

[0303] External boundary: bottom k=0, isothermal T wb The bounce boundary scheme is adopted; the top k=Nz-1 isothermal. T s All of these are implemented using an unbalanced extrapolation scheme.

[0304] Step 8-6: Solve the heat conduction equation in the solid domain using the finite difference method:

[0305] Heat conduction equation:

[0306]

[0307] In the formula:

[0308] —The thermal conductivity of a solid;

[0309] t —The current moment;

[0310] Discretizing the heat conduction equation using the explicit finite difference method:

[0311]

[0312] In the formula:

[0313] ——Current moment ( i , j , k Temperature of the grid points;

[0314] —The previous moment ( i , j , k The temperature of the grid points, ;

[0315] in, = =1, i , j , k These represent the grid point labels in the x, y, and z directions, respectively.

[0316] Step 9: Repeat steps 3 through 8 until the set total number of simulation steps is reached.

[0317] The change in droplet mass over time can be obtained through programming calculations using the above method, such as... Figure 4 As can be seen from the figure, the mass of the droplet is constantly decreasing, indicating that the droplet is gradually being evaporated;

[0318] By performing programming calculations using the above method, the droplet positions at different times can be obtained, as shown below. Figure 5 As shown, at step 0, the droplet is located at the center of the computational domain. At step 88, the droplet impacts the high-temperature chip. Because the chip temperature is higher than the Leidenfrost temperature point, the Leidenfrost effect occurs, creating a uniform layer between the bottom of the droplet and the chip. At step 151, the droplet jumps off the high-temperature chip surface. At step 281, the droplet reaches its maximum bounce height and begins to fall due to gravity. At step 507, the droplet impacts the high-temperature chip surface again. The chip heats the droplet, creating a vapor layer at its bottom, which again causes the droplet to leave the chip surface. Figure 5 It can be seen that this method can simulate the bouncing phenomenon of droplets when the Leidenfrost effect occurs;

[0319] By performing programming calculations using the above method, the chip temperature change along the y-direction at z=18 and x=75 is obtained, as follows: Figure 6 The temperature in the middle is lower than the temperatures on both sides, because the middle area is cooled by the droplets; by Figure 6 It can be seen that this method can successfully simulate the cooling effect of droplets on the chip, that is, realize fluid-structure conjugate heat transfer.

[0320] The above description is a further detailed explanation of the present invention in conjunction with specific preferred embodiments. It should not be considered that the specific embodiments of the present invention are limited to this. For those skilled in the art, several simple deductions or substitutions can be made without departing from the concept of the present invention, and all such deductions or substitutions should be considered to fall within the scope of patent protection determined by the submitted claims.

Claims

1. A three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method, characterized in that: Includes the following steps: Step 1: Set basic parameters: Step 1-1: Set initial physical property parameters: liquid phase density is set to... The gas phase density is set to The specific heat at constant volume is set to The liquid phase diffusion coefficient is and gas phase thermal diffusivity This is used for subsequent calculations of the fluid's thermal diffusivity; the thickness Z of the solid wall is set. solid Thermal diffusivity of solid walls Density of solids relaxation time of solids ; Steps 1-2: Set simulation parameters: Use the D3Q19 model as the discrete model for velocity and temperature. This model contains 19 discrete velocity directions. The flow field uses a multi-relaxation time collision operator with higher numerical stability, and the temperature field uses a single-relaxation time collision operator. Construct a three-dimensional computational domain of size Nx×Ny×Nz, i.e., Nx grid points in the x-direction, Ny grid points in the y-direction, and Nz grid points in the z-direction, for a total of Nx×Ny×Nz grid points. Set the total number of simulation steps N. step ; Step 2: Initialize the grid type, fluid domain density field, velocity field, temperature field, and distribution function: Step 3: Calculate the fluid domain density, pressure, and thermal diffusivity: Step 4: Calculate the resultant force of intermolecular interaction forces and volume forces between different phases: Step 5: Computation of the fluid domain collision process: Step 6: Calculate the migration process: Step 7: Handle wall boundary conditions: Step 7-1: The bottom boundary of the fluid domain, i.e., the fluid-solid conjugate heat transfer boundary, is handled using the standard half-step bounce scheme to achieve a no-slip boundary condition on the solid wall: In the formula: --exist At time, located at the boundary fluid node At that point, along with discrete velocity i The input distribution function in the opposite direction, where, t For the current moment, For time step; --exist At time, located at the boundary fluid node At that point, along the discrete velocity i The distribution function of the directional output; Step 7-2: Top boundary of the fluid domain. Pressure boundary conditions are achieved using a non-equilibrium extrapolation scheme. This method involves setting a target density on the virtual boundary. And assuming that the non-equilibrium distribution function at the boundary is the same as that of the adjacent fluid internal grid points, the unknown distribution function flowing into the computational domain is determined. : In the formula: ——Grid Input distribution function at grid point It is a boundary lattice point at k = Nz-1; --use and The calculated equilibrium distribution function uses the pressure equilibrium boundary in this method, assuming the velocity gradient at the upper boundary is 0, and uses the velocities of adjacent internal fluid grid points. To approximate, It is the density at lattice point k=Nz-1; —Inner fluid grid points immediately adjacent to the upper boundary Input distribution function at grid point It is a boundary lattice point at k = Nz-2; --use and The calculated equilibrium distribution function, The lattice point represents the density at k = Nz - 2; Step 8: Solve for the temperature field: Calculate the temperature distribution function : Step 8-1: Calculate the relaxation time of the temperature distribution function: Relaxation time of the fluid domain , The fluid thermal diffusivity; Step 8-2: Calculate the temperature distribution function using the source term-corrected lattice Boltzmann equation. The collision process: In the formula: —Temperature source item; Discrete direction i Temperature field output distribution function; Discrete direction i Temperature field input distribution function; —Temperature is T When, discrete direction i The temperature field equilibrium distribution function; Step 8-3: Calculate the migration process: The migration process of the temperature distribution function will involve the output distribution function after collision relaxation. Along discrete velocity Direction, migrate to its adjacent grid points At this point, it becomes the input distribution function for the next time step. For the x and y directions, due to the use of periodic boundary conditions, the migration is cyclical; for the z direction, when migrating to the boundary or fluid-structure interface, its direction is unknown. The value will be handled in the boundary conditions step; The migration equation is as follows: In the formula: —— time, Location, discrete direction i Temperature field input distribution function; —— t time, Location, discrete direction i Temperature field output distribution function; Step 8-4: Update macroscopic temperature parameters: Update the temperature of the entire fluid domain : In the formula: Discrete direction i Temperature field distribution function; Step 8-5: Calculate the fluid-structure conjugate heat transfer boundary conditions: On the solid side of the fluid-solid interface k s =Z solid -1 grid points fluid side of the fluid-solid interface k f =Z solid Grid points A non-equilibrium extrapolation method is used to achieve fluid-solid conjugate heat transfer, ensuring the continuity of temperature and heat flux density at the interface: Fluid-structure interaction interface temperature ,in, It is a grid. temperature, It is a grid. Temperature; Fluid-side temperature field input distribution function ,in It is the internal lattice point k=Z of the fluid. solid +1, For grid points Discrete direction i The corresponding temperature field input distribution function, Discrete directions at the fluid-structure interaction interface i The corresponding temperature field equilibrium distribution function, For grid points Discrete direction i The corresponding temperature field input distribution function, For grid points Discrete direction i The corresponding temperature field equilibrium distribution function; Solid side: ,in The internal lattice point of the solid is k=Z solid -2, For grid points Discrete direction i The corresponding temperature field input distribution function, For grid points Discrete direction i The corresponding temperature field input distribution function, For grid points Discrete direction i The corresponding temperature field equilibrium distribution function; External boundary: bottom k=0, isothermal T wb The bounce boundary scheme is adopted; the top k=Nz-1 isothermal. T s All of these are implemented using an unbalanced extrapolation scheme. Step 8-6: Solve the heat conduction equation in the solid domain using the finite difference method: Heat conduction equation: —The thermal conductivity of a solid; T —Grid point temperature; Discretizing the heat conduction equation using the explicit finite difference method: In the formula: ——Current moment ( i , j , k Temperature of the grid points; —The previous moment ( i , j , k The temperature of the grid points, ; in, = =1, i , j , k These represent the grid point labels in the x, y, and z directions, respectively. Step 9: Repeat steps 3 through 8 until the set total number of simulation steps is reached.

2. The three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method according to claim 1, characterized in that: Step 2 is as follows: Step 2-1: Initialize grid type: The grid points constructed in steps 1-2 are classified using the grid identifier pointsflag, dividing the entire computational domain into internal fluid domain grid points, fluid domain boundary grid points, solid domain grid points, and virtual fluid domain grid points: Step 2-2: Initialize the fluid domain density field: The droplet radius is set to R0, and the droplet center position is set to (x0, y0, z0); The gas-liquid interface width is set to ; Calculate the distance between the grid point (x, y, z) and the droplet center position (x0, y0, z0). d ; The initial grid density of the fluid domain is calculated using the following formula. : The densities of solid virtual grid points and boundary grid points are not included in the calculation, so they do not need to be set. Step 2-3: Initialize the velocity field: Determine the distance between the grid point (x, y, z) and the center of the droplet (x0, y0, z0). d Is the initial velocity of the lattice point (x, y, z) less than or equal to the droplet radius R0? If so, then... Set to (0, 0, -0.2), otherwise set to 0; Step 2-4: Initialize the temperature field: Set the temperature of the bottom solid domain grid point to T b Initialize the grid point temperature within the fluid domain to T s =0.86 T c ,in T c It is the critical temperature; Step 2-5: Initialize the distribution function: Initialize the temperature field equilibrium distribution function across the entire computational domain. Temperature field input distribution function and temperature field output distribution function ; Initial velocity field equilibrium distribution function in the fluid domain Velocity field input distribution function and velocity field output distribution function ; Initial Moment Space Equilibrium Distribution Function in Fluid Domain Moment space input distribution function Sum-moment space output distribution function .

3. The three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method according to claim 1, characterized in that: Step 3 is as follows: Step 3-1: Calculate the density of the fluid domain : Step 3-2: Calculate pressure using the equation of state The following formula gives the calculation formula for the PR state equation, or it can be replaced with other state equations for calculation: In the formula; R —Gas constant; T --temperature; —Constant, take ; b --constant, ; —Critical pressure; K — Scaling factor, used to adjust the interface thickness during simulation; Step 3-3: Calculate the thermal diffusivity: Within the fluid domain, the thermal diffusivity is modeled using a density-dependent linear interpolation model, and based on... and diffusion coefficient at time and Dynamic calculation: For fluid thermal diffusivity, in the solid domain, the solid thermal diffusivity is... Maintain its preset fixed value.

4. The three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method as described in claim 1, characterized in that: Step 4 is as follows: Step 4-1: Calculate the potential function of the fluid domain : In the formula: G —Inter-particle interaction coefficient, G =-1; c —Speed ​​of sound; —Grid speed of sound, ; Step 4-2: Calculate the intermolecular interaction forces between different phases : In the formula: —The i-th discrete velocity vector in the D3Q19 model; —D3Q19 model The corresponding weighting coefficients; —Time step; —Current grid point position; Step 4-3: Calculate volume forces : In the formula: —The average density of the fluid domain; g — acceleration due to gravity; Step 4-4: Calculate the resultant force F: Steps 4-5: Calculation speed : 。 5. The three-dimensional numerical simulation method for the Leiden-Frost effect based on the dual-distribution function lattice Boltzmann method coupled with the finite difference method according to claim 1, characterized in that: Step 5 is as follows: Step 5-1: Obtain the input distribution function of the velocity field from initialization step 2. By projecting the vector onto the moment space through an orthogonal transformation matrix M, the moment vector before the collision is obtained. : In the formula: M — a 19×19 orthogonal transformation matrix; Step 5-2: Based on the density calculated in step 3 Calculate the corresponding equilibrium moment vector based on the velocity u calculated in step 4. These moment vectors are analytical expressions of macroscopic physical quantities in moment space, representing the local equilibrium state that the system tends to. Step 5-3: Calculate the moment of force term: To incorporate macroscopic forces into the evolution equations of the lattice Boltzmann method, it is necessary to transform the forces into moment space and calculate the moment of the force terms. ; This force term moment originates from a higher-order analysis of the evolution equations of the lattice Boltzmann method, and aims to compensate for the errors introduced by discrete effects and force terms; Step 5-4: Calculate the collision in the moment space: In the formula: S – Diagonal relaxation matrix; —A regulating term that enhances mechanical stability; S – Diagonal relaxation matrix; Step 5-5: Inverse transformation from moment space back to velocity space: 。

Images