CN113779904B - Icing phase change calculation method based on coupling of continuous liquid film and discrete liquid film - Google Patents

Icing phase change calculation method based on coupling of continuous liquid film and discrete liquid film Download PDF

Info

Publication number
CN113779904B
CN113779904B CN202111084586.0A CN202111084586A CN113779904B CN 113779904 B CN113779904 B CN 113779904B CN 202111084586 A CN202111084586 A CN 202111084586A CN 113779904 B CN113779904 B CN 113779904B
Authority
CN
China
Prior art keywords
liquid film
liquid
icing
unit
discrete
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111084586.0A
Other languages
Chinese (zh)
Other versions
CN113779904A (en
Inventor
王强
易贤
刘宇
任靖豪
李维浩
石雅楠
周靓
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Low Speed Aerodynamics Institute of China Aerodynamics Research and Development Center
Publication of CN113779904A publication Critical patent/CN113779904A/en
Application granted granted Critical
Publication of CN113779904B publication Critical patent/CN113779904B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention is suitable for the technical field of aircraft deicing, and provides an icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film, which comprises modeling a flow field and constructing a water drop motion equation; tracking the track of liquid drops, and calculating the area of each fused liquid film unit by fusing the different liquid film units into a larger circular liquid film unit when different liquid film units meet; traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unit critical Is a relationship of (2); when S is>S critical Calculating the icing thickness by adopting a continuous liquid film icing phase calculation method; when S is<S critical And calculating the icing thickness by adopting a discrete liquid film icing phase calculation method, updating the solid wall surface, and finally obtaining the simulated ice shape of the complex structure. Compared with the traditional icing phase change calculation method based on coupling of the continuous liquid film and the discrete liquid film, the icing phase change calculation method based on coupling of the continuous liquid film and the discrete liquid film only adopts continuous liquid film hypothesis to calculate, so that the simulated ice shape is more accurate, and accurate simulation of the complex ice shape can be realized.

Description

Icing phase change calculation method based on coupling of continuous liquid film and discrete liquid film
Technical Field
The invention relates to the technical field of aircraft deicing, in particular to an icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film.
Background
Icing is an important contributor to the safety of aircraft flight. The aerocraft is extremely easy to generate icing phenomenon when crossing the cloud layer containing supercooled water drops. Ice accumulation attached to the surfaces of core components such as wings not only can reduce the lift force and increase the resistance of an airplane, but also can cause important potential safety hazards such as pitching and course stability loss of the airplane. At present, the icing research means of the aircraft mainly comprise a real flight test, a spraying machine flight test, a ground icing wind tunnel test, icing numerical calculation and the like, wherein the icing numerical calculation has the advantages of rapidness, high efficiency, economy and the like, and is an important means which cannot be obtained by the current icing safety evaluation of the aircraft.
At present, the icing numerical calculation flow field mainly comprises four aspects of flow field calculation, liquid water collection calculation, icing phase change calculation and grid deformation calculation. The icing phase transformation calculation is to realize the simulation of the growth of the ice accumulation on the surface of the aircraft through the calculation of the heat and mass transfer process of the liquid film on the surface of the aircraft after the liquid water collection characteristic of the surface of the aircraft is obtained. The icing phase change calculation model of the current main stream mainly comprises: a messenger model, a Myers model, a ShallowWater model, and the like. The Messinger model simplifies the icing process of the aircraft surface into a simple mass and energy conservation equation, and solves the problem by algebraic equation calculation of a solid wall surface unit; the Myers model further considers the normal energy transport of the liquid film and the ice layer on the basis of the Messinger model by introducing the Stefen condition, so that the calculation accuracy is improved; from the point of numerical calculation, the Shallow Water model develops a more comprehensive liquid film transportation equation, realizes the simulation of the liquid film icing process by using compatible conditions, and is successfully applied to ANSYS FENSAP-ICE software. The formed icing phase change calculation model provides important conditions for icing calculation of the aircraft.
Despite the rapid development of current icing phase change computational models, accurate simulation of highly complex ice shapes, such as three-dimensional shell ice typical of swept wings, is still not possible.
Disclosure of Invention
In order to solve the technical problem that the accurate simulation of the highly complex ice shape cannot be realized in the prior art, the invention provides an ice shape simulation method, which is an icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film.
The invention discovers in the long-term research process that the existing icing phase change model is based on continuous liquid film hypothesis, so that the liquid film behavior under the characteristic condition of highly complex discrete ice shape can not be accurately represented, and the actual complex ice shape has both continuous liquid film and discrete liquid film. Therefore, the invention breaks the traditional continuous liquid film assumption, calculates the liquid water collection characteristic based on the Monte Carlo method (Monte Carlo Method), and converts the discrete supercooled water drop particle packet into a liquid film unit which can move discretely on the solid wall surface. And on the basis of tracking the track of the discrete liquid film unit, calculating the heat and mass transfer process inside the liquid film unit, and finally realizing the simulation of the complex icing phase change process.
An icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film comprises the following steps:
1-1, setting a calculated time starting point t=0, wherein t is time;
1-2, modeling a flow field, constructing a water drop motion equation, tracking a water drop track, and calculating a liquid water collection coefficient beta;
1-3, tracking the movement process of a discrete liquid film unit formed after liquid drops strike a wall surface, and calculating the area S of each fused liquid film unit by fusing the liquid film units into a larger circular liquid film unit when different liquid film units meet;
1-4, traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unit critical Is a relationship of (2);
when S is>S critical Calculating the icing thickness by adopting a continuous liquid film icing phase calculation method; and updating the fixed wall surface;
when S is<S critical By usingCalculating the icing thickness by a discrete liquid film icing phase calculation method; and updating the fixed wall surface;
1-5 let t=t+1, if t<t end Returning to the step 1-2; if t=t end And (5) ending the calculation.
Further, in the step 1-3, the method for calculating the area S of the fusion liquid film unit is as follows:
Figure BDA0003264620850000031
wherein m is add To fuse the mass sum of the liquid film, h ave To average thickness of the fused liquid film ρ w Is the density of the liquid film.
Further, in the step 1-5, the continuous liquid film icing phase calculation method is a mass conservation equation, a momentum conservation equation and an energy conservation equation of a simultaneous continuous liquid film icing model, and the liquid film icing mass is solved
Figure BDA0003264620850000041
Further obtaining the icing thickness h of the fusion liquid film unit ice ，/>
Figure BDA0003264620850000042
ρ ice Is ice density.
Further, the mass conservation equation is:
Figure BDA0003264620850000043
Figure BDA0003264620850000044
h is a wall fixing unit C S Is used for the liquid film thickness of the liquid crystal display panel,
Figure BDA0003264620850000045
for the average speed in the liquid film, +.>
Figure BDA0003264620850000046
Is a wall fixing unit C S Boundary normal vector of S Cs Is C S L is C S Boundary line length, < >>
Figure BDA0003264620850000047
For the mass of the droplets striking the wall-fixing unit, +.>
Figure BDA0003264620850000048
Evaporation mass for liquid film and->
Figure BDA0003264620850000049
For the icing quality of the liquid film, the subscript i takes a solid wall unit C S η is the normal coordinate in the boundary layer of the liquid film;
Figure BDA00032646208500000410
Figure BDA00032646208500000411
Figure BDA0003264620850000051
wherein U is ∞ LWC is the liquid water content, h, for the free incoming wind speed c R is a gas constant, C is a convective heat transfer coefficient p To fix the pressure specific heat capacity, L e Is the number of Lewis, p sat Is saturated vapor pressure, T s And T ∞ Respectively the surface of the liquid film and the infinite incoming flow temperature,
Figure BDA0003264620850000052
is a wall fixing unit C S Boundary normal vector of boundary surface i, hi is fixed wall unit C S Liquid film thickness of boundary surface i>
Figure BDA0003264620850000053
Is a wall fixing unit C S Liquid film flow velocity of boundary surface i, L i Is a wall fixing unit C S Boundary length of boundary surface i;
the conservation of momentum equation is:
Figure BDA0003264620850000054
Figure BDA0003264620850000055
Figure BDA0003264620850000056
Figure BDA0003264620850000057
Figure BDA0003264620850000058
Figure BDA0003264620850000059
wherein ζ is a shape factor, p, generated by curve coordinate transformation w Is the pressure in the liquid film, and is the pressure in the liquid film,
Figure BDA0003264620850000061
gravitational acceleration vector, < >>
Figure BDA0003264620850000062
And->
Figure BDA0003264620850000063
Shear forces acting on the free and bottom surfaces of the liquid film, respectively,/-, are applied to the liquid film>
Figure BDA0003264620850000064
For vector cross-multiplying symbols, p a Is the ambient pressure, σ is the surface tension coefficient, +.>
Figure BDA0003264620850000065
Is a no-pull operator; A. b, C is an intermediate variable.
The energy conservation equation is:
Figure BDA0003264620850000066
wherein the method comprises the steps of
Figure BDA0003264620850000067
Kinetic energy and internal energy for impinging droplets, +.>
Figure BDA0003264620850000068
Is the convective heat transfer between the liquid film and the airflow field,
Figure BDA0003264620850000069
Energy carried away by liquid film evaporation, +.>
Figure BDA00032646208500000610
Latent heat generated for freezing liquid film, +.>
Figure BDA00032646208500000611
And->
Figure BDA00032646208500000612
The energy input item and the energy output item related to the liquid film transportation are specifically shown as follows:
Figure BDA00032646208500000613
Figure BDA00032646208500000614
subscripts in and out characterize the liquid film inflow and outflow boundaries of the solid wall unit, respectively.
1. The method for calculating icing phase change based on coupling of continuous liquid film and discrete liquid film according to claim 4, wherein in step 1-6, the method for calculating icing phase of discrete liquid film is to solve icing quality of liquid film by combining the following equations
Figure BDA0003264620850000071
Further obtaining the icing thickness h of the fusion liquid film unit ice ；/>
Figure BDA0003264620850000072
Figure BDA0003264620850000073
F＝F contact +F between +F impingement
Figure BDA0003264620850000074
Wherein F is contact The solid-liquid-gas three-phase line acting force of discrete liquid film F between For discrete liquid film interaction force, F impingement The water drops impact the discrete liquid film to generate acting force, F is resultant force and Cp is Cp w Is the constant pressure specific heat capacity of water, T is the temperature of a discrete liquid film,
Figure BDA0003264620850000075
is the mutual energy transport between the discrete liquid films.
Further, in the step 1-2, the equation of motion of the water droplets is:
Figure BDA0003264620850000076
wherein m is d Is the mass of liquid drop,
Figure BDA0003264620850000077
Is a droplet position vector, t is time, ρ w 、ρ a Respectively the liquid film and the air density, V d For the volume of the drop>
Figure BDA0003264620850000078
Gravitational acceleration vector, < >>
Figure BDA0003264620850000081
For the air flow speed>
Figure BDA0003264620850000082
For droplet velocity, A d Is the cross-sectional area of the liquid drop, C d For the droplet drag coefficient, the following empirical solution is followed:
Figure BDA0003264620850000083
Figure BDA0003264620850000084
wherein Re is the Reynolds number of the liquid drop, re d Is the relative reynolds number of the droplet.
Further, the calculation of the liquid water collection coefficient beta is as follows:
Figure BDA0003264620850000085
wherein the method comprises the steps of
Figure BDA0003264620850000086
Is infinitely distant mass flow flux, < >>
Figure BDA0003264620850000087
Is impacted on the fixed wall unit C in unit time S The liquid water mass flux contained by the upper droplet particle package,
Figure BDA0003264620850000088
Figure BDA0003264620850000089
/>
wherein the method comprises the steps of
Figure BDA0003264620850000091
For striking against the fixed wall unit C S The mass of the liquid drops in the wall fixing unit, N is the total number of the liquid drops which strike the wall fixing unit in dt time, j is the number of the liquid drops, S Cs Is a wall fixing unit C S Area of U ∞ Is the free incoming flow velocity.
Compared with the traditional method for calculating by adopting continuous liquid film hypothesis, the icing phase change calculation method based on the coupling of the continuous liquid film and the discrete liquid film provided by the invention has the advantages that the simulated ice shape is more accurate, and the accurate simulation of the complex ice shape can be realized; the method can provide a powerful calculation way for simulating the sweepback wing discrete shell ice.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following description will briefly explain the embodiments of the present invention or the drawings used in the description of the prior art, and it is obvious that the drawings described below are only some embodiments of the present invention, and other drawings can be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of an icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film according to an embodiment of the present invention.
Detailed Description
The following description provides many different embodiments, or examples, for implementing different features of the invention. The elements and arrangements described in the following specific examples are presented for purposes of brevity and are provided only as examples and are not intended to limit the invention.
The icing phase change calculation method based on coupling of the continuous liquid film and the discrete liquid film comprises the following steps of:
1-1, setting a time starting point t, wherein t=0;
1-2, modeling a flow field, constructing a water drop motion equation, tracking a water drop track, and calculating a liquid water collection coefficient beta;
the present invention introduces the concept of a supercooled water droplet particle packet, i.e., it is believed that the typical characteristic supercooled water droplet dynamics in a particle packet can be representative of the water droplet properties within the particle packet. In order to obtain the motion trail of the particle packet, a water drop motion equation is constructed as follows:
Figure BDA0003264620850000101
wherein m is d Is the mass of liquid drop,
Figure BDA0003264620850000102
Is a droplet position vector, t is time, ρ w 、ρ a Respectively the liquid film and the air density, V d For the volume of the drop>
Figure BDA0003264620850000103
Gravitational acceleration vector, < >>
Figure BDA0003264620850000104
For the air flow speed>
Figure BDA0003264620850000105
For droplet velocity, A d Is the cross-sectional area of the liquid drop, C d For the droplet drag coefficient, the following empirical solution is followed:
Figure BDA0003264620850000106
Figure BDA0003264620850000107
wherein Re is the Reynolds number of the liquid drop, re d Is the relative reynolds number of the droplet; on the basis, the liquid water quality brought by different particle packages collected in a unit time in a solid wall unit is counted, the liquid water collection characteristic is obtained, and the liquid water collection coefficient beta is further obtained by calculating:
Figure BDA0003264620850000111
wherein, the liquid crystal display device comprises a liquid crystal display device,
Figure BDA0003264620850000112
is infinitely distant mass flow flux, < >>
Figure BDA0003264620850000113
Is impacted on the fixed wall unit C in unit time S The liquid water mass flux contained by the upper droplet particle package,
Figure BDA0003264620850000114
Figure BDA0003264620850000115
wherein the method comprises the steps of
Figure BDA0003264620850000116
For striking against the fixed wall unit C S The mass of the liquid drops in the wall fixing unit, N is the total number of the liquid drops which strike the wall fixing unit in dt time, j is the number of the liquid drops, S Cs Is a wall fixing unit C S Area of U ∞ Is the free incoming flow velocity.
1-3, tracking the movement process of a discrete liquid film unit formed after a liquid drop impacts a wall surface by utilizing a water drop movement equation, and supposing that the mass of a particle packet impacting on a solid wall unit is so small that a round unit attached to the solid wall can be considered to be formed after the solid wall is impacted, so that the liquid drop impacts on the solid wall to form the discrete liquid film, and when different discrete liquid films are fused, the discrete liquid film is converted into a continuous liquid film.
Assuming that different discrete liquid films meet, the discrete liquid films are fused into larger circular liquid film units, the center positions of the fused liquid film units are located at the average positions of the centroids of the different liquid film units, the liquid film mass is the sum of the mass of each liquid film before meeting, and the area S of each fused liquid film unit is calculated:
Figure BDA0003264620850000121
wherein m is add To fuse the mass sum of the liquid film, h ave To average thickness of the fused liquid film ρ w Is the density of the liquid film.
1-4, traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unit critical Is a relationship of (2);
when S is>S critical Calculating the icing thickness by adopting a continuous liquid film icing phase calculation method; and updating the fixed wall surface;
the continuous liquid film icing phase calculation method is to solve the liquid film icing quality by using a mass conservation equation, a momentum conservation equation and an energy conservation equation of a simultaneous continuous liquid film icing model
Figure BDA0003264620850000122
Further obtaining the icing thickness h of the fusion liquid film unit ice ，
Figure BDA0003264620850000123
ρ ice Is ice density.
Wherein, the mass conservation equation is:
Figure BDA0003264620850000131
Figure BDA0003264620850000132
h is a wall fixing unit C S Is used for the liquid film thickness of the liquid crystal display panel,
Figure BDA0003264620850000133
for the average speed in the liquid film, +.>
Figure BDA0003264620850000134
Is a wall fixing unit C S Boundary normal vector of S Cs Is C S L is C S Boundary line length, < >>
Figure BDA0003264620850000135
For the mass of the droplets striking the wall-fixing unit, +.>
Figure BDA0003264620850000136
Evaporation mass for liquid film and->
Figure BDA0003264620850000137
For the icing quality of the liquid film, the subscript i takes a solid wall unit C S η is the normal coordinate in the boundary layer of the liquid film;
Figure BDA0003264620850000138
Figure BDA0003264620850000139
Figure BDA00032646208500001310
wherein U is ∞ For free-running wind speed, LWC is the liquid water content,h c r is a gas constant, C is a convective heat transfer coefficient p To fix the pressure specific heat capacity, L e Is the number of Lewis, p sat Is saturated vapor pressure, T s And T ∞ Respectively the surface of the liquid film and the infinite incoming flow temperature,
Figure BDA0003264620850000141
is a wall fixing unit C S Boundary normal vector, h, of boundary surface i i Is a wall fixing unit C S Liquid film thickness of boundary surface i>
Figure BDA0003264620850000142
Is a wall fixing unit C S Liquid film flow velocity of boundary surface i, L i Is a wall fixing unit C S Boundary length of boundary surface i;
the conservation of momentum equation is:
Figure BDA0003264620850000143
Figure BDA0003264620850000144
Figure BDA0003264620850000145
Figure BDA0003264620850000146
Figure BDA0003264620850000147
Figure BDA0003264620850000148
wherein ζ is a shape generated by curve coordinate transformationFactor, p w Is the pressure in the liquid film, and is the pressure in the liquid film,
Figure BDA0003264620850000149
gravitational acceleration vector, < >>
Figure BDA00032646208500001410
And->
Figure BDA00032646208500001411
Shear forces acting on the free and bottom surfaces of the liquid film, respectively,/-, are applied to the liquid film>
Figure BDA00032646208500001412
For vector cross-multiplying symbols, p a Is the ambient pressure, σ is the surface tension coefficient, +.>
Figure BDA0003264620850000151
Is that of a no operator.
The energy conservation equation is:
Figure BDA0003264620850000152
wherein the method comprises the steps of
Figure BDA0003264620850000153
Kinetic energy and internal energy for impinging droplets, +.>
Figure BDA0003264620850000154
Is the convective heat transfer between the liquid film and the airflow field,
Figure BDA0003264620850000155
Energy carried away by liquid film evaporation, +.>
Figure BDA0003264620850000156
Latent heat generated for freezing liquid film, +.>
Figure BDA0003264620850000157
And->
Figure BDA0003264620850000158
The energy input item and the energy output item related to the liquid film transportation are specifically shown as follows:
Figure BDA0003264620850000159
Figure BDA00032646208500001510
subscripts in and out respectively represent liquid film inflow and outflow boundaries of the solid wall unit;
the icing thickness h obtained by calculation is calculated ice And updating to the current solid wall unit to form a new solid wall surface.
When S is<S critical Calculating the icing thickness by adopting a discrete liquid film icing phase calculation method; and updating the fixed wall surface;
the method for calculating the icing phase of the discrete liquid film is that the icing quality of the liquid film is solved by the following equation simultaneously
Figure BDA00032646208500001511
Further obtaining the icing thickness h of the fusion liquid film unit ice ；
Figure BDA0003264620850000161
Figure BDA0003264620850000162
F＝F contact +F between +F impingement
Figure BDA0003264620850000163
Wherein F is contact The solid-liquid-gas three-phase line acting force of discrete liquid film F between For discrete liquid film interaction force, F impingement Generating forces for water droplets striking a discrete liquid film, cp w Is the constant pressure specific heat capacity of water, T is the temperature of a discrete liquid film,
Figure BDA0003264620850000164
is the mutual energy transport between the discrete liquid films.
1-5 let t=t+1, if t<t end Returning to the step 1-2; if t=t end And (5) ending the calculation.
Therefore, the calculation method of the invention determines whether the fused liquid film unit is a discrete liquid film or a continuous liquid film by tracking the movement track of water drops and comparing the fused liquid film unit after the liquid film fusion with a set value, and then calculates the icing thickness according to a corresponding discrete liquid film icing phase calculation method or a continuous liquid film icing phase calculation method. The invention breaks the traditional continuous liquid film assumption, and independently calculates the conditions of the discrete liquid film, so that the simulated ice shape is more accurate, and the accurate simulation of the complex ice shape can be realized.
The foregoing description of the preferred embodiments of the invention is not intended to be limiting, but rather is intended to cover all modifications, equivalents, and alternatives falling within the spirit and principles of the invention.

Claims (1)

1. The icing phase change calculation method based on coupling of the continuous liquid film and the discrete liquid film is characterized by comprising the following steps of:
1-1.t =0, t being time;
1-2, modeling a flow field, constructing a water drop motion equation, tracking a water drop track, and calculating a liquid water collection coefficient beta;
1-3, tracking the movement process of a discrete liquid film unit formed after liquid drops strike a wall surface, and calculating the area S of each fused liquid film unit by fusing the liquid film units into a larger circular liquid film unit when different liquid film units meet;
1-4, traversing each fused liquid film unit and judging the area of each fused liquid film unitS and set point S critical Is a relationship of (2);
when S is>S critical Calculating the icing thickness by adopting a continuous liquid film icing phase calculation method; and updating the fixed wall surface;
when S is<S critical Calculating the icing thickness by adopting a discrete liquid film icing phase calculation method; and updating the fixed wall surface;
1-5.t =t+1, if t<t end Returning to the step 1-2; if t>t end Ending the calculation;
in the step 1-3, the calculation method of the fused liquid film unit area S is as follows:
Figure FDA0004130887000000011
wherein m is add To fuse the mass sum of the liquid film, h ave To average thickness of the fused liquid film ρ w Is the density of the liquid film;
in the step 1-4, the continuous liquid film icing phase calculation method is to solve the mass conservation equation, the momentum conservation equation and the energy conservation equation of the simultaneous continuous liquid film icing model, and solve the liquid film icing mass
Figure FDA0004130887000000012
Further obtaining the icing thickness h of the fusion liquid film unit ice ，/>
Figure FDA0004130887000000013
ρ ice Is ice density;
the mass conservation equation is:
Figure FDA0004130887000000021
Figure FDA0004130887000000022
h is a wall fixing unit C S Is used for the liquid film thickness of the liquid crystal display panel,
Figure FDA0004130887000000023
for the average speed in the liquid film, +.>
Figure FDA0004130887000000024
Is a wall fixing unit C S Boundary normal vector of S Cs Is C S L is C S Boundary line length, < >>
Figure FDA0004130887000000025
For the mass of the droplets striking the wall-fixing unit, +.>
Figure FDA0004130887000000026
For evaporating mass and liquid film
Figure FDA0004130887000000027
For the icing quality of the liquid film, the subscript i takes a solid wall unit C S η is the normal coordinate in the boundary layer of the liquid film;
Figure FDA0004130887000000028
Figure FDA0004130887000000029
Figure FDA00041308870000000210
wherein U is ∞ LWC is the liquid water content, h, for the free incoming wind speed c R is a gas constant, C is a convective heat transfer coefficient p To fix the pressure specific heat capacity, L e Is the number of Lewis, p sat Is saturated vapor pressure, T s And T ∞ Respectively the surface of the liquid film and the infinite incoming flow temperature,
Figure FDA00041308870000000211
is a wall fixing unit C S Boundary normal vector, h, of boundary surface i i Is a wall fixing unit C S Liquid film thickness of boundary surface i>
Figure FDA00041308870000000212
Is a wall fixing unit C S Liquid film flow velocity of boundary surface i, L i Is a wall fixing unit C S Boundary length of boundary surface i;
the conservation of momentum equation is:
Figure FDA00041308870000000213
Figure FDA00041308870000000214
Figure FDA00041308870000000215
wherein ζ is a shape factor, p, generated by curve coordinate transformation w Is the pressure in the liquid film, and is the pressure in the liquid film,
Figure FDA00041308870000000216
the force of gravity is applied to the acceleration vector,
Figure FDA0004130887000000031
and->
Figure FDA0004130887000000032
Shear forces acting on the free and bottom surfaces of the liquid film, respectively,/-, are applied to the liquid film>
Figure FDA00041308870000000317
For vector cross-multiplying symbols, p a For ambient pressure, σ is the surface tension coefficient, +.v is that of the Don't operator;
the energy conservation equation is:
Figure FDA0004130887000000033
wherein the method comprises the steps of
Figure FDA0004130887000000034
Kinetic energy and internal energy for impinging droplets, +.>
Figure FDA0004130887000000035
Is convective heat transfer between the liquid film and the air flow field, < >>
Figure FDA0004130887000000036
Energy carried away by liquid film evaporation, +.>
Figure FDA0004130887000000037
Latent heat generated for freezing liquid film, +.>
Figure FDA0004130887000000038
And->
Figure FDA0004130887000000039
The energy input item and the energy output item related to the liquid film transportation are specifically shown as follows:
Figure FDA00041308870000000310
Figure FDA00041308870000000311
subscripts in and out respectively represent liquid film inflow and outflow boundaries of the solid wall unit;
in the step 1-4, the method for calculating the icing phase of the discrete liquid film is that the icing quality of the liquid film is solved by combining the following equations
Figure FDA00041308870000000312
Further obtaining the icing thickness h of the fusion liquid film unit ice ；
Figure FDA00041308870000000313
Figure FDA00041308870000000314
Figure FDA00041308870000000315
Wherein F is contact The solid-liquid-gas three-phase line acting force of discrete liquid film F between For discrete liquid film interaction force, F impingement Generating forces for water droplets striking a discrete liquid film, cp w Is the constant pressure specific heat capacity of water, T is the temperature of a discrete liquid film,
Figure FDA00041308870000000316
the mutual energy transportation between the discrete liquid films is realized;
in the step 1-2, the water drop motion equation is as follows:
Figure FDA0004130887000000041
wherein m is d Is the mass of liquid drop,
Figure FDA0004130887000000042
Is the position of the liquid dropVector, t is time, ρ w 、ρ a Respectively the liquid film and the air density, V d For the volume of the drop>
Figure FDA0004130887000000043
Gravitational acceleration vector, < >>
Figure FDA0004130887000000044
For the air flow speed>
Figure FDA0004130887000000045
For droplet velocity, A d Is the cross-sectional area of the liquid drop, C d For the droplet drag coefficient, the following empirical solution is followed:
Figure FDA0004130887000000046
Figure FDA0004130887000000047
wherein Re is the Reynolds number of the liquid drop, re d Is the relative reynolds number of the droplet;
the liquid water collection coefficient beta is calculated as follows:
Figure FDA0004130887000000048
wherein the method comprises the steps of
Figure FDA0004130887000000049
Is infinitely distant mass flow flux, < >>
Figure FDA00041308870000000410
Is impacted on the fixed wall unit C in unit time S The liquid water mass flux contained by the upper droplet particle package,
Figure FDA00041308870000000411
Figure FDA00041308870000000412
wherein the method comprises the steps of
Figure FDA00041308870000000413
For striking against the fixed wall unit C S The mass of the liquid drops in the wall fixing unit, N is the total number of the liquid drops which strike the wall fixing unit in dt time, j is the number of the liquid drops, S Cs Is a wall fixing unit C S Area of U ∞ Is the free incoming flow velocity. />
CN202111084586.0A 2021-06-09 2021-09-15 Icing phase change calculation method based on coupling of continuous liquid film and discrete liquid film Active CN113779904B (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2021106450242 2021-06-09
CN202110645024 2021-06-09

Publications (2)

Publication Number Publication Date
CN113779904A CN113779904A (en) 2021-12-10
CN113779904B true CN113779904B (en) 2023-04-25

Family

ID=78844410

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111084586.0A Active CN113779904B (en) 2021-06-09 2021-09-15 Icing phase change calculation method based on coupling of continuous liquid film and discrete liquid film

Country Status (1)

Country Link
CN (1) CN113779904B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114169077B (en) * 2021-12-13 2023-04-07 南京航空航天大学 Strong-coupling three-dimensional numerical simulation method for hot gas anti-icing of aircraft engine inlet part
CN114180072B (en) * 2022-02-16 2022-04-12 中国空气动力研究与发展中心低速空气动力研究所 Icing thickness detection method
CN117408053B (en) * 2023-10-18 2024-05-07 中国空气动力研究与发展中心计算空气动力研究所 Method for establishing low-temperature flat plate drying mode frosting characteristic curve under strong convection condition

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104298886A (en) * 2014-10-20 2015-01-21 上海电机学院 Icing 3-D numerical simulation method of aeroengine rotating part
CN112345196A (en) * 2021-01-11 2021-02-09 中国空气动力研究与发展中心低速空气动力研究所 Super-cooled large water drop splash simulation test device and method

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682145B (en) * 2011-11-30 2013-11-27 天津空中代码工程应用软件开发有限公司 Numerical simulation method of flight icing
US8696151B1 (en) * 2013-03-12 2014-04-15 Tcp Reliable, Inc. Monitoring shipment of biological products to determine remaining refrigerant quantity
CN105277485B (en) * 2015-09-24 2017-08-25 空气动力学国家重点实验室 It is a kind of to determine ice and the device of object plane cohesive force
CN108460217B (en) * 2018-03-13 2021-10-01 西北工业大学 Unsteady three-dimensional icing numerical simulation method
CN109376403B (en) * 2018-09-29 2023-06-23 南京航空航天大学 Two-dimensional icing numerical simulation method based on Cartesian self-adaptive reconstruction technology
CN111291505B (en) * 2020-05-08 2020-10-09 中国空气动力研究与发展中心低速空气动力研究所 Wing-type icing shape prediction method and device based on depth confidence network

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104298886A (en) * 2014-10-20 2015-01-21 上海电机学院 Icing 3-D numerical simulation method of aeroengine rotating part
CN112345196A (en) * 2021-01-11 2021-02-09 中国空气动力研究与发展中心低速空气动力研究所 Super-cooled large water drop splash simulation test device and method

Also Published As

Publication number Publication date
CN113779904A (en) 2021-12-10

Similar Documents

Publication Publication Date Title
CN113779904B (en) Icing phase change calculation method based on coupling of continuous liquid film and discrete liquid film
Villedieu et al. Glaciated and mixed phase ice accretion modeling using ONERA 2D icing suite
CN104268399B (en) The computational methods of model parameter in icing wind tunnel test under the conditions of big water droplet are subcooled
Wright et al. Semi-empirical modelling of SLD physics
CN111396269B (en) Multi-time-step unsteady icing calculation method and system and storage medium
Cao et al. Numerical simulation of supercooled large droplet icing phenomenon: a review
CN104354867B (en) A kind of method for designing of supercool big water droplet icing detector
CN114896906B (en) Ice accretion simulation method considering heat conduction in ice layer and solid wall surface
CN114139393B (en) Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer
CN114169077B (en) Strong-coupling three-dimensional numerical simulation method for hot gas anti-icing of aircraft engine inlet part
CN112678189B (en) Improved icing sensor installation position determining method
CN111967201A (en) Method for analyzing critical icing type based on numerical simulation model
Cao et al. Numerical simulation of rime ice accretions on an aerofoil using an Eulerian method
CN114398844B (en) Steady-state anti-icing simulation method based on continuous water film flow
Wang et al. Lagrangian approach for simulating supercooled large droplets’ impingement effect
Kollár et al. Modeling the evolution of droplet size distribution in two-phase flows
CN107798198A (en) Physical-based melting phenomenon realistic simulation method
Vahab et al. An adaptive coupled level set and moment-of-fluid method for simulating droplet impact and solidification on solid surfaces with application to aircraft icing
Hann UAV Icing: Challenges for computational fluid dynamic (CFD) tools
Shen et al. Effects of upstream component and air injection on water droplet impingement characteristics for downstream surfaces
Liu et al. Icing performance of stratospheric airship in ascending process
GB2572352A (en) Methods and apparatus for simulating liquid collection on aerodynamic components
Hann Numerical Simulation of In-Flight Icing of Unmanned Aerial Vehicles
CN106682289A (en) Method and device for predicting ice coating disaster development trend of electric transmission line
Szilder et al. Three-dimensional numerical simulation of ice accretion using a discrete morphogenetic approach

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant