A method for simulating the freezing process of a droplet impacting a cold surface

Abstract

The present invention relates to the technical field of computational fluid dynamics, and specifically to a method for simulating the freezing process of a droplet impacting a cold surface. The specific technical solution is: using the obtained distribution function to perform migration calculations; calculating the macroscopic physical quantities in the flow field based on the mass distribution function and the energy distribution function, and calculating the pressure and pseudo-potential energy of the fluid; calculating the pseudo-potential energy of the virtual fluid, and calculating the interaction forces between fluid particles and the interaction forces between fluid particles and solid walls; calculating the real macroscopic velocity, enthalpy value and solid phase volume fraction of the fluid, the collision term of the distribution function, the force source phase, and the heat source term; performing collision calculations on the mass distribution function and the energy distribution function, and judging whether the convergence condition or the preset number of iterations is met, and outputting the calculation results if it is met; otherwise, continuing the iterative process. The present invention can provide a means for exploring the freezing characteristics and mechanisms of tiny droplets impacting low-temperature surfaces, and provide ideas for the design and development of anti-icing surfaces.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular to a method for simulating a freezing process of a droplet impacting a cold surface based on a lattice Boltzmann method. Background Art

[0002] The freezing of droplets impacting low-temperature surfaces is common in many natural and industrial scenarios, such as aircraft, ships, and wind turbines operating in extremely cold environments. The icing phenomenon caused by droplet impact and freezing seriously threatens the safety and heat exchange performance of related equipment, leading to significant economic losses and safety accidents. Superhydrophobic anti-icing surface technology has significant advantages in reducing the contact time between impacting droplets and cold surfaces, delaying droplet freezing, and weakening the adhesion of ice layers, and has now been widely used. However, it still has problems such as poor stability, weak mechanical strength, and difficulty in coping with complex freezing conditions. Therefore, there is an urgent need to develop more efficient and reliable anti-icing surfaces. Numerical simulation of the icing process is one of the important technical means for the design and development of anti-icing surfaces.

[0003] The freezing process of droplet impacting cryogenic surface involves the coupling of complex multiphase flow motion and heat transfer with solidification phase change, so the simulation of this process is always a challenging task. The main means of simulating the freezing process of droplets at present can be divided into macroscopic methods (finite element method, finite volume method, etc.), microscopic methods (molecular simulation), and mesoscopic methods (lattice Boltzmann method). Macroscopic methods are currently more mature simulation calculation means, but they are difficult to capture the microscopic motion data of droplets and can ignore the important microscopic features of the freezing process of droplets; the microscopic methods headed by molecular simulation can well illustrate the microscopic mechanism of freezing process, but are limited by existing computer hardware technology, and their calculation results are difficult to be extended to large-scale freezing process. Therefore, the present invention uses mesoscopic particle beam as basic simulation object by lattice Boltzmann method, and its calculation results can both consider the interaction between particles and be extended to general macroscopic scale. Summary of the Invention

[0004] In response to the shortcomings of the existing technology, the present invention provides a method for simulating the freezing process of a droplet impacting a cold surface. The method can simulate the movement processes of a droplet impacting a low-temperature surface, such as diffusion, contraction, and rebound, as well as the freezing phase transition process coupled therewith.

[0005] To achieve the above objectives, the present invention is implemented through the following technical solutions:

[0006] The present invention discloses a method for simulating the freezing process of a droplet impacting a cold surface, which can simulate the coupled process of the motion and freezing of a droplet impacting a low-temperature surface, comprising the following steps:

[0007] (1) Obtaining computational area grid information and fluid medium information;

[0008] (2) Determine the macroscopic physical quantity distribution of the initial flow field and calculate the equilibrium distribution function f according to the macroscopic physical quantity of the flow field i eq and As the mass distribution function f in the initial state i * (x, 0) and the energy distribution function

[0009] (3) Using the obtained distribution function to perform migration calculations, calculate the mass distribution function and energy distribution function after migration;

[0010] (4) calculating macroscopic physical quantities in the flow field based on the mass distribution function and the energy distribution function, wherein the macroscopic physical quantities are density ρ, fluid velocity u during migration, and fluid temperature T;

[0011] (5) Based on the calculated density and temperature distribution, the pressure p and pseudo potential energy ψ(x) of the fluid are calculated using the PR state equation; the solid wall boundary of the flow field in the calculation area is regarded as a virtual fluid, which has pseudo potential energy and generates an interaction force with the surrounding flow field. Its pseudo potential energy is the weighted average of the pseudo potential energy of the surrounding real fluid; based on the potential energy distribution of the fluid, the interaction force between the fluid particles and the interaction force between the fluid particles and the solid wall are calculated, and the force F acting on the fluid is calculated;

[0012] (6) According to the formula U=u+0.5F / δ t Calculate the real macroscopic velocity U of the fluid, the enthalpy H and the solid volume fraction f s , calculate the collision terms of the mass distribution function and energy distribution function and Force source phase Δf i (x, t), heat source term φ;

[0013] (7) Perform collision calculation on the mass distribution function and energy distribution function, and calculate the mass distribution function f after collision i * and energy distribution function At the same time, the mass distribution function and energy distribution function at the boundary of the flow field area after the collision are calculated according to the boundary conditions.

[0014] (8) Determine whether the convergence condition is met or whether the preset number of iterations is reached. If so, output the calculation result and end the calculation program; otherwise, use the f obtained in step (7) i * (x, t) and Return to step (3) to continue calculation.

[0015] Preferably, in step (1), the calculation area includes the flow field area and the solid wall area, and the fluid medium information includes the viscosity ν, thermal conductivity k, and specific heat capacity c of the simulated fluid medium. p , phase transition temperature T m and the latent heat of phase change ΔH ls .

[0016] Preferably, in step (2), the macroscopic physical quantities include the density ρ, velocity U, temperature T, solid fraction f of the flow field in the initial state s distributed.

[0017] Preferably, in step (2), the equilibrium distribution function f i eq and The calculation formula is:

[0018]

[0019] Preferably, in step (3), the migration calculation formula is:

[0020] f i (x+e i δ t , t+δ t )=f i * (x, t)

[0021]

[0022] Among them, f i (x+e i δ t , t+δ t ) and g i (x+e i δ t , t+δ t ) is the distribution function corresponding to the current moment, f i * (x, t) and is the distribution function corresponding to the previous moment. When step (3) is executed for the first time, f i * (x, t) and Corresponding to f in step (2) i * (x, 0) and

[0023] Preferably, in step (4), the calculation formulas for density ρ, fluid velocity u during migration, and fluid temperature T are:

[0024] ρ=∑ i fi

[0025] T=∑ i g i

[0026] ρu=∑ i e i f i .

[0027] Preferably, in step (5), the calculation formulas for the pressure p and pseudo potential energy ψ(x) of the fluid are:

[0028]

[0029] in, b=0.0778RT cr / p cr , the expression of ε(T0) is:

[0030] It should be noted that T0 in the state equation is a set value, and its value depends on the selected liquid-gas density ratio;

[0031] The formula for calculating the pseudo potential energy of the virtual fluid is:

[0032] ψ(x s )=∑w i ψ(x′)(1-s(x′)) / ∑w i (1-s(x′))

[0033] Where s(x′) is an indicator function. When the region where the x′ grid point is located is a solid wall, its value is 1, otherwise it is 0.

[0034] The interaction force between fluid particles is calculated as:

[0035]

[0036] Where G(x, x′) represents the interaction force between fluid particles, and x′ represents the surrounding fluid grid. Due to the different wettability between the stationary droplet and the solid surface, there will be an additional force between the fluid and the solid surface, as shown in the following formula:

[0037] F s (x)=-ψ 2 (x)∑ i G s w i s(x′)(x′-x)

[0038] Among them, G s The value of G depends on the wettability of the surface. s >0 means the surface is hydrophilic, G s<0 indicates a hydrophobic surface.

[0039] Preferably, in step (6), the enthalpy value H and the solid phase volume fraction f s The calculation formulas are:

[0040]

[0041] During the calculation process, the enthalpy value H is calculated first, and then H is used to update the solid phase volume fraction f s ;

[0042] Mass distribution function and energy distribution function collision terms and The calculation formulas are:

[0043]

[0044] Force source phase Δf i The formula for (x, t) is:

[0045] Δf i (x, t) = f i eq (ρ,u+Δu)-f i eq (ρ, u)

[0046] Where Δu=Fδ t / ρ, F is the cooperative force, including the force between fluid particles, the force between the fluid and the wall, and gravity: F = F int (x)+F g (x)+F s (x);

[0047] The heat source term is calculated as: The heat source term is divided into the sum of two parts. The former is the correction term of the heat transfer equation to the lattice Boltzmann equation, and the latter represents the temperature change caused by the solidification phase transition.

[0048] Preferably, in step (7), the collision formulas of the mass distribution function and the energy distribution function are:

[0049]

[0050] Among them, τ f , τ g For relaxation time, is the no-slip equivalent rebound collision term, which represents the particle collision mode when the fluid medium is in a solidified state. The expression is:

[0051]

[0052] Among them, u sis the velocity of the solidified state, Represents the discrete velocity vector number in opposite direction;

[0053] B is the weight coefficient, and its value depends on the phase state of the fluid medium. Its calculation formula is:

[0054]

[0055] The present invention has the following beneficial effects:

[0056] 1. The present invention is based on a particle beam that is much larger than the microscopic scale and much smaller than the macroscopic control body. By solving a double-distributed lattice Boltzmann equation group with a pseudo-potential energy model, a collision migration process, and a phase change process heat source term, the multiphase flow field and temperature field changes in the process of droplets impacting the cold surface are calculated. In order to make the simulation of the gas-liquid interface more realistic and reliable, the state equation is introduced into the calculation for the calculation of the fluid pseudo-potential energy. At the same time, the immersed moving boundary format is used to calculate the collision and migration of particles at the gas-solid and liquid-solid boundaries, and then the movement of the phase interface is calculated. The present invention uses mesoscopic particle clusters as the basic object to solve the complex physical field changes in the process of droplets impacting and freezing on supercooled surfaces, captures the droplet movement and phase interface change data that are difficult to obtain in the experimental process, and provides ideas for the design and development of anti-icing surfaces.

[0057] 2. The coupled motion and freezing of impacting droplets on cryogenic surfaces is a complex physical process. This invention provides a mesoscopic solution for simulating this complex multiphase flow physics. By using mesoscopic particle clusters as the primary research object, it is possible to capture data on the microscopic motion and solidification phase changes of droplets, which are difficult to obtain in macroscopic simulations and experiments. This invention provides a means to explore the freezing characteristics and mechanisms of tiny droplets impacting cryogenic surfaces, offering insights into the design and development of anti-icing surfaces. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 is the relationship between the algorithm process and each mathematical model, where for

[0059] Figure 2 Schematic diagram of droplet impact in Example 1;

[0060] Figure 3 The droplet of Example 1 hits the surface of the low-temperature array structure (τ * is dimensionless time ) simulation results (T w =-18°C);

[0061] Figure 4 The droplet of Example 1 hits the surface of the low-temperature array structure (τ * is dimensionless time ) simulation results (T w =-26°C);

[0062] Figure 5 The droplet of Example 1 hits the surface of the low-temperature array structure (τ * is dimensionless time ) simulation results (T w =-34°C). DETAILED DESCRIPTION

[0063] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0064] Unless otherwise specified, the technical means used in the embodiments are conventional means well known to those skilled in the art.

[0065] Example 1

[0066] The present invention discloses a method for simulating the freezing process of a droplet hitting a cold surface. Figure 1 As shown, we now take the simulation of the freezing process of a droplet impacting a supercooled surface with an array structure as an example. The droplet impact model is as follows Figure 2 As shown, the present invention is used as a simulation method to explore the freezing characteristics and mechanism of tiny droplets impacting a cryogenic surface. The model construction process includes the following steps:

[0067] (1) Constructing the calculation model: the calculation area is as follows: Figure 2 As shown in the figure, the solid wall area below the calculation area is the target surface for simulating droplet impact. It is a super-hydrophobic surface with an array of conical structures. The solid wall boundary condition is calculated using the no-slip rebound format. The open boundary condition above the calculation area adopts the non-equilibrium extrapolation boundary format. The periodic boundary conditions are used around the calculation area. The calculation area is divided into regular hexahedral grids, and the grid cells can be divided into fluid nodes and solid nodes. The calculation model adopts the D3Q19 scheme, and its discrete velocity and weighting coefficient are:

[0068]

[0069] (2) Initialization of the calculation model: At the initial moment, the droplet is located above the target surface. In the flow field area where the droplet is located, its initial density is set to the liquid phase density ρ l , the remaining flow field areas are set to gas phase density ρ g The initial temperatures of the fluid and solid wall are set to Tstr , the freezing phase transition temperature of the droplet is set to T m (i.e., freezing point). The solid fraction f at the initial moment s is 0, the flow field is in a static state and the velocity is 0. Through the equilibrium distribution function f i eq and As the mass distribution function f in the initial state i * (x, 0) and the energy distribution function The calculation formula of the equilibrium distribution function is:

[0070]

[0071] In the D3Q19 scheme, u is the fluid velocity, c s is the lattice sound speed, and its value is

[0072] (3) In the pseudo-potential energy model, the density gradient and flow at the gas-liquid interface cannot be directly set artificially. Therefore, after initialization, the calculation model needs to be set in a zero-gravity constant temperature (T * =T str ) under the condition of evolution iteration for 1000 steps, so that the gas-liquid interface of the droplet in the flow field tends to be stable. At the 1001st moment, a downward force is applied to the droplet to make it obtain an initial vertical downward velocity, and the temperature of the target surface is set to T w , and remains unchanged in the subsequent calculation process (T w <T m <T str ).

[0073] The specific steps of model iterative evolution are as follows:

[0074] (1) Using the obtained mass distribution function and energy distribution function to perform migration calculation, the mass distribution function and energy distribution function after migration are calculated. The formula is:

[0075] f i (x+e i δ t , t+δ t )=f i * (x, t)

[0076]

[0077] Among them, f i (x+e i δ t , t+δ t ) and g i (x+ei δ t , t+δ t ) is the distribution function corresponding to the current moment, f i * (x, t) and is the distribution function corresponding to the previous moment. When the migration calculation is performed for the first time, f i * (x, t) and Corresponding to the mass distribution function and energy distribution function at the initial moment, δ t is the unit time step.

[0078] (2) Calculate the macroscopic physical quantities in the flow field based on the distribution function. The macroscopic physical quantities that need to be calculated are density ρ, fluid velocity u during migration, and fluid temperature T. The calculation formulas are:

[0079] ρ=∑ i f i

[0080] T=∑ i g i

[0081] ρu=∑ i e i f i .

[0082] (3)(3.1) Based on the calculated density and temperature distribution, the PR state equation is used to calculate the fluid pressure p and pseudo potential energy ψ(x).

[0083] (a) The PR state equation is as follows:

[0084]

[0085] The parameters in the state equation are b=0.0778RT cr / p cr , T cr and p cr represent the critical temperature and critical pressure respectively, R is a thermodynamic constant, and the expression of ε(T0) is:

[0086]

[0087] In the embodiment, a=2 / 49, b=2 / 21, R=1, and ω=0.344 are set. It should be noted that T0 in the state equation is a set value, and its value depends on the selected liquid-gas density ratio. In this embodiment, the liquid-gas density ratio of the multiphase flow field is 15:1, so T0 is 0.87T. cr .

[0088] (b) The calculation formula of pseudo potential energy is:

[0089]

[0090] c0 is a model parameter, and its value is only related to the discrete scheme of particle velocity. In the D3Q19 scheme, c0 is 6; G is a formal parameter, and its value has no effect on the calculation.

[0091] (3.2) In order to calculate the interaction force between the fluid and the wall, it is necessary to regard the solid wall boundary of the flow field in the calculation area as a virtual fluid, which has pseudo potential energy and generates an interaction force with the surrounding flow field. Its pseudo potential energy is the weighted average of the pseudo potential energy of the surrounding real fluid, that is:

[0092] ψ(x s )=∑w i ψ(x′)(1-s(x′)) / ∑w i (1-s(x′))

[0093] Where s(x′) is an indicator function. When the area where the x′ grid point is located is a solid wall, its value is 1, otherwise it is 0. Where x′=x+e i δ t , that is, the coordinates of the centers of the grid points around the grid point centered on x.

[0094] (3.2) Based on the potential energy distribution of the fluid, the interaction forces between the fluid particles and the interaction forces between the fluid particles and the solid wall are calculated, and the force F acting on the fluid is calculated.

[0095] The calculation formula for the interaction force between fluid particles is:

[0096]

[0097] Among them, G(x, x′) is the Green's function, which represents the interaction force between fluid particles and its expression is:

[0098]

[0099] For the D3Q19 scheme, G1=G, G2=G / 2.

[0100] In addition, due to the different wettability of stationary droplets and solid surfaces, in order to achieve different droplet contact angles, the fluid will require additional forces with the solid surface:

[0101] F s (x)=-ψ 2 (x)∑ i G s w i s(x′)(x′-x)

[0102] Among them, G s The value of G depends on the wettability of the surface. s >0 means the surface is hydrophilic, G s <0 indicates a hydrophobic surface.

[0103] In this embodiment, G s The value is -1.25, and the corresponding surface contact angle is 160°.

[0104] (4) According to the formula U=u+0.5F / δ t Calculate the true macroscopic velocity U of the fluid. During the simulation process, the fluid velocity u during the migration process is the intermediate variable of the calculation. Only the true macroscopic velocity can correctly reflect the true velocity of the flow field.

[0105] (4.1) Update the enthalpy value H and solid volume fraction f of the flow field node s , which is calculated as follows:

[0106]

[0107]

[0108] In each iteration, first pass f s Calculate the enthalpy value H of the node, and then use the enthalpy value to update the solid phase volume fraction of the node. is the liquid fraction, is the gas phase fraction, ΔH ls is the latent heat of solid-liquid phase change, H s ,H l ,H g are the enthalpy values ​​of the solid, liquid, and gas phases, respectively, and ΔH ls =H l -H s .

[0109] (4.2) Mass distribution function and energy distribution function collision terms and The calculation formula is:

[0110]

[0111]

[0112] Note that when calculating the collision term of the mass distribution function f, the equilibrium distribution function requires the use of the fluid velocity u during the migration process, while when calculating the collision term of the energy distribution function g, the true velocity U is required.

[0113] (4.3) In this embodiment, the force source phase Δf i(x, t) uses the EDM format, and the formula is:

[0114] Δf i (x, t) = f i eq (ρ,u+Δu)-f i eq (ρ, u)

[0115] Where Δu=Fδ t / ρ, F is the cooperative force, including the force between fluid particles, the force between the fluid and the wall, and gravity: F = F int (x)+F g (x)+F s (x).

[0116] (4.4)The calculation of the heat source term φ can be divided into two parts, as shown below:

[0117]

[0118] The former is the correction term of the heat transfer equation to the lattice Boltzmann equation, and the latter represents the temperature change caused by the solidification phase change. s0 is the solid volume fraction of the node at the previous moment.

[0119] (5) Perform collision calculation on the mass distribution function and energy distribution function, and calculate the mass distribution function f after the collision i * and energy distribution function The formula is:

[0120]

[0121] Among them, τ f , τ g For relaxation time, is the no-slip equivalent rebound collision term, which represents the particle collision mode when the fluid medium is in a solidified state. The expression is:

[0122]

[0123] Among them, u s is the movement speed in the solidified state, and its value is 0 in this embodiment. Represents the discrete velocity vector number in the opposite direction.

[0124] B is the weight coefficient, and its value depends on the phase state of the fluid medium. Its calculation formula is:

[0125]

[0126] When f s= 0, B = 0, and the collision mode of the mass distribution function f is the same as that of the conventional lattice Boltzmann method, that is:

[0127]

[0128] (5) Finally, determine whether the convergence condition is met or whether the preset number of iterations is reached. If so, output the calculation result and end the calculation program. Otherwise, return to (1) to continue the iterative process.

[0129] The simulation results of this embodiment are as follows Figure 3-5 As shown in the figure, at lower surface supercooling, the droplet can rebound after contacting the cold surface, while increasing the supercooling makes it more difficult for the droplet to rebound, and the droplet freezes and adheres to the surface.

[0130] The embodiments described above are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by persons skilled in the art should fall within the scope of protection defined by the claims of the present invention.

Claims

1. A method for simulating the freezing process of a droplet impacting a cold surface, characterized in that: The following steps are involved: (1) Obtaining computational area grid information and fluid medium information; (2) Determine the macroscopic physical quantity distribution of the initial flow field and calculate the equilibrium distribution function based on the macroscopic physical quantity of the flow field and As the mass distribution function f in the initial state i * (x, 0) and the energy distribution function (3) Using the obtained distribution function to perform migration calculations, calculate the mass distribution function and energy distribution function after migration; (4) calculating macroscopic physical quantities in the flow field based on the mass distribution function and the energy distribution function, wherein the macroscopic physical quantities are density ρ, fluid velocity u during migration, and fluid temperature T; (5) Based on the calculated density and temperature distribution, the pressure p and pseudo potential energy ψ(x) of the fluid are calculated using the PR state equation; the solid wall boundary of the flow field in the calculation area is regarded as a virtual fluid, which has pseudo potential energy and generates an interaction force with the surrounding flow field. Its pseudo potential energy is the weighted average of the pseudo potential energy of the surrounding real fluid; based on the potential energy distribution of the fluid, the interaction force between the fluid particles and the interaction force between the fluid particles and the solid wall are calculated, and the force F acting on the fluid is calculated; (6) According to the formula U=u+0.5F / δ t Calculate the real macroscopic velocity U of the fluid, the enthalpy H and the solid volume fraction f s , calculate the collision terms of the mass distribution function and energy distribution function and Force source phase Δf i (x, t), heat source term φ; (7) Perform collision calculation on the mass distribution function and energy distribution function, and calculate the mass distribution function f after collision i * and energy distribution function At the same time, the mass distribution function and energy distribution function at the boundary of the flow field area after the collision are calculated according to the boundary conditions; (8) Determine whether the convergence condition is met or whether the preset number of iterations is reached. If so, output the calculation result and end the calculation program; otherwise, use the f obtained in step (7) i * (x, t) and Return to step (3) to continue calculation.

2. A method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (1), the calculation area includes the flow field area and the solid wall area, and the fluid medium information includes the viscosity ν, thermal conductivity k, and specific heat capacity c of the simulated fluid medium. p , phase transition temperature T m and the latent heat of phase change ΔH ls .

3. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (2), the macroscopic physical quantities include the density ρ, velocity U, temperature T, solid fraction f of the flow field in the initial state s distributed.

4. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (2), the equilibrium distribution function and The calculation formula is:

5. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (3), the migration calculation formula is: Among them, f i (x+e i δ t , t+δ t ) and g i (x+e i δ t , t+δ t ) is the distribution function corresponding to the current moment, f i * (x, t) and is the distribution function corresponding to the previous moment. When step (3) is executed for the first time, f i * (x, t) and Corresponding to f in step (2) i * (x, 0) and 6. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (4), the calculation formulas for density ρ, fluid velocity u during migration, and fluid temperature T are: p=∑ i f i T=∑ i g i ρu=∑ i and i f i 。 7. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (5), the calculation formulas for the fluid pressure p and pseudo potential energy ψ(x) are: in, b=0.0778RT cr / p cr , the expression of ε(T0) is: The formula for calculating the pseudo potential energy of the virtual fluid is: ψ(x s )=∑w i ψ(x ′ )(1-s(x ′ )) / ∑w i (1-s(x ′ )) Among them, s(x ′ ) is the indicator function, when x ′ When the grid point is located in a solid wall, its value is 1, otherwise it is 0; The interaction force between fluid particles is calculated as: Among them, G(x, x ′ ) represents the interaction force between fluid particles, x ′ represents the surrounding fluid grid; due to the different wettability between the stationary droplet and the solid surface, there will be an additional force between the fluid and the solid surface, as shown in the following formula: F s (x)=-ψ 2 (x)∑ i G s w i s(x ′ )(x ′ -x) Among them, G s The value of G depends on the wettability of the surface. s >0 means the surface is hydrophilic, G s <0 indicates a hydrophobic surface.

8. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (6), the enthalpy value H and the solid phase volume fraction f s The calculation formulas are: During the calculation process, the enthalpy value H is calculated first, and then H is used to update the solid phase volume fraction f s ; Mass distribution function and energy distribution function collision terms and The calculation formulas are: Force source phase Δf i The formula for (x, t) is: Where Δu=Fδ t / ρ, F is the cooperative force, including the force between fluid particles, the force between the fluid and the wall, and gravity: F = F int (x)+F g (x)+F s (x); The heat source term is calculated as:

9. The method for simulating the freezing process of a droplet impacting a cold surface according to claim 1, characterized in that: In step (7), the collision formulas of mass distribution function and energy distribution function are: Among them, τ f , τ g For relaxation time, is the no-slip equivalent rebound collision term, which represents the particle collision mode when the fluid medium is in a solidified state. The expression is: Among them, u s is the velocity of the solidified state, Represents the discrete velocity vector number in opposite direction; B is the weight coefficient, and its value depends on the phase state of the fluid medium. Its calculation formula is:

Images