Method for calculating impact characteristics of surface-icing water droplets based on lattice Boltzmann method

Abstract

This paper proposes a method for calculating the impact characteristics of frozen water droplets on surfaces based on the lattice Boltzmann method. This method, based on mesoscopic fluid clusters, which are much smaller than the macroscopic control volume, solves the multicomponent lattice Boltzmann equation with large eddy simulation, a splitting collision model, and a differential lattice velocity format. The method uses equilibrium state degradation of air components and the introduction of the polymerization number of water vapor molecules to make the model more physically realistic, thereby obtaining macroscopic physical quantities of air and water droplets, and further obtaining the impact characteristics of water droplets. This method can solve the air-water mixed flow field using fluid molecular clusters as the research object, capturing the microscopic motion of water droplets that cannot be simulated using macroscopic methods. Under complex freezing conditions, this method is beneficial for resolving grid scale constraints faced by complex geometry when solving flow fields.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and in particular to a method for calculating the impact characteristics of ice droplets on a surface based on a lattice Boltzmann method. Background Art

[0002] When an aircraft passes through a cloud layer containing supercooled water droplets, icing will occur. In the process of numerical simulation of icing, the calculation of water droplet impact characteristics is the basis for predicting icing patterns. The impact area of water droplets on the surface, the impact amount, and the distribution of water droplets in the impact area are collectively referred to as the impact characteristics of water droplets on the surface, such as Figure 1 The water droplet impact limit and the local water droplet collection coefficient are important parameters of the impact characteristics.

[0003] The lattice Boltzmann method is a mesoscopic computational fluid dynamics method that builds a bridge between microscopic molecular motion and macroscopic fluid flow. It solves the particle distribution function (PDF) based on the collision and migration theory of fluid molecules at grid points, thereby obtaining macroscopic physical quantities at the grid points, such as velocity, density, and pressure. Figure 2 The following is a schematic diagram of the evolution of fluid molecules in two dimensions. It can be seen that the fluid molecules migrate in eight directions. The calculation scale of the lattice Boltzmann is mesoscopic. Figure 2 For example, the Lattice Boltzmann method can be understood as studying the motion of the entire flow field by analyzing the eight directional motion trends of fluid molecules on each grid in the flow field. These clusters are not individual fluid molecules, but rather clusters of fluid molecules with the same motion trend. These clusters are larger than molecules but much smaller than the control volume of the macroscopic model. The evolution of the velocity distribution function of microscopic particles allows for a more intuitive description of the interactions between fluid components. Therefore, the Lattice Boltzmann method has attracted considerable attention in complex multiphase and multicomponent flow phenomena.

[0004] Among existing technical solutions, the Euler method for solving water droplet impact characteristics uses the air-water two-phase flow method (reference: Lin Guiping, Bu Xueqin, Shen Xiaobin, Yu Jia. Aircraft Icing and Anti-icing Technology [M]. Beijing: Beijing University of Aeronautics and Astronautics Press, 2016). Two-phase flow refers to the gas and liquid phases. The Euler method treats water droplets as a continuous phase and introduces the concept of droplet volume fraction to solve the continuity and momentum equations for droplet motion. The droplet volume fraction is defined as the ratio of the volume occupied by the water droplet phase in the control volume to the total control volume. The Euler method calculates on a macroscopic scale. The macroscopic model can be understood as studying the flow field motion by studying fixed finite control volumes in the flow field. The finite control volume contains a certain number of fluid molecules and their disordered motion. In icing applications, the Euler method calculates water droplet impact characteristics based on the macroscopic control volume. The goal is to simulate the macroscopic ice pattern, focusing only on the main characteristics of ice. Therefore, the Euler method does not consider the effects of intermolecular interactions or the breakup of water droplets. Furthermore, in the Euler method's droplet collection rate algorithm, the water droplet flow field and the air flow field are two relatively independent computational systems. The air flow field, after iteration, applies forces to the water droplet flow field, thereby enabling the motion and iteration of the water droplet flow field. In practice, the interaction between water droplets of a few microns and air, as well as secondary collisions and breakup of water droplets, influences the droplet impact characteristics and thus ice formation. Furthermore, the motion of air and water droplets involves collisions and interactions between identical and different molecules, often leading to phenomena such as anomalous diffusion. To capture these characteristics, it is necessary to study the droplet impact characteristics at the scale of fluid molecular motion. Microscopic computations of molecular motion are limited to a very limited number of molecules due to computer limitations, making them unsuitable for simulating the motion of large numbers of molecules. Discrete velocity methods, which are used at the mesoscopic scale for multicomponent gas motion, are complex and can only be applied to low-dimensional, simple motions, making them inappropriate for the complex flows involved in ice formation. Summary of the Invention

[0005] This paper investigates the impact of frozen water droplets on surfaces at the mesoscopic scale of fluid molecular motion and provides a novel method for calculating the impact characteristics of frozen water droplets on surfaces using the lattice Boltzmann method. This method can calculate the motion of air and water droplets in mixed air-water droplet flow fields and the impact characteristics of surface water droplets.

[0006] The technical solution of the present invention is:

[0007] A method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method, comprising the following steps:

[0008] Step 1: Obtain flow field grid data, where the flow field is the flow field around the object for which the impact characteristics of ice droplets on the surface need to be calculated;

[0009] Step 2: Set the initial flow parameters, including the initial velocity u of the air 1，0 , the initial density of air ρ 1，0 , the initial velocity of water vapor u 2，0 , the initial density of water vapor ρ 2，0 , the initial velocity u of the mixed fluid 0 ;

[0010] Calculate the initial function f of the equilibrium distribution of air components using the set initial flow field parameters α 1(eq0) and initial function f of equilibrium distribution of water vapor components α 2(eq0) ;

[0011] The air component equilibrium distribution function and the water vapor component equilibrium distribution function are used to calculate the macroscopic quantities of the air component and the water vapor component respectively: air density ρ 1 , air speed u 1 , water vapor density ρ 2 , water vapor velocity u 2 ; And use the macroscopic quantities of air components and water vapor components to calculate the macroscopic quantities of the mixed fluid: density ρ, velocity u, pressure p;

[0012] Step 3: Based on the macroscopic quantities of air components and water vapor components and the macroscopic quantities of the mixed fluid, calculate the equilibrium velocity distribution function of the air components And the equilibrium velocity distribution function of water vapor components

[0013] Step 4: The equilibrium velocity distribution function of the air components obtained in step 3 And the equilibrium velocity distribution function of water vapor components Calculate the collision term:

[0014] Air component collision term in: is the air component velocity distribution function. In the first iteration, the velocity distribution function is obtained in step 2. τ is the self-collision relaxation time;

[0015] Water vapor component collision term in: is the velocity distribution function of the water vapor component. In the first iteration, the velocity distribution function is obtained in step 2. τ 12 is the mutual collision relaxation time;

[0016] Step 5: Use the air component collision term obtained in step 4 Collision term with water vapor component According to the formula

[0017]

[0018] Calculate the velocity distribution function after the collision, including the air velocity distribution function and the water vapor velocity distribution function at non-grid points

[0019] in is the velocity distribution function of the current calculation step; is the velocity distribution function of the previous time step. When it is calculated for the first time, is the initial function of the equilibrium distribution of the corresponding component; is the collision term of the corresponding component calculated in step 4;

[0020] Then, the water vapor velocity distribution function at the non-grid point Interpolation to obtain the water vapor velocity distribution function on the grid points

[0021] Step 6: Calculate the velocity distribution function at the flow field boundary;

[0022] Step 7: Based on the velocity distribution functions of the air component and the water vapor component obtained in step 5, calculate the density, macroscopic velocity and pressure of each component, as well as the density, velocity and pressure of the mixed fluid;

[0023] Step 8: Determine whether the convergence condition is met. If so, output the calculation result of the flow field; otherwise, return to step 3.

[0024] Step 9: Using the flow field calculation results obtained in step 8, use the formula

[0025]

[0026] Calculate the local droplet collection coefficient β, where ρ 2 is the density of water vapor output in step 8, ρ ref is the reference density of water vapor, MVD is the mean water droplet diameter, MVD air is the average diameter of gas phase water particles, n is the number of water molecules, The density of air at 0 degrees Celsius.

[0027] Furthermore, the flow field grid data obtained in step 1 includes grid characteristic parameters required for realizing flow field calculation, boundary point LBM grid data files, flow field point LBM grid data files, grid unit connection relationship data files and grid unit information files.

[0028] Furthermore, in step 2, the initial function of the equilibrium distribution of air components is The formula is:

[0029]

[0030] where c s,1 is the lattice sound velocity of the air component, c 1 is the lattice velocity of the air component; To calculate the model velocity distribution, w α is the weight of velocity distribution in each direction;

[0031] Initial function of equilibrium velocity distribution of water vapor components The formula is:

[0032]

[0033] in and The distribution is the same, c s,2 is the lattice sound velocity of the water vapor component where c 2 is the lattice velocity of the water vapor component, M1 is the molecular mass of air, and M2 is the molecular mass of water vapor clusters.

[0034] Furthermore, in step 2, D2Q9 is selected for calculating the model velocity distribution.

[0035] Furthermore, in step 3,

[0036] Equilibrium velocity distribution function of air components

[0037]

[0038] Equilibrium velocity distribution function of water vapor components

[0039]

[0040] Furthermore, in step 5, according to the formula

[0041]

[0042] The water vapor velocity distribution function f on the non-grid grid point α 2 Interpolation to obtain the water vapor velocity distribution function on the grid in is the water vapor velocity distribution function after collision migration, η is the velocity distribution The ratio of the component in the x-direction to the grid size Δx.

[0043] Furthermore, in step 6, the distribution function at the boundary is processed as follows: the far-field boundary adopts the non-equilibrium extrapolation format, the layered boundary adopts the virtual layered boundary high-precision processing format, and the object surface boundary adopts the surface boundary processing method.

[0044] Beneficial effects

[0045] The present invention can solve the problem of studying the impact characteristics of surface frozen water droplets at the mesoscopic scale of fluid molecular motion; it can solve the air-water mixed flow field with fluid molecular microclusters as the research object, and capture the microscopic motion of water droplets that cannot be simulated by macroscopic methods; under complex freezing conditions, it is helpful to solve the grid scale constraint problem faced due to complex geometry when solving the flow field.

[0046] Additional aspects and advantages of the present invention will be set forth in part in the description which follows and, in part, will be obvious from the description which follows, or may be learned by practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the following description of the embodiments with reference to the accompanying drawings, in which:

[0048] Figure 1 A schematic diagram of water drop impact is shown;

[0049] Figure 2 The evolution of fluid molecules on a lattice in two dimensions is shown;

[0050] Figure 3 The interpolation of distribution functions in the differential lattice velocity format is demonstrated;

[0051] Figure 4 Describes the algorithm flow and the relationship between various models;

[0052] Figure 5 Describes the flow field calculation data flow of the embodiment;

[0053] Figure 6 The grid diagram of the cylindrical example is shown; (a) Schematic diagram of the computational domain, (b) Schematic diagram of a local grid;

[0054] Figure 7 The mesh diagram of the airfoil example is shown; (a) Schematic diagram of the computational domain, (b) Schematic diagram of a local mesh;

[0055] Figure 8 Shows the motion trajectories of air and water droplets on a cylindrical surface; (a) air (b) water droplets;

[0056] Figure 9 The calculation results of the local water droplet collection rate on the cylindrical surface are presented;

[0057] Figure 10 Shows the trajectory of air and water droplets on the airfoil surface; (a) air (b) water droplets;

[0058] Figure 11 The calculation results of the water droplet collection coefficient on the airfoil surface are shown. DETAILED DESCRIPTION

[0059] In order to study the impact characteristics of surface icing water droplets at the mesoscopic scale of fluid molecular motion, the present invention combines the algorithm characteristics of the lattice Boltzmann method and the real physical characteristics of the flow field under icing application scenarios, proposes a flow field calculation model hypothesis, constructs a flow field calculation model, establishes a calculation model of the local water droplet collection coefficient, and correspondingly proposes a numerical simulation method of water droplet impact characteristics based on the lattice Boltzmann method for cylindrical surfaces and airfoil surfaces.

[0060] First, the flow field calculation model assumptions are given:

[0061] In the application of icing scenes, the air-water droplet mixed flow field is a cloud field, in which the diameter of the water droplets is at the micron level, and the water droplets are evenly distributed throughout the flow field. At this time, it is almost impossible to track the interface between the gas phase and the liquid phase at the molecular scale. A large number of tiny water droplets are diffused in the space, and the two are more like two mutually miscible fluids. At this time, it is necessary to introduce a multi-component model to reproduce this phenomenon. At the same time, the lattice Boltzmann method has instantaneous characteristics, which can control the calculation step time step to an extremely short physical time. Taking the calculation conditions of the initial Mach number Ma = 0.1732, the grid size dx = 0.001, and the characteristic length of 1 meter as an example, every time the LBM calculation model evolves 10 million steps, the corresponding real physical world only takes one second. During this extremely short period of time, if the water droplet moves at the same speed as the incoming flow, the water droplet has only moved 0.06 microns, which is 0.6% of the diameter of the water droplet. That is to say, the water droplet has not completely left its original position, but has completed this process by deformation of its phase interface. From a macroscopic perspective, it can be approximately considered that the water droplet has hardly moved in situ. However, if this microscopic displacement process is observed from a macroscopic perspective of the overall flow field, the displacement process is too small and difficult to capture. The displacement of the water droplet cluster is a microscopic concept and cannot be accurately described on the macroscopic scale of the global flow field. On the basis of the above, in order to correctly describe the motion state of the liquid in the cloud field, based on the lattice Boltzmann method, the present invention proposes a calculation model for the single-phase multi-component mixed flow field of water vapor and air, and solves the above problems by approximating the water droplet flow field by the water flow field. The following are the assumptions of the flow field calculation model:

[0062] (1) Gas phase water tends to undergo phase transformation into liquid phase water;

[0063] In reality, gaseous water in a flow field tends to undergo a phase transition to liquid water due to its own temperature, density, or pressure. This tendency is similar to the equilibrium state in a chemical reaction: water can reach equilibrium under certain temperature, pressure, and stress conditions, forming a three-phase equilibrium system of ice, water, and vapor (reference: Liu Fei. Molecular Dynamics Simulation of Phase Transitions of Water and Methanol in Two-Dimensional Confinement [D]. University of Science and Technology of China, 2011). This is the fundamental basis for the transformation of gaseous water into liquid water in a flow field.

[0064] (2) There are sufficient equivalent water droplets distributed at each node in the flow field;

[0065] High-density water vapor is difficult to observe. This state will exist briefly and is a special transition state that can only be maintained for a short period of time. However, this state objectively exists in the physical process of the transition from gas phase water to liquid phase water. When water vapor condenses, steam (gaseous state) changes to water (liquid state) and releases latent heat at the same time (Reference: Bao Heming, Gao Shuning, Cheng Siyuan, Guan Xin. Research on the radiation mechanism of water phase change process [J]. Energy Research and Information, 2021, 37(01): 40-45. DOI: 10.13259 / j.cnki.eri.2021.01.007.). Therefore, this instantaneous flow field is used as the theoretical basis for the numerical calculation of the flow field, and the small time step of the lattice Boltzmann method is fully utilized to establish a numerical model for simulation.

[0066] (3) The equivalent water droplet does not undergo phase change during the process calculation based on the lattice Boltzmann method;

[0067] During the flow field calculation process, any equivalent water droplet is assumed to be in the gas phase at the current moment and will not undergo phase change. In other words, the equivalent water droplet will not be converted into liquid water or ice. This assumption is consistent with the Euler method (see Lee T, Lin C L. An Eulerian description of the streaming process in the lattice Boltzmann equation [J]. Journal of Computational Physics, 2003, 185(2): 445-471.). Water droplets moving in the air will not sublimate or freeze, resulting in mass or energy loss or exchange.

[0068] (4) The phase change process is short enough;

[0069] In order to ensure that the overall flow characteristics of the flow field do not change significantly, it is assumed that all phase change processes can be completed within a time length less than one flow field calculation time step. This can ensure that the water droplets complete the phase change process at the moment of contact with the surface to form the surface water film required for ice formation, and can also ensure that the flow field will not gradually shift in physical properties over time.

[0070] Combining assumptions (3) and (4), the result in the calculation process is that the initial phase of water in the flow field is always equal to the equilibrium phase of water, and the phase of water does not change before it contacts the surface, but maintains a state that tends to change its phase.

[0071] (5) The influence of air flow field on water droplet flow field is unidirectional;

[0072] Since the diameter of water droplets during the freezing process is generally in the range of 1-100 microns (reference: Yi Xian, Zhu Guolin. Numerical calculation of water droplet impact characteristics on ice surface. Collection of Aerodynamic Research, Vol. 15, 2006), it is approximately assumed that it has no effect on the air flow field. This idea is reflected in the construction process of the Euler method (reference Durst F, Miloievic D, Schnung B. Eulerian and Lagrangian predictions of particulate two-phase flows: a numerical study [J]. Applied Mathematical Modelling, 1984, 8 (2): 101-115.).

[0073] (6) There is no direct heat exchange between the water droplet flow field and the air flow field;

[0074] Since the energy exchange between the air and the water droplets during the flow process is relatively small, its impact on the overall flow field is minimal. If this process is ignored, the accuracy of the numerical calculation will also be low. Therefore, it is assumed that there is no heat exchange between the water droplets and the air (see Bragg M B. Rime ice accretion and its effect on airfoil performance / [J]. Dissertation Abstracts International, Volume: 42-07, Section: B, page: 2913. 1982.).

[0075] Based on the above assumptions, it can be concluded that the phase change of the water flow field and the lattice Boltzmann numerical simulation process are short enough. The method of approximating the water flow field to the equivalent cloud field can simulate the real physical motion through the interaction between components in single-phase multi-component flow.

[0076] Secondly, build the flow field calculation model:

[0077] In order to calculate the impact characteristics of ice droplets on surfaces in the lattice Boltzmann framework, when studying the flow of multi-component fluids in the same phase, the following four factors are mainly considered when building the flow field calculation model:

[0078] (1) In the atmosphere, the viscosity and mass of water mist and dry air are different. The model must be able to simulate flow phenomena in which the flow of components with different viscosities can affect the mixed flow.

[0079] (2) Aircraft in flight often need to face a high Reynolds number. In this state, a reasonable turbulence model must be selected.

[0080] (3) The liquid water content of the cloud field under icing application is about 1g / m 3 That is, the mass of water droplets in every cubic meter of air is 1 gram. The effect of water droplets on air flow is very small and can be ignored, so this physical property needs to be reflected in the model.

[0081] (4) In the air-water droplet mixed flow field under icing conditions, from a macroscopic perspective, the mass and inertia of water droplets are greater than those of air. From the perspective of molecular motion, the average molecular mass of air is 29, and the molecular mass of water molecules is 18. Therefore, when solving multi-component flows, actual physical phenomena must also be considered to ensure that the inertia of water vapor is greater than that of the air component.

[0082] Integrating factors 1 and 2 above, this paper selects a multi-component lattice Boltzmann model with a splitting collision model, large eddy simulation, and the Differential Lattice Speed Scheme (DLS) for the dynamic modeling of air-water mixed flow fields suitable for solving icing applications. Furthermore, addressing factor 3, this paper proposes a method for treating the equilibrium degradation of the base component, i.e., the air component. Addressing factor 4, this paper incorporates the polymerization number of water vapor molecules, combining several water vapor molecules into fluid microclusters to achieve a greater inertia than air.

[0083] For multicomponent flows, the lattice Boltzmann equation for a component i (refer to Joshi AS, Peracchio AA, Grew KN, et al. Lattice Boltzmann method for continuum, multi-component mass diffusion in complex 2D geometries[J]. Journal of Physics D: Applied Physics, 2007, 40(9): 2961) is as follows:

[0084]

[0085] in the formula is the discrete velocity of component i Velocity distribution function at position x and time t. is the discrete velocity of the particle in the αth direction. Figure 2 The model shown can obtain the discrete velocity of component i:

[0086]

[0087] where c i is the lattice velocity of component i. The differential lattice velocity format (refer to McCracken ME, Abraham J. Lattice Boltzmann methods for binary mixtures with different molecular weights [J]. Physical Review E, 2005, 71 (4): 046704) means that the lattice velocities of components with different masses are different. This method is used to solve the problem of multi-component flow with different masses. Here, it is assumed that component 1 is the lightest component, then the lattice velocities of other components c i for Differences in lattice velocities lead to different migration distances between components within a time step. Therefore, the velocity distribution of a component at a grid point will migrate to non-grid nodes. Here, interpolation is used to obtain the distribution function of the grid point, which is the core of the differential lattice velocity format. Figure 3 is a schematic diagram of the distribution function interpolation, with Figure 2 Taking the velocity component f5 of the model as an example, its migration process is as follows Figure 3 The distribution function at the grid point O is obtained from the migrated distribution function, where η is the ratio of the component of the velocity e5 in the x direction to the grid scale Δx. It can be calculated as follows:

[0088]

[0089] The collision term on the right side of equation (1) is expanded as follows:

[0090]

[0091] The item on the right represents the self-collision item, which is approximated by BGK. The calculation formula for this item is:

[0092]

[0093] In the formula, τ is the self-collision relaxation time, is the lattice sound velocity. The equilibrium distribution function It is given by the following formula:

[0094]

[0095] Among them Figure 2 For example, the weight w0=4 / 9, w 1-4 =1 / 9, w 5-8 =1 / 36. According to the distribution function, the density of each component ρ i and equilibrium speed u i(eq) It can be calculated by the following formula:

[0096]

[0097]

[0098] According to the Stefan-Maxwell model, the total density ρ and macroscopic velocity u satisfy the following formula:

[0099]

[0100]

[0101] In order to improve the calculation of Reynolds number, large eddy simulation is introduced. After adding Smagorinsky eddy viscosity model, the eddy viscosity coefficient ν t It can be expressed as:

[0102]

[0103]

[0104] Among them S c is the Smagorinsky constant, which is taken as 0.16 here; Δ is the filter scale, which is related to the grid scale; τ0 is the laminar relaxation time of the mixed gas, which can be obtained by τ0 = ν0 / c s,1 2Δt+1 / 2 is calculated, ν0 is the kinematic viscosity of the gas mixture; is the filtered average momentum flux, the second moment of the non-equilibrium velocity distribution It can be obtained by the following formula:

[0105]

[0106] Turbulent relaxation time τ t Equal to 3ν t , so the relaxation time in formula (5) can be expressed as:

[0107] τ=τ0+τ t (14)

[0108] The mutual collision term in equation (4) is calculated as follows:

[0109]

[0110] For the gas-liquid two components in this paper, the mutual collision relaxation time τ ij It is given by the following formula:

[0111]

[0112] Among them, M i 、M j is the molecular mass of each component, n is the total molar coefficient, D ij is the diffusion coefficient,

[0113] p=ρ1c s,1 2 +ρ2c s,2 2 is the total pressure.

[0114] Considering that the impact of water droplets on air is negligible, the present invention performs equilibrium degradation on the air components. The equilibrium degradation method is: the mutual collision term in the equilibrium equation of the air components is zero, and the equation is degenerated into a single-component flow equilibrium equation. The velocity u of the mixed fluid in the equilibrium function is equal to the velocity of the air component.

[0115] Therefore, the specific expression of the equilibrium degenerate lattice Boltzmann equation of the air component is:

[0116]

[0117] Equilibrium distribution function f α i(eq) for:

[0118]

[0119] Then the local water droplet collection coefficient calculation model is established:

[0120] In the framework of the lattice Boltzmann model, water vapor molecules are taken as the research object. The expression of the local water droplet collection coefficient β that characterizes the water droplet impact characteristics in the present invention is:

[0121]

[0122] The above formula is the ratio of the actual amount of water collected on the micro-element surface to the maximum possible amount of water collected on the micro-element surface. Therefore, it is a parameter that characterizes the water collection ability of the micro-element surface. water is the density of water vapor, ρ ref is the reference density of water vapor, MVD is the mean water droplet diameter, MVD air is the average diameter of gas phase water particles, LWC is the liquid water content, n is the polymerization number of water molecules, The density of air at 0 degrees Celsius.

[0123] In summary, the algorithm flow and the relationship between each model are as follows: Figure 4 As shown in the figure, the first step of flow field initialization can obtain the initial macroscopic quantities and velocity distribution function; the second step is to calculate the equilibrium distribution function based on the initial macroscopic quantities and velocity distribution function; the third step is to calculate the collision term, which is the reason for the change of the velocity distribution function; the fourth step is the migration and interpolation of the distribution function. The value of the distribution function is not changed in this process. It only migrates on the grid nodes, and for the components that deviate from the grid nodes, the distribution function on the nodes is solved by interpolation.

[0124] The present invention is further explained below by taking the calculation of water droplet impact characteristics on cylindrical surfaces and airfoil surfaces as an example:

[0125] Example 1: Numerical simulation of water droplet impact characteristics on cylindrical surface:

[0126] To calculate the impact characteristics of water droplets on a surface, the motion state of the flow field must first be obtained. In the Euler method, based on the macroscopic control body, the motion state of the air is first solved, and then applied to the water droplets as a known condition to solve the continuity equation and momentum equation of the water droplets. Finally, the volume distribution and velocity distribution of the water droplets are obtained, and then the impact characteristics of the water droplets are obtained. In the present invention, based on the fluid micro-clusters at a mesoscopic scale that are much smaller than the macroscopic control body, by solving the multi-component lattice Boltzmann equation with large eddy simulation, split collision model, and differential lattice velocity format, the equilibrium state of the air components is degraded and the polymerization number of water vapor molecules is introduced to make the model more consistent with physical reality, thereby obtaining the macroscopic physical quantities of air and water droplets, and then obtaining the impact characteristics of the water droplets.

[0127] Table 1 shows the calculation conditions for the cylindrical example

[0128] Table 1 Calculation conditions for cylindrical example

[0129] V(m / s) AOA(°) MVD(um) <![CDATA[LWC(g / m 3 )]]> Ts(℃ / K) Reference conditions 31.0 0 20 1.0 -15 / 258.3

[0130] This example is referenced from (Chen Jinping, Ice wind tunnel test study on two-dimensional cylindrical icing and calculation of water droplet impact characteristics [D]. Nanjing University of Aeronautics and Astronautics, 2013). Figure 5 The data flow diagram of the calculation process is shown in the figure. Figure 5 The dashed line represents the flow of water vapor data, and the dotted line represents the flow of air data. The data that needs external input during the entire flow field calculation process include grid information and initial flow distribution parameters. Figure 6 The grid diagram of the cylindrical example is shown in the figure. As shown in the figure, the grid adopts a non-uniform adaptive Cartesian grid, where the layer boundary is the intersection of the non-uniform grid, the surface boundary is the intersection of the grid and the surface, and the far field boundary is the boundary of the external flow field. The initial information of the fluid (air and water vapor) in the grid is generated at the same time. Figure 4 and Figure 5 The specific calculation steps are as follows:

[0131] Step (1) inputs the flow field mesh data for the cylindrical model. The mesh data includes the mesh characteristic parameters required for flow field calculation, boundary point LBM mesh data files, flow field point LBM mesh data files, mesh cell connection relationship data files, and all mesh cell information files. The above mesh information can determine the number of mesh nodes, mesh boundary information and geometry information, and the type of mesh points (boundary points or flow field points).

[0132] Step (2), flow field initialization. The parameters of flow field initialization are given manually, including: the initial velocity u of the air 1，0 , the initial density of air ρ 1，0 , the initial velocity of water vapor u 2，0 , the initial density of water vapor ρ 2，0 , the initial velocity u of the mixed fluid 0 The equilibrium distribution function of the entire basin is solved using the following formula using the initial parameters, and the macroscopic quantities are calculated using the equilibrium distribution function to complete the flow field initialization.

[0133] Solve the equilibrium state degenerate form of the air component in each velocity direction equilibrium velocity distribution function f α 1(eq0) :

[0134]

[0135] When initialized: u 1，0 is the initial velocity of the air; ρ 1，0 is the initial density of air; c s,1 is the lattice sound velocity of the air component, c1 is the lattice velocity of the air component, generally taken as 1. To calculate the model velocity distribution, w α is the weight of velocity distribution in each direction, and the velocity distribution model is D2Q9, such as Figure 2 As shown, w α The expression is as follows:

[0136]

[0137] w0=4 / 9,w 1-4 =1 / 9, w 5-8 =1 / 36.

[0138] Equilibrium velocity distribution function of water vapor components The calculation is as follows:

[0139]

[0140] in w α Calculate the same w α ;u 2，0 is the initial velocity of water vapor; ρ 2，0 is the initial density of water vapor; c s,2 is the lattice sound velocity of the water vapor component in is the lattice velocity of the water vapor component, M1 is the molecular mass of air, and M2 is the molecular mass of the water vapor cluster; u 0 is the initial velocity of the mixed fluid.

[0141] Equilibrium distribution function of air and water vapor and Solve the macroscopic quantities of air and water components for iterative calculation, density ρ i , speed u i , pressure P i It can be calculated by the following formula:

[0142]

[0143]

[0144]

[0145] The macroscopic quantities of the mixed fluid: density ρ, velocity u, and pressure p, can be obtained by the following formula:

[0146]

[0147]

[0148] p=ρ 1 c s,1 2 +ρ 2 c s,2 2

[0149] After the initialization process, the macroscopic quantities and equilibrium distribution functions required for subsequent iterative settlement are obtained.

[0150] Step (3), evolution and iterative calculation, calculate the equilibrium distribution function. The calculation steps and formula of the equilibrium distribution function are the same as step (2), where the required macroscopic quantity is: air density ρ 1 , air speed u 1 , water vapor density ρ 2 , water vapor velocity u 2 , the velocity u of the mixed fluid, which is the macroscopic quantity calculated by step (2). Based on the above formula of step (2), the equilibrium velocity distribution function of the air component used for step (4) is obtained: And the equilibrium velocity distribution function of water vapor components

[0151] Equilibrium velocity distribution function of air components in each velocity direction

[0152]

[0153] Equilibrium velocity distribution function of water vapor components The calculation is as follows:

[0154]

[0155] Step (4), calculate the collision term. The equilibrium distribution function obtained from step (3) is The collision term is calculated by the following formula. The collision is the reason for the change of the velocity distribution function.

[0156] Collision term for air components The calculation is as follows:

[0157]

[0158] in: is the velocity distribution function of the air component. Since the collision has not occurred at the time of initialization, the velocity distribution function is the equilibrium velocity distribution function obtained in step (2) during the first iteration. is the equilibrium distribution function obtained in step (3). τ is the self-collision relaxation time.

[0159] For the calculation of collision terms of water vapor components:

[0160]

[0161] in is the velocity distribution function of the water vapor component. Since the collision has not occurred at the time of initialization, the velocity distribution function is the equilibrium velocity distribution function obtained in step (2) during the first iteration. is the equilibrium velocity distribution function obtained in step (3). τ is the self-collision relaxation time, τ 12 is the mutual collision relaxation time. u is the velocity of the mixed flow, which is obtained by initialization during the first calculation; u 2 、u 1 are the velocities of water vapor and air, obtained by initialization during the first iteration.

[0162] Step (5), calculate the velocity distribution function of the air and the water vapor velocity distribution function at non-grid points Using the air component collision term obtained in the previous step Collision term with water vapor component The velocity distribution function after the collision is calculated using the following formula.

[0163]

[0164] in is the velocity distribution function of the current calculation step; is the velocity distribution function of the previous time step. When the calculation is for the first time, this value is the initial equilibrium distribution function, that is, the air component has The components of water vapor are is the collision term calculated in step (4).

[0165] Solve the distribution function of water vapor velocity at the grid points It is composed of the water vapor velocity distribution function on the non-grid grid point Interpolation is obtained, the schematic diagram is as follows Figure 3 :

[0166]

[0167] in is the water vapor velocity distribution function after collision migration, η is the velocity distribution The ratio of the component in the x-direction to the grid size Δx.

[0168] Step (6), calculate the distribution function at the boundary. The calculation in step (5) obtains the velocity distribution function of all flow field points, and then calculates the velocity distribution function at the boundary. Processing of distribution functions at boundaries: the far-field boundary adopts the non-equilibrium extrapolation format (refer to Guo ZL, Zheng CG, Shi B C. Non-equilibrium extrapolation method for velocity and boundary conditions in the lattice Boltzmann method. Chinese Physics, 2002, 11(4): 0366-0374.); the layered boundary adopts the virtual layered boundary high-precision processing format (refer to Northwestern Polytechnical University. A high-precision processing method for virtual layered boundaries in the tree grid lattice Boltzmann method: CN202210219613.9[P]. 2022-06-07.); the object surface boundary adopts the surface boundary processing method (refer to Northwestern Polytechnical University. Processing method for surface boundaries in numerical simulation of lattice Boltzmann method: CN202210219614.3[P]. 2022-06-07.).

[0169] Step (7), calculate the macroscopic quantity. The velocity distribution function of the air and water vapor components obtained in step (5) is and To find the density of each component ρ i and macroscopic speed u i and pressure p i , calculated by the following formula:

[0170]

[0171]

[0172]

[0173] The macroscopic quantity of the mixed fluid is obtained by the following formula:

[0174]

[0175]

[0176]

[0177] Where ρ is the density of the mixed flow, u is the velocity of the mixed flow, p is the pressure of the mixed flow, and u i is the macroscopic velocity of a component, ρ i is the macroscopic density of component i, c s,i Lattice sound velocity of component i.

[0178] Step (8) After reaching the convergence standard, the calculated data of the flow field and the resulting image processing are output, including the calculated results of density, velocity, and pressure.

[0179] Step (9) calculates the local water droplet collection coefficient. The local water droplet collection coefficient can be expressed in the formula as follows:

[0180]

[0181] Where β is the local droplet collection coefficient, ρ 2 is the density of water vapor, ρ ref is the reference density of water vapor, where the initial density, MVD is the average water droplet diameter, and MVD air is the average diameter of gas phase water particles (taken as 0.0004), n is the polymerization number of water molecules (taken as 80.556), The density of air at 0 degrees Celsius.

[0182] Figure 8 The motion trajectories of air (a) and water droplets (b) for a numerical simulation of water droplet impact characteristics on a cylindrical surface are presented. The figure shows that the streamlines of the air bend sharply as it approaches the cylindrical surface, with the airflow sticking to the upper and lower surfaces, exhibiting obvious flow-around characteristics. At the trailing edge, the airflow recovers. At the same time, the collection characteristics of the water droplets can be clearly seen. Due to their large mass and inertia, they are not prone to changing their state of motion. Therefore, the water droplets cannot bypass the leading edge of the cylinder and instead impact the wall. Due to the collection characteristics of the water droplets, a water droplet shielding area appears at the trailing edge, meaning no water droplets exist in this area.

[0183] Figure 9 Comparison results of the local water droplet collection coefficient on a cylindrical surface are presented. It can be seen that the local water droplet collection coefficient ranges from 0 to 1, is maximized at the leading edge stagnation point, and gradually decreases along the upper and lower surfaces until it reaches zero at the droplet impact limit. The calculation results of the present invention are consistent with those calculated using the Euler method.

[0184] Example 2: Numerical simulation of water droplet impact characteristics on airfoil surface:

[0185] Table 2 shows the calculation conditions of the airfoil example

[0186] Table 2 Calculation conditions for airfoil example

[0187]

[0188] This example is based on (AIAA. Results of an icing test on a NACA0012 airfoil in the NASA Lewis Icing Research Tunnel-30th Aerospace Sciences Meeting and Exhibit (AIAA). 1992). When performing the numerical simulation of the water droplet impact characteristics on the airfoil surface, the calculation model and calculation process used are the same as those in Example 1, where Figure 7 Schematic diagram of the mesh for the airfoil example. Figure 10 The motion trajectories of air (a) and water droplets (b) are given in the numerical simulation example of water droplet impact characteristics on the airfoil surface. It can also be found that the water droplets impact and collect on the leading edge of the airfoil. Figure 11 The comparison chart of the local water droplet collection rate of the airfoil example shows that the calculation results are generally in good comparison with the experimental results.

[0189] Although the embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and are not to be construed as limitations on the present invention. A person skilled in the art may change, modify, replace and modify the above embodiments within the scope of the present invention without departing from the principles and purpose of the present invention.

Claims

1. A method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method, characterized by: The following steps are involved: Step 1: Obtain flow field grid data, where the flow field is the flow field around the object for which the impact characteristics of ice droplets on the surface need to be calculated; Step 2: Set the initial flow parameters, including the initial velocity u of the air 1，0 , the initial density of air ρ 1，0 , the initial velocity of water vapor u 2，0 , the initial density of water vapor ρ 2，0 , the initial velocity u of the mixed fluid 0 ; Calculate the initial function of the equilibrium distribution of air components using the set initial flow field parameters and initial function of equilibrium distribution of water vapor components The air component equilibrium distribution function and the water vapor component equilibrium distribution function are used to calculate the macroscopic quantities of the air component and the water vapor component respectively: air density ρ 1 , air speed u 1 , water vapor density ρ 2 , water vapor velocity u 2 ; The macroscopic quantities of the mixed fluid are calculated using the macroscopic quantities of the air component and the water vapor component: density ρ, velocity u, and pressure p; Step 3: Based on the macroscopic quantities of air components and water vapor components and the macroscopic quantities of the mixed fluid, calculate the equilibrium velocity distribution function of the air components and the equilibrium velocity distribution function of water vapor components Step 4: The equilibrium velocity distribution function of the air components obtained in step 3 and the equilibrium velocity distribution function of water vapor components Calculate the collision term: Air component collision term in: is the air component velocity distribution function. In the first iteration, the velocity distribution function is obtained in step 2. τ is the self-collision relaxation time; Water vapor component collision term in: is the velocity distribution function of the water vapor component. In the first iteration, the velocity distribution function is obtained in step 2. τ 12 is the mutual collision relaxation time; Step 5: Use the air component collision term obtained in step 4 Collision term with water vapor component According to the formula Calculate the velocity distribution function after the collision, including the air velocity distribution function and the water vapor velocity distribution function at non-grid points in is the velocity distribution function of the current calculation step; is the velocity distribution function of the previous time step. When it is calculated for the first time, is the initial function of the equilibrium distribution of the corresponding component; is the collision term of the corresponding component calculated in step 4; Then, the water vapor velocity distribution function at the non-grid point Interpolation to obtain the water vapor velocity distribution function on the grid Step 6: Calculate the velocity distribution function at the flow field boundary; Step 7: Based on the velocity distribution functions of the air component and the water vapor component obtained in step 5, calculate the density, macroscopic velocity and pressure of each component, as well as the density, velocity and pressure of the mixed fluid; Step 8: Determine whether the convergence condition is met. If so, output the calculation result of the flow field; otherwise, return to step 3. Step 9: Using the flow field calculation results obtained in step 8, use the formula Calculate the local droplet collection coefficient β, where ρ 2 is the density of water vapor output in step 8, ρ ref is the reference density of water vapor, MVD is the mean water droplet diameter, MVD air is the average diameter of gas phase water particles, n is the number of water molecules, The density of air at 0 degrees Celsius.

2. The method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method according to claim 1, characterized in that: The flow field grid data obtained in step 1 includes the grid characteristic parameters required for realizing flow field calculation, boundary point LBM grid data files, flow field point LBM grid data files, grid unit connection relationship data files and grid unit information files.

3. The method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method according to claim 1, characterized in that: In step 2, the initial function of the equilibrium distribution of air components is The formula is: where c s,1 is the lattice sound velocity of the air component, c 1 is the lattice velocity of the air component; To calculate the model velocity distribution, w α is the weight of velocity distribution in each direction; Initial function of equilibrium velocity distribution of water vapor components The formula is: in and The distribution is the same, c s,2 is the lattice sound velocity of the water vapor component where c 2 is the lattice velocity of the water vapor component, M 1 is the molecular mass of air, M 2 is the molecular mass of water vapor clusters.

4. The method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method according to claim 3, characterized in that: In step 2, select D2Q9 to calculate the model velocity distribution.

5. The method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method according to claim 2, characterized in that: In step 3, Equilibrium velocity distribution function of air components Equilibrium velocity distribution function of water vapor components 6. The method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method according to claim 1, characterized in that: In step 5, according to the formula The water vapor velocity distribution function at the non-grid point Interpolation to obtain the water vapor velocity distribution function on the grid in is the water vapor velocity distribution function after collision migration, η is the velocity distribution The ratio of the component in the x-direction to the grid size Δx.

7. The method for calculating the impact characteristics of ice droplets on a surface based on the lattice Boltzmann method according to claim 1, characterized in that: In step 6, the distribution function at the boundary is processed as follows: the far-field boundary adopts the non-equilibrium extrapolation format, the layered boundary adopts the virtual layered boundary high-precision processing format, and the object surface boundary adopts the surface boundary processing method.

Images