Three-phase flow and heat transfer coupling simulation method and system based on IBM-voset-dem

Abstract

Provided in the present invention are a three-phase flow and heat transfer coupling simulation method and system based on an IBM-VOSET-DEM. The method comprises: determining a computational domain and dividing the grid of the computational domain; initializing parameters; reconstructing a liquid-liquid interface; determining a fluid volume fraction at an (n+1)th moment; determining a boundary grid, and determining a signed distance function of a grid cell at an nth moment; determining an interface curvature of the grid cell of the interface; smoothing abruptly changing physical quantities of the interface; determining the surface tension of the interface; determining a fluid intermediate velocity; obtaining fluid volume fractions and fluid intermediate velocities of all of the grid cells in the computational domain; determining the velocity of a particle; determining the displacement of a pth solid particle at the (n+1)th moment and a solid volume fraction of the grid cell at the (n+1)th moment; determining the final velocity of the computational domain at the (n+1)th moment; and solving for a temperature and an interface heat flux density. In the present invention, a DEM soft-sphere model is used to calculate a collision force of a solid, such that microscopic motion information of a solid particle can be accurately obtained and tracked; and an interface heat flux density model is used to take the conjugate heat transfer effect of an interface into consideration.

Description

A three-phase flow and heat transfer coupling simulation method and system based on IBM-VOSET-DEM Technical Field

[0001] The present invention relates to the field of fluid mechanics simulation, and in particular to a three-phase flow and heat transfer coupling simulation method and system based on IBM-VOSET-DEM. Background Art

[0002] Flow and heat transfer problems involving interface changes in two-phase flows containing solid particles are widespread in the energy and chemical industries. These three-phase flow and heat transfer problems involve complex interactions among liquid-liquid, solid-liquid, and solid-solid. Experimental studies are severely limited in their ability to obtain microscopic information about multiphase flow and heat transfer. Numerical simulations address this deficiency and have become a crucial tool for studying multiphase flow and heat transfer.

[0003] Regarding numerical simulations of multiphase flow and heat transfer, for liquid-liquid two-phase flow, the volume of fluid (VOF), level set (LS), and their coupled models (CLSVOF, VOSET) are commonly used to capture the two-phase interface. For particle multiphase flow models, particles are typically treated as a continuous phase, but this does not accurately capture their microscopic motion and their effects on the fluid. Alternatively, particles are treated as discrete phases, using methods such as DPM and the discrete element method (DEM) to obtain information on particle motion.

[0004] However, for three-phase systems, there are currently few flow and heat transfer mathematical models that simultaneously consider liquid-liquid interaction, solid-liquid interaction, and solid-solid interaction.

[0005] Summary of the Invention

[0006] In response to the shortcomings of the prior art, the present invention provides a three-phase flow and heat transfer coupled simulation method and system based on IBM-VOSET-DEM. The VOSET method is used to capture the liquid-liquid interface, overcoming the disadvantage that the Level Set method is very likely to cause fluid mass non-conservation in the calculation. At the same time, it overcomes the disadvantages of the VOF method in calculating the large error of the interface curvature and the poor effect of sudden changes in physical quantities near the smooth interface. The immersed boundary method IBM (Immersed boundary method) is used to obtain the solid particle interface, avoiding the generation of body-fitting grids at the particle boundary, thereby reducing the calculation cost. The DEM soft ball model is used to calculate the solid collision force, which can accurately obtain and track the motion microscopic information of the solid particles. The interface heat flux density model is used to consider the conjugate heat transfer effect of the liquid-solid interface.

[0007] The present invention achieves the above technical objectives through the following technical means.

[0008] A three-phase flow and heat transfer coupled simulation method based on IBM-VOSET-DEM includes the following steps:

[0009] S01: Determine the computational domain and the grids that divide the computational domain: Determine the computational domain size, the number of particles N and the particle size D in the computational domain; any grid in the computational domain is represented by (N x ,N y ), where x∈[1, N1]; y∈[1, N2]; set the calculation time step Δt, the total number of time steps N t ; Where N1 is the total number of grids in the horizontal direction; N2 is the total number of grids in the vertical direction;

[0010] S02: Initialization parameters: Determine the velocity of the fluid in all grid cells at the initial moment The temperature of the fluid in all grid cells at the initial moment and the volume fraction F of the fluid in all grid cells at the initial moment 0 (x,y); determine the center position of all solid particles at the initial moment (x s ,y s ), the velocity of all grid unit solids at the initial moment The temperature of the solid in all grid cells at the initial moment and the volume fraction of solid in all grid cells at the initial moment Determine all thermophysical parameters, including: density of the first fluid ρ1, density of the second fluid ρ2, density of solid particles ρ s , viscosity of the first fluid μ1, viscosity of the second fluid μ2, thermal conductivity of the first fluid λ1, thermal conductivity of the second fluid λ2, thermal conductivity of the solid particles λ s ;

[0011] S03: Reconstruction of liquid-liquid interface:

[0012] When the fluid volume fraction F n (x,y)∈(0,1), then the network unit (N x ,N y ) is the liquid-liquid interface grid; when F n When (x,y)=0, the network unit (N x ,N y ) is filled with the second fluid; when F n When (x,y)=1, the network unit (N x ,N y ) is filled with a first fluid;

[0013] When the grid (x, y) is a liquid-liquid interface grid, the piecewise linear interpolation technique PLIC is used to reconstruct the interface shape of all grid cells at the nth moment;

[0014] S04: Determine the fluid volume fraction at time n+1:

[0015] According to the grid unit (N x ,N y ) and the speed of the grid unit at the nth moment (N x ,N y ) fluid volume fraction F n (x,y) to obtain the inflow and outflow of the grid cell (N x ,N y ) of the fluid volume, and then obtain the fluid volume fraction F of the network unit at time n+1 n+1 (x,y);

[0016] S05: Determine the boundary grid. When the grid (x, y) is a liquid-liquid interface grid, determine the grid unit (N) at the nth moment. x ,N y )'s signed distance function

[0017] S06: Determine the grid unit of the interface (N x ,N y )’s interface curvature k(x,y):

[0018] S07: Physical quantity of smooth interface mutation: Use Heaviside function to determine the grid unit (N at the nth moment) x ,N y )Density and viscosity near the smooth interface:

[0019] The grid unit at the nth moment (N x ,N y )The smooth density calculation formula is:

[0020] The grid unit at the nth moment (N x ,N y ) The formula for calculating smooth viscosity is:

[0021] Smooth Heaviside function for:

[0022] Where, ∈ is the thickness of the smooth interface;

[0023] S08: Determine the surface tension of an interface

[0024] Solve the surface tension of the liquid-liquid interface, the grid unit (N x ,N y )Surface tension F s n (x,y) is calculated using the continuous surface force model. The specific calculation formula is as follows:

[0025] Where σ is the surface tension coefficient,

[0026] When there is no interface, there is no surface tension, and the grid unit (N x ,N y ) surface tension F s n (x,y)=0;

[0027] S09: Determine the intermediate velocity of the fluid

[0028] S10: Repeat steps S03-S9 to obtain the fluid volume fraction F of all grid cells in the calculation domain n (x,y) and intermediate speeds

[0029] S11: Determine particle velocity:

[0030] The IBM and DEM methods are used to determine the fluid forces acting on solid particles. and collision force

[0031] Determine the velocity of the particle: Based on the total fluid force and total collision force on the pth particle, use Newton's momentum equation and angular momentum equation to solve the translational velocity of the pth particle. The rotation speed of the pth particle and the velocity vector of the pth particle The calculation formula is as follows:

[0032] Among them, m sp is the mass of the pth particle, I sp is the inertia tensor of the pth particle, is the radius vector from the center of the pth particle to the boundary of the pth particle; is the total fluid force on the pth particle; is the total collision force on the pth particle at the nth moment

[0033] Obtain the velocity of all particles, and then obtain the velocity vectors of all particles in the computational domain, and the fluid forces of all grid cells occupied by the particles

[0034] S12, determine the displacement d of the pth solid particle at the n+1th moment n+1 and the grid unit at the n+1th moment (N x ,N y ) solid volume fraction

[0035] S13: Determine the final velocity u of the computational domain at the n+1th moment n+1 (x,y);

[0036] S14: Solve the temperature and interface heat flux density. The specific steps are as follows:

[0037] A. Solve the particle interface heat flux density:

[0038] A-1. Based on the solid volume fraction at the n+1th moment Get all solid interface meshes and set all solid interface meshes as a collection Among them, H g For collection The gth interface grid unit in the set G is The total number of interface grid cells within;

[0039] A-2. On the interface grid H g The temperature gradient is divided into the normal direction and tangent direction Where i and j are the grid numbers in the horizontal and vertical directions respectively, T n (H g ) is the interface grid H at the nth moment g temperature, is the unit vector in the normal direction of the particle interface, and I is the unit matrix;

[0040] A-3. The thermal conductivity in the normal and tangential directions is calculated using the following formula:

[0041] In the formula is the interface grid unit H e Thermal conductivity in the normal direction at the n+1th moment, is the tangential thermal conductivity at the n+1th moment; λ s is the thermal conductivity of the solid; λ f The thermal conductivity of the fluid, when the solid interface F n+1 (H g )>0, then λ f The value is λ f1 , otherwise the value is λ f2 ;

[0042] A-4. Heat flux at the interface grid Calculated by the following formula:

[0043] B. Solution of fluid interface heat flux density:

[0044] B-1. According to the fluid volume fraction F at the n+1th moment n+1 (x,y) determines all liquid-liquid interface meshes; all liquid-liquid interface meshes are set as a collection Among them, F h For the collection The hth interface grid unit in, H is the set The total number of interface grid cells within;

[0045] B-2. On the interface grid F h The temperature gradient is divided into the normal direction and tangent direction Where i and j are the grid numbers in the horizontal and vertical directions respectively, and T is T n (F h ) is the interface grid F at the nth moment h The temperature, is the unit vector in the normal direction of the liquid-liquid interface, and I is the unit matrix;

[0046] B-3. ​​The thermal conductivity in the normal and tangential directions is calculated using the following formula:

[0047] In the formula is the interface grid unit F h Thermal conductivity in the normal direction at the n+1th moment, is the tangential thermal conductivity at the n+1th moment; λ f1 and λ f2 are the thermal conductivities of the first fluid and the second fluid, respectively;

[0048] B-4. Heat flux at the interface mesh Calculated by the following formula:

[0049] C. Use the Crank-Nicolson method to update the diffusion term of the energy equation in time and solve the temperature field on the Euler system;

[0050] S15, determine whether the calculation is finished: if the current time step number n <N t , then repeat steps S03-S14; if the current time step number n≥N t , the calculation terminates.

[0051] Furthermore, the piecewise linear interpolation technique PLIC is used to reconstruct the interface shapes of all grid cells at the nth moment, specifically:

[0052] According to the grid unit (N x ,N y ) of the fluid volume fraction F n (x,y) determines the unit normal vector of the interface of the grid cell at the nth moment in, is the Hamiltonian operator, i represents the horizontal component, j represents the vertical component; where n is the number of time steps, n∈[1, N t ];

[0053] According to the grid unit (N x ,N y ) of the grid unit normal vector and the grid unit at the nth moment (N x ,N y ) of the fluid volume fraction F n (x,y) determines the grid cell (N x ,N y )The interface shape within the grid.

[0054] Further, determine the boundary grid and solve the signed distance function Specifically:

[0055] At the nth moment, the fluid volume fraction F of the grid unit n When (x,y) is 0.5, the grid cell is a boundary grid; determine all boundary grids of the computational domain;

[0056] Computational grid cells (N x ,N y ) to each boundary grid, thus obtaining the grid cell (N x ,N y ) to the minimum distance d among all boundary grids, that is, the grid unit (N x ,N y ) the shortest distance d from the two-phase interface;

[0057] According to the grid unit (N x ,N y ) of the fluid volume fraction F n (x,y) determines the grid cell (N x ,N y )'s signed distance function

[0058] Furthermore, the mesh unit of the interface (N x ,N y )’s interface curvature k(x,y), specifically:

[0059] When the liquid-liquid interface appears, according to the grid unit at the nth moment (N x ,N y )'s signed distance function Determine the interface curvature k(x,y) using the following calculation formula:

[0060] When the grid cell (N x ,N y ) The fluid volume fraction F n (x,y)∈(0,1) and solid volume fraction When the grid unit (N x ,N y ) is a three-phase contact unit; when three phases are in contact, the contact angle boundary condition is used on the contact line to correct the interface unit normal vector of the mesh unit. The correction equation is:

[0061] Where θ is the contact angle, is the grid unit (N x ,N y ) solid wall normal vector, is the grid unit (N x ,N y ) solid wall tangent vector;

[0062] Determine the grid cell (N x ,N y )’s interface curvature k(x,y):

[0063] Further, determine the fluid intermediate velocity Specifically:

[0064] According to the grid unit (N x ,N y ) to solve the Navier-Stokes equations by using the velocity, pressure and surface tension inside the grid and obtain the grid cell (N x ,N y ) of the intermediate speed The specific calculation formula is as follows:

[0065] Where Δt is the time step, is the convection term, is a sticky item;

[0066] Based on the fast Fourier transform method, the pressure Poisson equation is solved to obtain the grid unit (N at the n+1th moment) x ,N y ) pressure correction value P n+1 (x,y);

[0067] Using the grid unit (N at the n+1th moment x ,N y ) pressure correction value P n+1 (x,y) corrects the grid cell (N) at the n+1th moment x ,N y ) of the intermediate speed Get the corrected intermediate speed Specifically:

[0068] Furthermore, IBM and DEM methods were used to determine the fluid forces acting on the solid particles. and collision force Specifically:

[0069] Determining the fluid forces acting on solid particles using the IBM method

[0070] The fluid force on the pth particle at the nth moment is determined by vector summing the fluid forces on the grid occupied by the particle; based on the center position and radius of the pth particle at the nth moment, the set of all grid cells occupied by the pth particle is determined Among them, P e For the collection The e-th grid cell in the set, E is the set The total number of grid cells within;

[0071] Solve the set Q according to the IBM method p Medium P e Fluid forces acting on the mesh This method can also make the solid-liquid interface meet the no-slip boundary condition. The specific calculation formula is as follows:

[0072] in:

[0073] P is the n+1th moment e The median speed of the grid; P is the nth moment e Solid velocity of the grid; P is the nth moment e The solid particle volume fraction of the grid;

[0074] Get a collection The fluid force on all grid cells at the nth moment is the fluid force on all grid cells occupied by the pth particle at the nth moment. The total fluid force on the pth particle is obtained by summing them up.

[0075] in, is the total fluid force on the pth particle;

[0076] Determining the impact force on solid particles using the DEM method

[0077] Determine the set of particles except the pth particle in Q q For collection The qth particle in the calculation domain, N is the number of particles in the calculation domain;

[0078] Determine the collection of walls Where w represents different wall surfaces, w = 1, 2, 3, 4;

[0079] The DEM soft ball collision model is used to solve the collision force between the pth particle and the qth particle at the nth moment. and the collision force between the pth particle and the wth wall at the nth moment

[0080] When the distance between the center of the pth particle and the center of the qth particle is greater than D, the pth particle and the qth particle do not collide, and the collision force between the pth particle and the qth particle at the nth moment is That is 0;

[0081] When the distance between the center of the pth particle and the w wall is greater than D / 2, the pth particle does not collide with the w wall, and the collision force between the pth particle and the w wall at the nth moment is That is 0;

[0082] When the distance between the center of the pth particle and the center of the qth particle is less than D, the pth particle collides with the qth particle, and the collision force between the pth particle and the qth particle is The calculation formula is as follows:

[0083] Where: ε is the elastic coefficient, η is the viscous dissipation coefficient, and are the displacements in the tangential and normal directions, is the relative velocity between the pth and qth particles, and are the tangential and normal velocities of the pth particle, is the unit normal vector from the center of the pth particle to the center of the qth particle;

[0084] When the distance between the center of the pth particle and the w wall is less than D / 2, the pth particle collides with the w wall, and the collision force between the pth particle and the w wall is The calculation formula is as follows:

[0085] in: is the relative velocity between the pth particle and the wth wall, is the unit normal vector pointing from the center of the pth particle to the wth wall;

[0086] The total collision force on particle p at time n is the sum of the forces of collision with other particles or walls:

[0087] Furthermore, the displacement d of the pth solid particle at the n+1th moment is determined n+1 , specifically:

[0088] The velocity vector of the pth particle Decompose into horizontal and vertical velocity components and According to the pth particle speed and Solve the horizontal and vertical displacements of the pth particle in the next time step Δt respectively, and then obtain the position of the pth particle at the n+1th moment:

[0089] in, and are the horizontal and vertical coordinates of the pth particle at the nth moment, and is the component velocity of the pth particle in the horizontal and vertical directions at the nth moment, is the displacement of the pth particle from the nth moment to the n+1th moment.

[0090] Further, determine the grid unit (N x ,N y ) solid volume fraction Specifically:

[0091] When the grid cell (N x ,N y ) from the center of the pth particle When the shortest distance is less than the particle radius D / 2, the grid unit (N x ,N y ) solid particle volume fraction

[0092] When the grid cell (N x ,N y ) from the center of the pth particle When the maximum distance is greater than the particle radius, the grid unit (N x ,N y ) solid particle volume fraction

[0093] When the grid cell (N x ,N y ) from the center of the pth particle When the shortest distance is greater than D / 2 and the longest distance is less than D / 2, then the grid unit (N x ,N y ) is the solid-liquid boundary, and its volume fraction is calculated as follows: γ=0.05(1-β 2 )+0.3

[0094] Where h is the distance from the solid-liquid interface to the grid unit (N x ,N y ) is the distance between the centers; Δx is the grid size; Solid-liquid interface normal vector The component in the horizontal direction, Solid-liquid interface normal vector The component in the vertical direction,

[0095] According to the displacement of all particles and the solid volume fraction of all particles Determine the grid cell (N) at the n+1th moment in the computational domain. x ,N y ) solid volume fraction

[0096] if but

[0097] Furthermore, determine the final velocity u of the computational domain at the n+1th moment n+1 (x,y), specifically:

[0098] The grid cell at the nth moment (N x ,N y The force exerted on the fluid by the particle is according to Get the fluid velocity at time n+1:

[0099] The final velocity u of the entire computational domain at the n+1th moment is determined using the following formula: n+1 (x,y):

[0100] A three-phase flow heat transfer coupling simulation system based on IBM-VOSET-DEM comprises a storage medium; the storage medium stores a program written using the three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM.

[0101] The beneficial effects of the present invention are:

[0102] 1. The IBM-VOSET-DEM-based three-phase flow and heat transfer coupled simulation method and system described in this invention can simultaneously capture the three-phase interfaces, obtain momentum exchange between fluid-fluid, fluid-solid, and solid-solid, and obtain the heat flux density near the interfaces, thereby more accurately obtaining the flow and heat transfer phenomena within the three-phase system.

[0103] 2. The IBM-VOSET-DEM-based three-phase flow and heat transfer coupled simulation method and system described in this invention uses the VOSET method to capture the liquid-liquid interface, overcoming the shortcomings of the Level Set method, which is prone to fluid mass non-conservation during calculations. It also overcomes the shortcomings of the VOF method, which suffers from large errors in calculating interface curvature and poor results in smoothing sudden changes in physical quantities near the interface.

[0104] 3. The IBM-VOSET-DEM-based three-phase flow and heat transfer coupled simulation method and system described in this invention uses the IBM method of marking the fluid-solid interface with a series of Lagrange points, avoiding the body-fitting mesh generation at particle boundaries required in traditional methods, thereby reducing computational costs. Furthermore, by adding a force source term to the momentum equation, it not only solves for momentum exchange between the fluid and solid, but also satisfies the no-slip boundary condition at the fluid-solid interface.

[0105] 4. The IBM-VOSET-DEM-based three-phase flow and heat transfer coupled simulation method and system described in the present invention uses the DEM soft-sphere model to calculate solid collision forces, enabling accurate acquisition and tracking of solid particle motion information.

[0106] 5. The IBM-VOSET-DEM-based three-phase flow heat transfer coupled simulation method and system described in the present invention adopts an interfacial heat flux model to consider the conjugate heat transfer effects of the solid-liquid and liquid-liquid interfaces. This model considers the interface direction and solves the temperature field in an Euler system, which is highly compatible with the immersed solid method. This calculation method is more consistent with actual heat transfer and can more accurately simulate heat transfer at the boundary. BRIEF DESCRIPTION OF THE DRAWINGS

[0107] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. The drawings described below are some embodiments of the present invention. For ordinary technicians in this field, it is obvious that other drawings can be obtained based on these drawings without paying any creative work.

[0108] FIG1 is a flow chart of a three-phase flow and heat transfer coupled simulation based on IBM-VOSET-DEM according to the present invention;

[0109] FIG2 is a calculation flow of the VOSET method of the present invention;

[0110] FIG3 is a schematic diagram of a calculation domain according to an embodiment of the present invention;

[0111] FIG4 is a schematic diagram of a geometric calculation method for calculating a signed distance function according to the present invention.

[0112] FIG5 is a schematic diagram of the contact angle boundary conditions used in the present invention;

[0113] FIG6 is a flow chart of fluid motion calculations according to the present invention;

[0114] FIG7 is a solid motion calculation process of the present invention;

[0115] FIG8 is a schematic diagram of the calculation of thermal conductivity of the solid-liquid interface grid in the normal and tangential directions of the present invention;

[0116] FIG9 is a diagram showing simulation results of a liquid droplet impacting a spherical solid particle according to an embodiment of the present invention;

[0117] FIG10 is a comparison of the ratio of the minimum height of a droplet to the initial radius of the droplet after the droplet impacts a spherical solid particle in the embodiment of the present invention and in the experiment;

[0118] FIG11 is a graph showing the change in volume over time during the evaporation process of a liquid droplet on a hot spherical solid particle in an embodiment of the present invention and an experiment. DETAILED DESCRIPTION

[0119] The following describes embodiments of the present invention in detail, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and are not to be construed as limiting the present invention.

[0120] In the description of the present invention, it should be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "up", "down", "axial", "radial", "vertical", "horizontal", "inside", "outside" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the present invention. In addition, the terms "first" and "second" are used for descriptive purposes only, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined as "first" and "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, "multiple" means two or more, unless otherwise clearly and specifically defined.

[0121] In the present invention, unless otherwise expressly specified or limited, the terms "mounted," "connected," "connect," "fixed," etc. should be understood broadly. For example, they may refer to fixed, detachable, or integral connections; mechanical or electrical connections; direct or indirect connections through an intermediary; or internal communication between two components. Those skilled in the art will understand the specific meanings of the above terms in the present invention based on specific circumstances.

[0122] As shown in FIG1 , the three-phase flow heat transfer coupled simulation method based on IBM-VOSET-DEM of the present invention includes the following steps:

[0123] S01. Determine the computational domain and the grids for dividing the computational domain: Determine the computational domain size L*H (length*height), determine the number of particles N and particle size D in the computational domain; divide the computational domain into grids (N1, N2) under the Euler system, where N1 is the total number of grids in the horizontal direction; N2 is the total number of grids in the vertical direction; any grid in the computational domain is represented by (N x ,N y ), where x∈[1, N1]; y∈[1, N2]; set the calculation time step Δt, the total number of time steps N t .

[0124] FIG3 is the calculation domain of Example 1. The calculation domain size (length L = 1 mm, height H = 1 mm) contains solid particle phase, gas phase and liquid phase; the number of particles N = 1, the particle size D = 300 μm; the diameter of the droplet is 300 μm; the grid (N) under the Euler system is divided i ,N j ), where N i =128 and N j =128, which are the number of grids in the horizontal and vertical directions respectively; set the calculation time step Δt = 1×10 -5 s, total number of time steps N t =1×10 6 .

[0125] S02, initialization parameters: determine the velocity of all grid cell fluids at the initial moment The temperature of the fluid in all grid cells at the initial moment and the volume fraction F of the fluid in all grid cells at the initial moment 0 (x,y); determine the center position of all solid particles at the initial moment (x s ,y s ), the velocity of all grid unit solids at the initial moment The temperature of the solid in all grid cells at the initial moment and the volume fraction of solid in all grid cells at the initial moment Determine all thermophysical parameters, including: density of the first fluid ρ1, density of the second fluid ρ2, density of solid particles ρ s , viscosity of the first fluid μ1, viscosity of the second fluid μ2, thermal conductivity of the first fluid λ1, thermal conductivity of the second fluid λ2, thermal conductivity of the solid particles λ s .

[0126] In Example 1, the initialization parameters are: the first fluid is a droplet, the second fluid is a gas, and the center of the droplet is (N x =64, N y =75), at the center, the fluid volume fraction of the grid cells within a radius of 150 μm is 1, and that of other locations is 0; the velocity of the grid cells where the droplet and gas are located at the initial moment is temperature The center of the solid particle is (N x =64, N y =35), at the center, the solid volume fraction of the grid cells within a radius of 100 μm is 1, and the solid volume fraction of other locations is 0; the velocity of the solid grid cells Temperature of the grid cell Droplet density ρ1 = 1000 kg / m 3 , gas density ρ2=1kg / m3 , the density of solid particles ρ s =1000kg / m 3 , the viscosity of the droplet μ1=3×10 -3 Pa·s, viscosity of the second fluid μ2=2×10 -6 Pa·s, thermal conductivity of the droplet λ1 = 0.59 W / m·k, thermal conductivity of the second fluid λ2 = 0.02 W / m·k, thermal conductivity of the solid particles λ s =59W / m·k.

[0127] S03: Reconstruct the liquid-liquid interface, as shown in Figure 2:

[0128] When the fluid volume fraction F n (x,y)∈(0,1), then the network unit (N x ,N y ) is the liquid-liquid interface grid; when F n When (x,y)=0, the network unit (N x ,N y ) is filled with the second fluid; when F n When (x,y)=1, the network unit (N x ,N y ) is filled with a first fluid;

[0129] When the grid (x, y) is a liquid-liquid interface grid, the interface shape at that moment needs to be reconstructed. The piecewise linear interpolation technique PLIC (Particle Line-segment Interface Calculation) is used to reconstruct the interface shape of all grid cells at the nth moment. Specifically,

[0130] According to the grid unit (N x ,N y ) of the fluid volume fraction F n (x,y) determines the unit normal vector of the interface of the grid cell at the nth moment in, is the Hamiltonian operator, i represents the horizontal component, j represents the vertical component; where n is the number of time steps, n∈[1, N t ];

[0131] According to the grid unit (N x ,N y ) of the grid unit normal vector and the grid unit at the nth moment (N x ,N y ) of the fluid volume fraction F n(x,y) determines the grid cell (N x ,N y )Interface shape within the grid;

[0132] S04: Solve the fluid volume fraction at time n+1:

[0133] According to the grid unit (N x ,N y ) and the speed of the grid unit at the nth moment (N x ,N y ) fluid volume fraction F n (x,y) to obtain the inflow and outflow of the grid cell (N x ,N y ) of the fluid volume, and then obtain the fluid volume fraction F of the network unit at time n+1 n+1 (x,y), the specific solution equation is the VOF transport equation:

[0134] S05: When the grid (x, y) is a liquid-liquid interface grid, solve the signed distance function, as shown in Figure 4, specifically:

[0135] At the nth moment, the fluid volume fraction F of the grid unit n When (x,y) is 0.5, the grid cell is a boundary grid; determine all boundary grids of the computational domain;

[0136] Computational grid cells (N x ,N y ) to each boundary grid, thus obtaining the grid cell (N x ,N y ) to the minimum distance d among all boundary grids, that is, the grid unit (N x ,N y ) is the shortest distance d from the two-phase interface; as shown in Figure 4, the grid unit (N x ,N y ) to each boundary grid is represented by d 1b d 1b ,……d 8b ; and d 1b d 1b ,……d 8b Medium 4b Minimum, then the grid unit (N x ,N y ) to the minimum distance d from all boundary grids 4b .

[0137] According to the grid unit (N x ,N y) of the fluid volume fraction F n (x,y) determines the grid cell (N x ,N y )'s signed distance function

[0138] S06: Solve the interface curvature:

[0139] There are three situations:

[0140] The first type: when the liquid-liquid interface is formed, the grid unit (N x ,N y )'s signed distance function Determine the interface curvature k(x,y) using the following calculation formula:

[0141] The second type: when the grid unit (N x ,N y ) The fluid volume fraction F n (x,y)∈(0,1) and solid volume fraction When , it can be considered that the grid unit (N x ,N y ) is a three-phase contact unit, as shown in Figure 5. When the three phases are in contact, the contact angle boundary condition is used on the contact line to correct the interface unit normal vector of the mesh unit. The correction equation is:

[0142] Where θ is the contact angle, is the grid unit (N x ,N y ) solid wall normal vector, is the grid unit (N x ,N y ) solid wall tangent vector; in Example 1, the contact angles θ are 30°, 60°, 90°, 120° and 150° respectively.

[0143] Determine the grid cell (N x ,N y )’s interface curvature k(x,y):

[0144] The third type: when there is no interface, there is no curvature, that is, k(x,y)=0.

[0145] S07: Physical quantity of smooth interface mutation: Use Heaviside function to determine the grid unit (N at the nth moment) x ,N y) Density and viscosity near the smooth interface:

[0146] The grid unit at the nth moment (N x ,N y )The smooth density calculation formula is:

[0147] The grid unit at the nth moment (N x ,N y ) The formula for calculating smooth viscosity is:

[0148] Smooth Heaviside function for:

[0149] Where, ∈ is the thickness of the smoothed interface, which is generally taken as 1.5 or 2 times the grid size;

[0150] S08: Solving the surface tension of the interface

[0151] Solve the surface tension of the liquid-liquid interface described in S03. The grid unit (N x ,N y )Surface tension F s n (x,y) is calculated using the continuous surface force model. The specific calculation formula is as follows:

[0152] Where σ is the surface tension coefficient,

[0153] When there is no interface, there is no surface tension, and the grid unit (N x ,N y ) surface tension F s n (x,y)=0.

[0154] S09: Solving for the intermediate velocity of the fluid As shown in Figure 6:

[0155] According to the grid unit (N x ,N y ) to solve the Navier-Stokes equations by using the velocity, pressure and surface tension inside the grid and obtain the grid cell (N x ,N y ) of the intermediate speed The specific calculation formula is as follows:

[0156] Where Δt is the time step, is the convection term, is a sticky item;

[0157] Based on the fast Fourier transform method, the pressure Poisson equation is solved to obtain the grid unit (N at the n+1th moment) x ,N y ) pressure correction value P n+1 (x,y);

[0158] Using the grid unit (N at the n+1th moment x ,N y ) pressure correction value P n+1 (x,y) corrects the grid cell (N) at the n+1th moment x ,N y ) of the intermediate speed Get the corrected intermediate speed Specifically:

[0159] S10: Repeat steps S03-S9 to solve the fluid volume fraction F of all grid cells in the computational domain n (x,y) and intermediate speeds

[0160] S11: Solve for particle velocity:

[0161] Since particles are rigid bodies, the velocity of each particle at all positions is the same. The overall velocity of each particle can be obtained by solving the net force acting on each particle. The net force acting on each particle includes the fluid force and collision force Fluid forces on solid particles and collision force The IBM and DEM methods are used for calculation, as shown in Figure 7.

[0162] A. Fluid force calculation:

[0163] The fluid force on the pth particle at time n is determined by vector summing the fluid forces on the grid occupied by the particle. Based on the center position and radius of the pth particle at time n, the set of all grid cells occupied by the pth particle is determined. Among them, P e For collection The e-th grid cell in the set, E is the set The total number of grid cells within .

[0164] Solve the set Q according to the IBM method p Medium P e Fluid forces acting on the mesh This method can also make the solid-liquid interface meet the no-slip boundary condition. The specific calculation formula is as follows:

[0165] in:

[0166] P is the n+1th moment e The median speed of the grid; P is the nth moment e Solid velocity of the grid; P is the nth moment e The solid particle volume fraction of the grid;

[0167] According to the above formula, we can solve the set The fluid force on all grid cells at the nth moment is the fluid force on all grid cells occupied by the pth particle at the nth moment. The total fluid force on the pth particle is obtained by summing them up.

[0168] in, is the total fluid force on the pth particle.

[0169] B. Collision force calculation:

[0170] Determine the set of particles except the pth particle in Q q For the collection where q is the qth particle, and N is the number of particles in the computational domain.

[0171] Determine the collection of walls Where w represents different wall surfaces, w = 1, 2, 3, 4;

[0172] The DEM soft ball collision model is used to solve the collision force between the pth particle and the qth particle at the nth moment. and the collision force between the pth particle and the wth wall at the nth moment

[0173] When the distance between the center of the pth particle and the center of the qth particle is greater than D, the pth particle and the qth particle do not collide, and the collision force between the pth particle and the qth particle at the nth moment is That is 0;

[0174] When the distance between the center of the pth particle and the w wall is greater than D / 2, the pth particle does not collide with the w wall, and the collision force between the pth particle and the w wall at the nth moment is That is 0;

[0175] When the distance between the center of the pth particle and the center of the qth particle is less than D, the pth particle collides with the qth particle, and the collision force between the pth particle and the qth particle is The calculation formula is as follows:

[0176] Where: ε is the elastic coefficient, η is the viscous dissipation coefficient, and are the displacements in the tangential and normal directions, is the relative velocity between the pth and qth particles, and are the tangential and normal velocities of the pth particle, is the unit normal vector from the center of the pth particle to the center of the qth particle;

[0177] When the distance between the center of the pth particle and the w wall is less than D / 2, the pth particle collides with the w wall, and the collision force between the pth particle and the w wall is The calculation formula is as follows:

[0178] in: is the relative velocity between the pth particle and the wth wall, is the unit normal vector pointing from the center of the pth particle to the wth wall;

[0179] The total collision force on particle p at time n is the sum of the forces of collision with other particles or walls:

[0180] C. Particle velocity solution: Based on the total fluid force and total collision force on the pth particle, use Newton's momentum equation and angular momentum equation to solve the translational velocity of the pth particle. The rotation speed of the pth particle and the velocity vector of the pth particle The calculation formula is as follows:

[0181] Among them, m sp is the mass of the pth particle, I sp is the inertia tensor of the pth particle, is the radius vector from the center of the pth particle to the boundary of the pth particle.

[0182] D. Repeat the above steps AC to solve the velocity of all particles, and then obtain the velocity vectors of all particles in the calculation domain, the fluid force of all grid cells occupied by particles

[0183] S12. Solve for the solid particle displacement d n+1 and volume fraction

[0184] A. Solve particle displacement:

[0185] The velocity vector of the pth particle Decompose into horizontal and vertical velocity components and According to the pth particle speed and Solve the horizontal and vertical displacements of the pth particle in the next time step Δt respectively, and then obtain the position of the pth particle at the n+1th moment:

[0186] in, and are the horizontal and vertical coordinates of the pth particle at the nth moment, and is the component velocity of the pth particle in the horizontal and vertical directions at the nth moment, is the displacement of the pth particle from the nth moment to the n+1th moment.

[0187] B. Solve for the solid volume fraction

[0188] When the grid cell (N x ,N y ) from the center of the pth particle When the shortest distance is less than the particle radius D / 2, the grid unit (N x ,N y ) solid particle volume fraction

[0189] When the grid cell (N x ,N y ) from the center of the pth particle When the maximum distance is greater than the particle radius, the grid unit (N x ,N y ) solid particle volume fraction

[0190] When the grid cell (N x ,N y) from the center of the pth particle When the shortest distance is greater than D / 2 and the longest distance is less than D / 2, then the grid unit (N x ,N y ) is the solid-liquid boundary, and its volume fraction is calculated as follows: γ=0.05(1-β 2 )+0.3

[0191] Where h is the distance from the solid-liquid interface to the grid unit (N x ,N y ) is the distance between the centers; Δx is the grid size; Solid-liquid interface normal vector The component in the horizontal direction, Solid-liquid interface normal vector The component in the vertical direction,

[0192] C. Obtain the displacement of all particles and the solid volume fraction of all particles Then, the grid cell (N) at the n+1th moment in the computational domain is obtained according to the following formula: x ,N y ) solid volume fraction:

[0193] if but

[0194] S13: Calculate the final velocity u of the computational domain at the n+1th moment n+1 (x,y):

[0195] According to the relationship between action and reaction, the grid unit at the nth moment (N x ,N y The force exerted by the particles on the fluid in According to this force, the fluid velocity at time n+1 is obtained:

[0196] The final velocity u of the entire computational domain at the n+1th moment is solved using the following formula: n+1 (x,y):

[0197] S14. Solve the temperature and interface heat flux density. The specific steps are as follows:

[0198] A. Solve the particle interface heat flux density:

[0199] A-1. Based on the solid volume fraction at the n+1th moment Get all solid interface meshes (i.e. when ) Among them, H g For the collection The gth interface grid unit in the set G is The total number of interface mesh cells within .

[0200] A-2. On the interface grid H g The temperature gradient is divided into the normal direction and tangent direction Where i and j are the grid numbers in the horizontal and vertical directions respectively, T n (H g ) is the interface grid H at the nth moment g temperature, is the unit vector in the normal direction of the particle interface, and I is the unit matrix.

[0201] A-3. The thermal conductivity in the normal and tangential directions is calculated using the following formulas, as shown in Figure 8:

[0202] In the formula is the interface grid unit H e Thermal conductivity in the normal direction at the n+1th moment, is the tangential thermal conductivity at the n+1th moment. s is the thermal conductivity of the solid. f The thermal conductivity of the fluid, when the solid interface F n+1 (H g )>0, then λ f The value is λ f1 , otherwise the value is λ f2 .

[0203] A-4. Heat flux at the interface grid Calculated by the following formula:

[0204] B. Solution of fluid interface heat flux density:

[0205] B-1. According to the fluid volume fraction F at the n+1th moment n+1 (x,y) to obtain all liquid-liquid interface grids (i.e. when F n+1 (x,y)∈(0,1)) gets the union Among them, F h For the collection The hth interface grid unit in, H is the set The total number of interface mesh cells within .

[0206] B-2. On the interface grid F h The temperature gradient is divided into the normal direction and tangent direction Where i and j are the grid numbers in the horizontal and vertical directions respectively, and T is T n (F h ) is the interface grid F at the nth moment h The temperature, is the unit vector in the normal direction of the liquid-liquid interface, and I is the unit matrix.

[0207] B-3. ​​The thermal conductivity in the normal and tangential directions is calculated using the following formulas, as shown in Figure 8:

[0208] In the formula is the interface grid unit F h Thermal conductivity in the normal direction at the n+1th moment, is the tangential thermal conductivity at the n+1th moment. f1 and λ f2 are the thermal conductivities of the first fluid and the second fluid, respectively.

[0209] B-4. Heat flux at the interface mesh Calculated by the following formula:

[0210] C. Use the Crank-Nicolson method to update the diffusion term of the energy equation in time and solve the temperature field on the Euler system.

[0211] S15, determine whether the calculation is completed: if the current time step number n <N t , then repeat steps S03-S14; if the current time step number n≥N t , the calculation terminates.

[0212] Figure 9 is a simulation diagram of a droplet impacting a spherical solid particle in an embodiment of the present invention, and Figure 10 is a comparison result of the ratio of the minimum height of the droplet to the initial radius of the droplet after the droplet impacts the spherical solid particle in the embodiment of the present invention and the experiment, where the contact angle θ is 30°, 60°, 90°, 120° and 150°, respectively. It can be seen from the figure that the accuracy of the simulation results obtained using the calculation method of the present invention compared with the experiment can reach more than 98%.

[0213] FIG11 is a graph showing the volume change over time during the evaporation process of a liquid droplet on a solid particle in an embodiment of the present invention and in an experiment when the temperature of the solid particle is 120° C. As can be seen from the graph, the accuracy of the heat transfer results obtained using the calculation method of the present invention compared with the experiment can reach more than 97%.

[0214] A three-phase flow heat transfer coupling simulation system based on IBM-VOSET-DEM comprises a storage medium; the storage medium stores a program written using the three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM.

[0215] It should be understood that although this specification is described according to various embodiments, not every embodiment contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

[0216] The series of detailed descriptions listed above are only specific descriptions of feasible embodiments of the present invention. They are not intended to limit the scope of protection of the present invention. Any equivalent embodiments or changes that do not deviate from the technical spirit of the present invention should be included in the scope of protection of the present invention.

Claims

1. A three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM, characterized in that: The following steps are involved: S01: Determine the computational domain and the grids that divide the computational domain: Determine the size of the computational domain, the number of particles N and the particle size D in the computational domain; any grid in the computational domain is represented by (N x ,N y ), where x∈[1, N1]; y∈[1, N2]; set the calculation time step Δt, the total number of time steps N t ; Where N1 is the total number of grids in the horizontal direction; N2 is the total number of grids in the vertical direction; S02: Initialization parameters: Determine the velocity of all grid cell fluids at the initial moment The temperature of the fluid in all grid cells at the initial moment and the volume fraction F of the fluid in all grid cells at the initial moment 0 (x,y); determine the center position of all solid particles at the initial moment (x s ,y s ), the velocity of all grid unit solids at the initial moment The temperature of the solid in all grid cells at the initial moment and the volume fraction of solid in all grid cells at the initial moment Determine all thermophysical parameters, including: density of the first fluid ρ1, density of the second fluid ρ2, density of the solid particles ρ s , the viscosity of the first fluid μ1, the viscosity of the second fluid μ2, the thermal conductivity of the first fluid λ1, the thermal conductivity of the second fluid λ2, the thermal conductivity of the solid particles λ s ; S03: Reconstruction of liquid-liquid interface: When the fluid volume fraction F n (x,y)∈(0,1), then the network unit (N x ,N y ) is the liquid-liquid interface grid; when F n When (x,y) = 0, the network unit (N x ,N y ) is filled with the second fluid; when F n When (x,y)=1, the network unit (N x ,N y ) is filled with a first fluid; When the grid (x, y) is a liquid-liquid interface grid, the piecewise linear interpolation technique PLIC is used to reconstruct the interface shape of all grid cells at the nth moment; S04: Determine the fluid volume fraction at time n+1: According to the grid unit (N x ,N y ) and the speed of the grid unit at the nth moment (N x ,N y ) of the fluid volume fraction F n (x, y) to obtain the inflow and outflow of the grid cell (N) within Δt time x ,N y ) to obtain the fluid volume fraction F of the network unit at time n+1 n+1 (x,y); S05: Determine the boundary grid. When the grid (x, y) is a liquid-liquid interface grid, determine the grid unit (N) at the nth moment. x ,N y ) S06: Determine the grid unit of the interface (N x ,N y )’s interface curvature k(x,y): S07: Physical quantity of smooth interface mutation: Use Heaviside function to determine the grid unit at the nth moment (N x ,N y )Density and viscosity near the smooth interface: The grid unit (N x ,N y )The smooth density calculation formula is: The grid unit (N x ,N y ) The formula for calculating smooth viscosity is: Smooth Heaviside function for: Among them, ∈ is the thickness of the smooth interface; S08: Determine the surface tension of an interface Solve the surface tension of the liquid-liquid interface. The grid unit (N x ,N y )Surface tension F s n (x,y) is calculated by the continuous surface force model. The specific calculation formula is as follows: Where σ is the surface tension coefficient, When there is no interface, there is no surface tension, and the grid unit (N x ,N y ) surface tension F s n (x,y)=0; S09: Determine the intermediate velocity of the fluid S10: Repeat steps S03-S9 to obtain the fluid volume fraction F of all grid cells in the calculation domain n (x,y) and intermediate speeds S11: Determine particle velocity: IBM and DEM methods are used to determine the fluid forces on solid particles. and collision force Determine the velocity of the particle: Based on the total fluid force and total collision force on the pth particle, use Newton's momentum equation and angular momentum equation to solve the translational velocity of the pth particle. The rotation speed of the pth particle and the velocity vector of the pth particle The calculation formula is as follows: Among them, m sp is the mass of the pth particle, I sp is the inertia tensor of the pth particle, is the radius vector from the center of the pth particle to the boundary of the pth particle; is the total fluid force on the pth particle; is the total collision force on the pth particle at the nth moment Obtain the velocity of all particles, and then obtain the velocity vector of all particles in the computational domain, the fluid force of all grid cells occupied by the particles S12, determine the displacement d of the pth solid particle at the n+1th moment n+1 and the grid unit at the n+1th moment (N x ,N y ) solid volume fraction S13: Determine the final velocity u of the computational domain at the n+1th moment n+1 (x,y); S14: Solve the temperature and interface heat flux. The specific steps are as follows: A. Solve the particle interface heat flux density: A-1. According to the solid volume fraction at the n+1th moment Get all solid interface meshes and set all solid interface meshes to a collection Among them, H g For collection The gth interface grid unit in the set The total number of interface grid cells within; A-2. On the interface grid H g The temperature gradient is divided into the normal direction and tangent direction Where i and j are the grid numbers in the horizontal and vertical directions respectively, T n (H g ) is the interface grid H at the nth moment g The temperature, is the unit vector in the normal direction of the particle interface, and I is the unit matrix; A-3. The thermal conductivity in the normal and tangential directions is calculated using the following formula: In the formula is the interface grid unit H e Thermal conductivity in the normal direction at the n+1th moment, is the thermal conductivity in the tangential direction at the n+1th moment; s is the thermal conductivity of the solid; f The thermal conductivity of the fluid is n+1 (H g )>0, then λ f The value is λ f1 , otherwise the value is λ f2 ; A-4. Heat flux at the interface grid Calculated by the following formula: B. Solution of fluid interface heat flux density: B-1. According to the fluid volume fraction F at the n+1th moment n+1 (x,y) Determine all liquid-liquid interface meshes; Set all liquid-liquid interface meshes to a collection Among them, F h For collection The hth interface grid unit in , H is the collection The total number of interface grid cells within; B-2. On the interface grid F h The temperature gradient is divided into the normal direction and tangent direction Where i and j are the grid numbers in the horizontal and vertical directions respectively, and T is T n (F h ) is the interface grid F at the nth moment h The temperature, is the unit vector in the normal direction of the liquid-liquid interface, and I is the unit matrix; B-3. ​​The thermal conductivity in the normal and tangential directions is calculated using the following formula: In the formula is the interface grid unit F h Thermal conductivity in the normal direction at the n+1th moment, is the thermal conductivity in the tangential direction at the n+1th moment; f1 and λ f2 are the thermal conductivities of the first fluid and the second fluid, respectively; B-4. Heat flux at the interface grid Calculated by the following formula: C. Use the Crank-Nicolson method to update the diffusion term of the energy equation in time and solve the temperature field on the Euler system; S15: If the current time step number is n <N t , then repeat steps S03-S14 until the current time step number n=N t .

2. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 1 is characterized in that: The piecewise linear interpolation technique PLIC is used to reconstruct the interface shape of all grid cells at the nth moment, specifically: According to the grid unit (N x ,N y ) of the fluid volume fraction F n (x,y) determines the unit normal vector of the interface of the grid cell at the nth moment in, is the Hamiltonian operator, i represents the horizontal component, j represents the vertical component; where n is the number of time steps, n∈[1, N t ]; According to the grid unit (N x ,N y ) of the grid unit normal vector and the grid unit at the nth moment (N x ,N y ) of the fluid volume fraction F n (x,y) determines the grid unit (N x ,N y )The interface shape within the grid.

3. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 1 is characterized in that: Determine the boundary grid and solve the signed distance function Specifically: At the nth moment, the fluid volume fraction F of the grid unit n When (x,y) is 0.5, the grid cell is a boundary grid; determine all boundary grids of the computational domain; Computational grid unit (N x ,N y ) to each boundary grid, thus obtaining the grid cell (N x ,N y ) to the minimum distance d among all boundary grids, that is, the grid unit (N x ,N y ) the shortest distance d from the interface between the two phases; According to the grid unit (N x ,N y ) of the fluid volume fraction F n (x,y) determines the grid unit (N x ,N y ) 4. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 3 is characterized in that: Determine the mesh unit (N) of the interface x ,N y )’s interface curvature k(x,y), specifically: When the liquid-liquid interface appears, according to the grid unit (N x ,N y ) Determine the interface curvature k(x,y), the calculation formula is as follows: When the grid unit (N x ,N y ) The volume fraction of the fluid in n (x,y)∈(0,1) and solid volume fraction When, it is the grid unit (N x ,N y ) is a three-phase contact unit; when the three phases are in contact, the contact angle boundary condition is used on the contact line to correct the interface unit normal vector of the grid unit, and the correction equation is: Where θ is the contact angle, is the grid unit (N x ,N y ) solid wall normal vector, is the grid unit (N x ,N y ) solid wall tangent vector; Determine the grid unit (N x ,N y )’s interface curvature k(x,y):

5. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 1 is characterized in that the intermediate velocity of the fluid is determined Specifically: According to the grid unit (N x ,N y ) to solve the Navier-Stokes equations for velocity, pressure and surface tension, and obtain the grid cell (N x ,N y ) of the intermediate speed The specific calculation formula is as follows: Where Δt is the time step, is the convection term, is a sticky item; Based on the fast Fourier transform method, the pressure Poisson equation is solved to obtain the grid unit (N) at the n+1th moment. x ,N y ) pressure correction value P n+1 (x,y); Using the grid unit (N at the n+1th moment x ,N y ) pressure correction value P n+1 (x,y) corrects the grid cell (N) at the n+1th moment x ,N y ) of the intermediate speed Get the corrected intermediate speed Specifically:

6. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 1 is characterized in that: IBM and DEM methods are used to determine the fluid forces on solid particles. and collision force Specifically: Determination of fluid forces on solid particles using the IBM method The fluid force on the pth particle at the nth moment is determined by vector summing the fluid forces on the grid occupied by the particle; according to the center position and radius of the pth particle at the nth moment, the set of all grid cells occupied by the pth particle is determined Among them, P e For collection The e-th grid cell in the set, E is the set The total number of grid cells within; Solve the set Q according to the IBM method p Medium e Fluid forces acting on the mesh This method can also make the solid-liquid interface meet the no-slip boundary condition. The specific calculation formula is as follows: in: P is the n+1th moment e The middle speed of the grid; P is the nth moment e Solid velocity of the grid; P is the nth moment e The solid particle volume fraction of the grid; Get a collection The fluid force on all grid cells at the nth moment is the fluid force on all grid cells occupied by the pth particle. The fluid force on each grid cell at the nth moment is summed up to obtain the total fluid force on the pth particle. in, is the total fluid force on the pth particle; Determining the impact force on solid particles using DEM method Determine the collection of particles other than the pth particle in Q q For collection The qth particle in the calculation domain, N is the number of particles in the calculation domain; Determine the collection of walls Where w represents different wall surfaces, w = 1, 2, 3, 4; The DEM soft ball collision model is used to solve the collision force between the pth particle and the qth particle at the nth moment. and the collision force between the pth particle and the wth wall at the nth moment When the distance between the center of the pth particle and the center of the qth particle is greater than D, the pth particle and the qth particle do not collide, and the collision force between the pth particle and the qth particle at the nth moment is That is 0; When the distance between the center of the pth particle and the w wall is greater than D / 2, the pth particle does not collide with the w wall. The collision force between the pth particle and the w wall at the nth moment is That is 0; When the distance between the center of the pth particle and the center of the qth particle is less than D, the pth particle collides with the qth particle, and the collision force between the pth particle and the qth particle is The calculation formula is as follows: Where: ε is the elastic coefficient, η is the viscous dissipation coefficient, and are the displacements in the tangential and normal directions, is the relative velocity between the pth and qth particles, and are the tangential and normal velocities of the pth particle, is the unit normal vector pointing from the center of the pth particle to the center of the qth particle; When the distance between the center of the pth particle and the w wall is less than D / 2, the pth particle collides with the w wall, and the collision force between the pth particle and the w wall is The calculation formula is as follows: in: is the relative velocity between the pth particle and the wth wall, is the unit normal vector pointing vertically from the center of the pth particle to the wth wall; The total collision force on the pth particle at time n It is the sum of the collision forces with other particles or walls:

7. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 6, characterized in that: Determine the displacement d of the pth solid particle at the n+1th moment n+1 , specifically: The velocity vector of the pth particle Decompose into horizontal and vertical velocity components and According to the pth particle speed and Solve the horizontal and vertical displacements of the pth particle in the next time step Δt respectively, and then obtain the position of the pth particle at the n+1th moment: in, and are the horizontal and vertical coordinates of the pth particle at the nth moment, and is the component velocity of the pth particle in the horizontal and vertical directions at the nth moment, is the displacement of the pth particle from the nth moment to the n+1th moment.

8. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 6, characterized in that: Determine the grid cell (N) at the n+1th moment x ,N y ) solid volume fraction Specifically: When the grid unit (N x ,N y ) from the center of the pth particle When the shortest distance is less than the particle radius D / 2, the grid unit (N x ,N y ) solid particle volume fraction When the grid unit (N x ,N y ) from the center of the pth particle When the maximum distance is greater than the particle radius, the grid unit (N x ,N y ) solid particle volume fraction When the grid unit (N x ,N y ) from the center of the pth particle The shortest distance is greater than D / 2 and the longest distance is less than D / 2 When the grid unit (N x ,N y ) is the solid-liquid boundary, and its volume fraction is calculated as follows: γ=0.05(1-β 2 )+0.3 Where h is the distance from the solid-liquid interface to the grid unit (N x ,N y ) is the distance from the center; Δx is the grid size; Solid-liquid interface normal vector The component in the horizontal direction, Solid-liquid interface normal vector The component in the vertical direction, According to the displacement of all particles and the solid volume fraction of all particles Determine the grid cell (N) at the n+1th moment in the computational domain x ,N y ) solid volume fraction if but 9. The three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM according to claim 6, characterized in that: Determine the final velocity u of the computational domain at the n+1th moment n+1 (x,y), specifically: The grid cell at the nth moment (N x ,N y The force exerted on the fluid by the particle is according to Get the fluid velocity at time n+1: The final velocity u of the entire computational domain at the n+1th moment is determined using the following formula: n+1 (x,y):

10. A three-phase flow heat transfer coupling simulation system based on IBM-VOSET-DEM, characterized in that: It comprises a storage medium; the storage medium stores a program written by the three-phase flow heat transfer coupling simulation method based on IBM-VOSET-DEM as described in any one of claims 1 to 9.

Images