CN113779904A - 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
CN113779904A
CN113779904A CN202111084586.0A CN202111084586A CN113779904A CN 113779904 A CN113779904 A CN 113779904A CN 202111084586 A CN202111084586 A CN 202111084586A CN 113779904 A CN113779904 A CN 113779904A
Authority
CN
China
Prior art keywords
liquid film
liquid
icing
discrete
unit
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.)
Granted
Application number
CN202111084586.0A
Other languages
Chinese (zh)
Other versions
CN113779904B (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

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 the steps of modeling a flow field and constructing a water drop motion equation; tracking the track of the liquid drop, fusing different liquid film units into a larger circular liquid film unit when the different liquid film units meet, and calculating the area of each fused liquid film unit; traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unitcriticalThe relationship of (1); when S is>ScriticalCalculating the icing thickness by adopting a continuous liquid film icing phase calculation method; when S is<ScriticalCalculating the icing thickness by adopting a discrete liquid film icing phase calculation method, and updating the solid wall profileAnd finally obtaining the simulated ice shape with a complex structure. Compared with the traditional method for calculating the icing phase change by only adopting the continuous liquid film hypothesis, the icing phase change calculation method based on the coupling of the continuous liquid film and the discrete liquid film enables the simulated ice shape to be more accurate and can realize the accurate simulation of the complex ice shape.

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 factor that threatens the flight safety of an aircraft. The phenomenon of icing of an aviation aircraft is very easy to occur when the aviation aircraft passes through a cloud layer containing supercooled water drops. The ice accretion attached to the surfaces of core components such as wings can not only reduce the lift force of the airplane and increase the resistance, but also cause important potential safety hazards such as loss of pitching and course stability of the airplane. Currently, the icing research means of the aviation aircraft mainly comprise real flight tests, sprinkler flight tests, ground icing wind tunnel tests, icing numerical calculation and the like, wherein the icing numerical calculation is an important means which cannot be obtained by the icing safety assessment of the current aviation aircraft due to the advantages of rapidness, high efficiency, economy and the like.
At present, an icing numerical value 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 transition calculation is to realize the simulation of the growth of the ice accretion on the surface of the aircraft through the calculation of the heat and mass transfer process of a 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 mainstream mainly comprises the following steps: the Messinger model, the Myers model, the ShallowWater model, and the like. The Messinger model simplifies the icing process of the surface of the airplane into a simple mass and energy conservation equation, and the solution is realized through the calculation of an algebraic equation of a fixed wall surface unit; the Myers model further considers the energy transport in the normal direction of a liquid film and an ice layer on the basis of the Messinger model by introducing a Stefen condition, so that the calculation precision is improved; from the angle of numerical calculation, the Shallow Water model develops a more comprehensive liquid film transport 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 aviation aircraft.
Despite the rapid development of current icing phase change calculation models, accurate simulation of highly complex ice shapes, such as three-dimensional conch ice typical of swept-back wings, still cannot be achieved.
Disclosure of Invention
In order to solve the technical problem that accurate simulation of highly complex ice shapes 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.
In the long-term research process, the existing icing phase change model is based on the assumption of a continuous liquid film, so that the behavior of the liquid film under the characteristic condition of a highly complex discrete ice shape cannot be accurately characterized, and the actual complex ice shape has both the continuous liquid film and the discrete liquid film. Therefore, the invention breaks the traditional continuous liquid film hypothesis, calculates the liquid water collection characteristic based on Monte Carlo Method (Monte Carlo Method), and converts the sub-packets of discrete supercooled water droplets into liquid film units which can discretely move on the solid wall surface. On the basis of tracking the track of the discrete liquid film unit, calculating the heat and mass transfer process in 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 to be 0, wherein t is time;
1-2, modeling a flow field, constructing a water drop motion equation, tracking a liquid drop track, and calculating a liquid water collection coefficient beta;
1-3, tracking the motion process of a discrete liquid film unit formed after liquid drops impact a wall surface, fusing different liquid film units into a larger circular liquid film unit when the different liquid film units meet, and calculating the area S of each fused liquid film unit;
1-4, traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unitcriticalThe relationship of (1);
when S is>ScriticalCalculating the icing thickness by adopting a continuous liquid film icing phase calculation method; and updating the shape of the solid wall;
when S is<ScriticalCalculating the icing thickness by adopting a discrete liquid film icing phase calculation method; and updating the shape of the solid wall;
1-5. let t be t +1, if t<tendReturning to the step 1-2; if t is tendAnd ending the calculation.
Further, in the step 1-3, the calculation method of the unit area S of the fused liquid film is as follows:
Figure BDA0003264620850000031
wherein m isaddTo fuse the mass sum of the liquid films, haveIs the average thickness of the liquid film after fusion, pwIs the liquid film density.
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 icing mass of the liquid film is solved
Figure BDA0003264620850000041
Further obtaining the icing thickness of the fused liquid film unithice，
Figure BDA0003264620850000042
ρiceIs the ice density.
Further, the mass conservation equation is:
Figure BDA0003264620850000043
Figure BDA0003264620850000044
h is a wall fixing unit CSThe thickness of the liquid film of (a),
Figure BDA0003264620850000045
is the average velocity within the liquid film,
Figure BDA0003264620850000046
is a wall-fixing unit CSBoundary normal vector of (S)CsIs CSL is CSThe length of the boundary line of (a),
Figure BDA0003264620850000047
the mass of the liquid drops impacting the solid wall unit,
Figure BDA0003264620850000048
For the liquid film evaporating mass and
Figure BDA0003264620850000049
the subscript i is taken as a fixed wall unit C for the icing quality of the liquid filmSEta is the normal coordinate in the boundary layer of the liquid film;
Figure BDA00032646208500000410
Figure BDA00032646208500000411
Figure BDA0003264620850000051
wherein U is∞For free incoming flow velocity, LWC is liquid water content, hcR is a gas constant, C for the convective heat transfer coefficientpSpecific heat capacity at constant pressure, LeIs the number of Liu Yi Si, psatTo saturated vapor pressure, TsAnd T∞Respectively the liquid film surface and the infinite incoming flow temperature,
Figure BDA0003264620850000052
is a wall-fixing unit CSBoundary normal vector of boundary surface i, hi is fixed wall unit CSThe thickness of the liquid film at the boundary surface i,
Figure BDA0003264620850000053
is a wall-fixing unit CSLiquid film flow velocity, L, of boundary surface iiIs a wall-fixing unit CSThe boundary length of boundary surface i;
the conservation of momentum equation is:
Figure BDA0003264620850000054
Figure BDA0003264620850000055
Figure BDA0003264620850000056
Figure BDA0003264620850000057
Figure BDA0003264620850000058
Figure BDA0003264620850000059
where ζ is the shape factor due to the curve coordinate transformation, pwIs the pressure within the liquid film and,
Figure BDA0003264620850000061
is a vector of the acceleration of gravity,
Figure BDA0003264620850000062
and
Figure BDA0003264620850000063
respectively acting on the free surface and the bottom surface of the liquid film,
Figure BDA0003264620850000064
for vector cross-multiplication of symbols, paIs the ambient pressure, σ is the surface tension coefficient,
Figure BDA0003264620850000065
is that operator; A. b, C is an intermediate variable.
The energy conservation equation is as follows:
Figure BDA0003264620850000066
wherein
Figure BDA0003264620850000067
The kinetic energy and the internal energy brought by the impact of the liquid drops,
Figure BDA0003264620850000068
Is the convection heat exchange between a liquid film and an airflow field,
Figure BDA0003264620850000069
The energy taken away for the evaporation of the liquid film,
Figure BDA00032646208500000610
The latent heat generated by the freezing of the liquid film,
Figure BDA00032646208500000611
and
Figure BDA00032646208500000612
the energy input item and the energy output item related to liquid film transportation are specifically as follows:
Figure BDA00032646208500000613
Figure BDA00032646208500000614
the subscripts in and out characterize the liquid film inflow and outflow boundaries of the solid-wall cell, respectively.
1. The method for calculating the icing phase change based on the coupling of the continuous liquid film and the discrete liquid film as claimed in claim 4, wherein in the steps 1 to 6, the method for calculating the icing phase of the discrete liquid film comprises the following equations to solve the icing quality of the liquid film
Figure BDA0003264620850000071
Further obtaining the icing thickness h of the fused liquid film unitice；
Figure BDA0003264620850000072
Figure BDA0003264620850000073
F＝Fcontact+Fbetween+Fimpingement
Figure BDA0003264620850000074
Wherein FcontactSolid-liquid-gas three-phase line force as a discrete liquid film, FbetweenFor discrete liquid film interaction force, FimpingementThe force generated by the impact of the water droplets on the discrete liquid film, F being the resultant force, CpwIs 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 step 1-2, the water drop motion equation is:
Figure BDA0003264620850000076
wherein m isdIs the mass of the liquid drop,
Figure BDA0003264620850000077
Is the droplet position vector, t is time, pw、ρaRespectively liquid film and air density, VdIs the volume of the liquid drop and is,
Figure BDA0003264620850000078
is a vector of the acceleration of gravity,
Figure BDA0003264620850000081
is the speed of the air flow and is,
Figure BDA0003264620850000082
as the droplet velocity, AdIs the cross-sectional area of the droplet, CdThe drop resistance coefficient is solved according to the following empirical formula:
Figure BDA0003264620850000083
Figure BDA0003264620850000084
where Re is the Reynolds number of the droplet, RedIs the relative reynolds number of the droplet.
Further, the liquid water collection coefficient β is calculated as:
Figure BDA0003264620850000085
wherein
Figure BDA0003264620850000086
The flux of mass flow at infinite distance,
Figure BDA0003264620850000087
Is impacted on the fixed wall unit C in unit timeSThe liquid water mass flux contained in the droplet particles above,
Figure BDA0003264620850000088
Figure BDA0003264620850000089
wherein
Figure BDA0003264620850000091
For impacting on the wall-fixing unit CSThe mass of the liquid drops in the chamber, N is the total number of liquid drops impacting on the solid wall unit in dt time, j is the number of liquid drops, SCsIs a wall-fixing unit CSArea of, U∞Is the free incoming flow velocity.
Compared with the traditional method for calculating the icing phase change based on the coupling of the continuous liquid film and the discrete liquid film, the method for calculating the icing phase change based on the coupling of the continuous liquid film and the discrete liquid film 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 discrete shell ice of the sweepback wing.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings needed to be used in the embodiments of the present invention or in the description of the prior art will be briefly described below, and it is obvious that the drawings described below are only some embodiments of the present invention, and it is obvious for those skilled in the art that other drawings can be obtained according to the drawings without creative efforts.
Fig. 1 is a flowchart 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 particular examples set forth below are illustrative only and are not intended to be limiting.
An icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film comprises the following steps, as shown in fig. 1:
1-1, setting a time starting point t, and enabling t to be 0;
1-2, modeling a flow field, constructing a water drop motion equation, tracking a liquid drop track, and calculating a liquid water collection coefficient beta;
the present invention introduces the concept of a packet of supercooled water droplets, i.e. it is believed that the typical characteristic supercooled water droplet dynamics in a packet of particles can be representative of the nature of the water droplets within the packet of particles. In order to obtain the motion trail of the particle package, a water drop motion equation is constructed as follows:
Figure BDA0003264620850000101
wherein m isdIs the mass of the liquid drop,
Figure BDA0003264620850000102
Being liquid dropletsPosition vector, t time, pw、ρaRespectively liquid film and air density, VdIs the volume of the liquid drop and is,
Figure BDA0003264620850000103
is a vector of the acceleration of gravity,
Figure BDA0003264620850000104
is the speed of the air flow and is,
Figure BDA0003264620850000105
as the droplet velocity, AdIs the cross-sectional area of the droplet, CdThe drop resistance coefficient is solved according to the following empirical formula:
Figure BDA0003264620850000106
Figure BDA0003264620850000107
where Re is the Reynolds number of the droplet, RedIs the relative reynolds number of the droplet; on the basis, the quality of liquid water brought by different particle packets collected in unit time in the solid wall unit is counted to obtain the liquid water collection characteristic, and the calculation of further obtaining the liquid water collection coefficient beta is as follows:
Figure BDA0003264620850000111
wherein the content of the first and second substances,
Figure BDA0003264620850000112
the flux of mass flow at infinite distance,
Figure BDA0003264620850000113
Is impacted on the fixed wall unit C in unit timeSThe liquid water mass flux contained in the droplet particles above,
Figure BDA0003264620850000114
Figure BDA0003264620850000115
wherein
Figure BDA0003264620850000116
For impacting on the wall-fixing unit CSThe mass of the liquid drops in the chamber, N is the total number of liquid drops impacting on the solid wall unit in dt time, j is the number of liquid drops, SCsIs a wall-fixing unit CSArea of, U∞Is the free incoming flow velocity.
1-3, tracking the motion process of a discrete liquid film unit formed after liquid drops impact a wall surface by using a water drop motion equation, and assuming that the mass of a particle packet impacting on a fixed wall unit is so small that a circular unit attached to the fixed wall can be formed after the particle packet impacts on the fixed wall, so that the liquid drops impact on the fixed wall to form a discrete liquid film, and when different discrete liquid films are fused, the discrete liquid film is converted into a continuous liquid film.
When different discrete liquid films meet, the discrete liquid films are fused into a larger circular liquid film unit, the center position of the fused liquid film unit is located at the centroid average position of different liquid film units, the liquid film quality is the sum of the quality of each liquid film before meeting, and the area S of each fused liquid film unit is calculated:
Figure BDA0003264620850000121
wherein m isaddTo fuse the mass sum of the liquid films, haveIs the average thickness of the liquid film after fusion, pwIs the liquid film density.
1-4, traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unitcriticalThe relationship of (1);
when S is>ScriticalFreezing with a continuous liquid filmCalculating the icing thickness by a phase calculation method; and updating the shape of the solid wall;
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 icing mass of the liquid film is solved
Figure BDA0003264620850000122
Further obtaining the icing thickness h of the fused liquid film unitice，
Figure BDA0003264620850000123
ρiceIs the ice density.
Wherein the mass conservation equation is as follows:
Figure BDA0003264620850000131
Figure BDA0003264620850000132
h is a wall fixing unit CSThe thickness of the liquid film of (a),
Figure BDA0003264620850000133
is the average velocity within the liquid film,
Figure BDA0003264620850000134
is a wall-fixing unit CSBoundary normal vector of (S)CsIs CSL is CSThe length of the boundary line of (a),
Figure BDA0003264620850000135
the mass of the liquid drops impacting the solid wall unit,
Figure BDA0003264620850000136
For the liquid film evaporating mass and
Figure BDA0003264620850000137
for freezing the liquid filmMass, subscript i, taken as wall-fixing unit CSEta is the normal coordinate in the boundary layer of the liquid film;
Figure BDA0003264620850000138
Figure BDA0003264620850000139
Figure BDA00032646208500001310
wherein U is∞For free incoming flow velocity, LWC is liquid water content, hcR is a gas constant, C for the convective heat transfer coefficientpSpecific heat capacity at constant pressure, LeIs the number of Liu Yi Si, psatTo saturated vapor pressure, TsAnd T∞Respectively the liquid film surface and the infinite incoming flow temperature,
Figure BDA0003264620850000141
is a wall-fixing unit CSBoundary normal vector, h, of boundary surface iiIs a wall-fixing unit CSThe thickness of the liquid film at the boundary surface i,
Figure BDA0003264620850000142
is a wall-fixing unit CSLiquid film flow velocity, L, of boundary surface iiIs a wall-fixing unit CSThe boundary length of boundary surface i;
the conservation of momentum equation is:
Figure BDA0003264620850000143
Figure BDA0003264620850000144
Figure BDA0003264620850000145
Figure BDA0003264620850000146
Figure BDA0003264620850000147
Figure BDA0003264620850000148
where ζ is the shape factor due to the curve coordinate transformation, pwIs the pressure within the liquid film and,
Figure BDA0003264620850000149
is a vector of the acceleration of gravity,
Figure BDA00032646208500001410
and
Figure BDA00032646208500001411
respectively acting on the free surface and the bottom surface of the liquid film,
Figure BDA00032646208500001412
for vector cross-multiplication of symbols, paIs the ambient pressure, σ is the surface tension coefficient,
Figure BDA0003264620850000151
is that not pull operator.
The energy conservation equation is:
Figure BDA0003264620850000152
wherein
Figure BDA0003264620850000153
The kinetic energy and the internal energy brought by the impact of the liquid drops,
Figure BDA0003264620850000154
Is the convection heat exchange between a liquid film and an airflow field,
Figure BDA0003264620850000155
The energy taken away for the evaporation of the liquid film,
Figure BDA0003264620850000156
The latent heat generated by the freezing of the liquid film,
Figure BDA0003264620850000157
and
Figure BDA0003264620850000158
the energy input item and the energy output item related to liquid film transportation are specifically as follows:
Figure BDA0003264620850000159
Figure BDA00032646208500001510
subscripts in and out represent liquid film inflow and outflow boundaries of the solid-wall unit respectively;
the calculated icing thickness hiceAnd updating to the current solid wall unit to form a new solid wall surface.
When S is<ScriticalCalculating the icing thickness by adopting a discrete liquid film icing phase calculation method; and updating the shape of the solid wall;
the calculation method of the icing phase of the discrete liquid film comprises the following equations simultaneously to solve the icing quality of the liquid film
Figure BDA00032646208500001511
Further obtaining the icing thickness h of the fused liquid film unitice；
Figure BDA0003264620850000161
Figure BDA0003264620850000162
F＝Fcontact+Fbetween+Fimpingement
Figure BDA0003264620850000163
Wherein FcontactSolid-liquid-gas three-phase line force as a discrete liquid film, FbetweenFor discrete liquid film interaction force, FimpingementGenerating a force for the water droplets to impact the discrete liquid film, CpwIs 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 be t +1, if t<tendReturning to the step 1-2; if t is tendAnd 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 the water drops and comparing the fused liquid film unit with the set value after the liquid films are fused, and then calculates the icing thickness according to the corresponding discrete liquid film icing phase calculation method or the continuous liquid film icing phase calculation method. The method breaks through the traditional continuous liquid film assumption, and performs independent calculation on the condition 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 above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.

Claims (7)

1. An icing phase change calculation method based on coupling of a continuous liquid film and a discrete liquid film is characterized by comprising the following steps of:
1-1, setting a time starting point t of calculation, and enabling t to be 0;
1-2, modeling a flow field, constructing a water drop motion equation, tracking a liquid drop track, and calculating a liquid water collection coefficient beta;
1-3, tracking the motion process of a discrete liquid film unit formed after liquid drops impact a wall surface, fusing different liquid film units into a larger circular liquid film unit when the different liquid film units meet, and calculating the area S of each fused liquid film unit;
1-4, traversing each fused liquid film unit, and judging the area S and the set value S of each fused liquid film unitcriticalThe relationship of (1);
when S is>ScriticalCalculating the icing thickness by adopting a continuous liquid film icing phase calculation method; and updating the shape of the solid wall;
when S is<ScriticalCalculating the icing thickness by adopting a discrete liquid film icing phase calculation method; and updating the shape of the solid wall;
1-5. let t be t +1, if t<tendReturning to the step 1-2; if t is tendAnd ending the calculation.
2. The method for calculating the phase change of ice formation based on the coupling of the continuous liquid film and the discrete liquid film as claimed in claim 1, wherein in the step 1-3, the calculation method of the unit area S of the fused liquid film is as follows:
Figure FDA0003264620840000011
wherein m isaddTo fuse the mass sum of the liquid films, haveIs the average thickness of the liquid film after fusion, pwIs the liquid film density.
3. The device of claim 2, based on a continuous liquid film and a discrete liquidThe method for calculating the icing phase change of the membrane coupling is characterized in that in the step 1-4, the method for calculating the icing phase of the continuous liquid membrane is a mass conservation equation, a momentum conservation equation and an energy conservation equation of a simultaneous continuous liquid membrane icing model, and the icing mass of the liquid membrane is solved
Figure FDA0003264620840000021
Further obtaining the icing thickness h of the fused liquid film unitice，
Figure FDA0003264620840000022
ρiceIs the ice density.
4. The icing phase transition calculation method based on the coupling of the continuous liquid film and the discrete liquid film as claimed in claim 3, wherein the mass conservation equation is as follows:
Figure FDA0003264620840000023
Figure FDA0003264620840000024
h is a wall fixing unit CSThe thickness of the liquid film of (a),
Figure FDA0003264620840000025
is the average velocity within the liquid film,
Figure FDA0003264620840000026
is a wall-fixing unit CSBoundary normal vector of (S)CsIs CSL is CSThe length of the boundary line of (a),
Figure FDA0003264620840000027
the mass of the liquid drops impacting the solid wall unit,
Figure FDA0003264620840000028
For the liquid film evaporating mass and
Figure FDA0003264620840000029
the subscript i is taken as a fixed wall unit C for the icing quality of the liquid filmSEta is the normal coordinate in the boundary layer of the liquid film;
Figure FDA0003264620840000031
Figure FDA0003264620840000032
Figure FDA0003264620840000033
wherein U is∞For free incoming flow velocity, LWC is liquid water content, hcR is a gas constant, C for the convective heat transfer coefficientpSpecific heat capacity at constant pressure, LeIs the number of Liu Yi Si, psatTo saturated vapor pressure, TsAnd T∞Respectively the liquid film surface and the infinite incoming flow temperature,
Figure FDA0003264620840000034
is a wall-fixing unit CSBoundary normal vector, h, of boundary surface iiIs a wall-fixing unit CSThe thickness of the liquid film at the boundary surface i,
Figure FDA0003264620840000035
is a wall-fixing unit CSLiquid film flow velocity, L, of boundary surface iiIs a wall-fixing unit CSThe boundary length of boundary surface i;
the conservation of momentum equation is:
Figure FDA0003264620840000041
Figure FDA0003264620840000042
Figure FDA0003264620840000043
Figure FDA0003264620840000044
Figure FDA0003264620840000045
Figure FDA0003264620840000046
where ζ is the shape factor due to the curve coordinate transformation, pwIs the pressure within the liquid film and,
Figure FDA0003264620840000047
is a vector of the acceleration of gravity,
Figure FDA0003264620840000048
and
Figure FDA0003264620840000049
respectively acting on the free surface and the bottom surface of the liquid film,
Figure FDA00032646208400000410
for vector cross-multiplication of symbols, paIs the ambient pressure, σ is the surface tension coefficient,
Figure FDA00032646208400000411
is that operator; A. b, C is an intermediate variable;
the energy conservation equation is as follows:
Figure FDA00032646208400000412
wherein
Figure FDA00032646208400000413
The kinetic energy and the internal energy brought by the impact of the liquid drops,
Figure FDA00032646208400000414
Is the convection heat exchange between a liquid film and an airflow field,
Figure FDA00032646208400000415
The energy taken away for the evaporation of the liquid film,
Figure FDA00032646208400000416
The latent heat generated by the freezing of the liquid film,
Figure FDA0003264620840000051
and
Figure FDA0003264620840000052
the energy input item and the energy output item related to liquid film transportation are specifically as follows:
Figure FDA0003264620840000053
Figure FDA0003264620840000054
the subscripts in and out characterize the liquid film inflow and outflow boundaries of the solid-wall cell, respectively.
5. The method for calculating the icing phase change based on the coupling of the continuous liquid film and the discrete liquid film as claimed in claim 4, wherein in the steps 1 to 4, the method for calculating the icing phase of the discrete liquid film comprises the following equations to solve the icing quality of the liquid film
Figure FDA0003264620840000055
Further obtaining the icing thickness h of the fused liquid film unitice；
Figure FDA0003264620840000056
Figure FDA0003264620840000057
F＝Fcontact+Fbetween+Fimpingement
Figure FDA0003264620840000058
Wherein FcontactSolid-liquid-gas three-phase line force as a discrete liquid film, FbetweenFor discrete liquid film interaction force, FimpingementThe force generated by the impact of the water droplets on the discrete liquid film, F being the resultant force, CpwIs the constant pressure specific heat capacity of water, T is the temperature of a discrete liquid film,
Figure FDA0003264620840000061
is the mutual energy transport between the discrete liquid films.
6. The icing phase transition calculation method based on the coupling of the continuous liquid film and the discrete liquid film as claimed in claim 1, wherein in the step 1-2, the water drop motion equation is as follows:
Figure FDA0003264620840000062
wherein m isdIs the mass of the liquid drop,
Figure FDA0003264620840000063
Is the droplet position vector, t is time, pw、ρaRespectively liquid film and air density, VdIs the volume of the liquid drop and is,
Figure FDA0003264620840000064
is a vector of the acceleration of gravity,
Figure FDA0003264620840000065
is the speed of the air flow and is,
Figure FDA0003264620840000066
as the droplet velocity, AdIs the cross-sectional area of the droplet, CdThe drop resistance coefficient is solved according to the following empirical formula:
Figure FDA0003264620840000067
Figure FDA0003264620840000068
where Re is the Reynolds number of the droplet, RedIs the relative reynolds number of the droplet.
7. The method according to claim 6, wherein the liquid water collection coefficient β is calculated as:
Figure FDA0003264620840000071
wherein
Figure FDA0003264620840000072
The flux of mass flow at infinite distance,
Figure FDA0003264620840000073
Is impacted on the fixed wall unit C in unit timeSThe liquid water mass flux contained in the droplet particles above,
Figure FDA0003264620840000074
Figure FDA0003264620840000075
wherein
Figure FDA0003264620840000076
For impacting on the wall-fixing unit CSThe mass of the liquid drops in the chamber, N is the total number of liquid drops impacting on the solid wall unit in dt time, j is the number of liquid drops, SCsIs a wall-fixing unit CSArea 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
CN202110645024 2021-06-09
CN2021106450242 2021-06-09

Publications (2)

Publication Number Publication Date
CN113779904A true CN113779904A (en) 2021-12-10
CN113779904B 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)

Cited By (3)

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

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682145A (en) * 2011-11-30 2012-09-19 天津空中代码工程应用软件开发有限公司 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
CN104298886A (en) * 2014-10-20 2015-01-21 上海电机学院 Icing 3-D numerical simulation method of aeroengine rotating part
CN105277485A (en) * 2015-09-24 2016-01-27 空气动力学国家重点实验室 Ice and object surface adhesion force measuring device
CN108460217A (en) * 2018-03-13 2018-08-28 西北工业大学 A kind of unstable state three-dimensional icing method for numerical simulation
CN109376403A (en) * 2018-09-29 2019-02-22 南京航空航天大学 A kind of two-dimentional icing method for numerical simulation based on Descartes's self-adapting reconstruction technology
CN111291505A (en) * 2020-05-08 2020-06-16 中国空气动力研究与发展中心低速空气动力研究所 Wing-type icing shape prediction method and device based on depth confidence network
CN112345196A (en) * 2021-01-11 2021-02-09 中国空气动力研究与发展中心低速空气动力研究所 Super-cooled large water drop splash simulation test device and method

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682145A (en) * 2011-11-30 2012-09-19 天津空中代码工程应用软件开发有限公司 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
CN104298886A (en) * 2014-10-20 2015-01-21 上海电机学院 Icing 3-D numerical simulation method of aeroengine rotating part
CN105277485A (en) * 2015-09-24 2016-01-27 空气动力学国家重点实验室 Ice and object surface adhesion force measuring device
CN108460217A (en) * 2018-03-13 2018-08-28 西北工业大学 A kind of unstable state three-dimensional icing method for numerical simulation
CN109376403A (en) * 2018-09-29 2019-02-22 南京航空航天大学 A kind of two-dimentional icing method for numerical simulation based on Descartes's self-adapting reconstruction technology
CN111291505A (en) * 2020-05-08 2020-06-16 中国空气动力研究与发展中心低速空气动力研究所 Wing-type icing shape prediction method and device based on depth confidence network
CN112345196A (en) * 2021-01-11 2021-02-09 中国空气动力研究与发展中心低速空气动力研究所 Super-cooled large water drop splash simulation test device and method

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JOSE M. MOLINAR-MONTERRUBIO等: "Internal combustion engine parametric identification scheme for misfire fault detection: Experimental results" *
任靖豪等: "复杂构型水滴收集率的拉格朗日计算方法" *
李胜超: "三维笛卡尔网格的结冰数值模拟" *

Cited By (5)

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

Also Published As

Publication number Publication date
CN113779904B (en) 2023-04-25

Similar Documents

Publication Publication Date Title
CN113779904A (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
CN111396269B (en) Multi-time-step unsteady icing calculation method and system and storage medium
CN114896906B (en) Ice accretion simulation method considering heat conduction in ice layer and solid wall surface
CN114169077B (en) Strong-coupling three-dimensional numerical simulation method for hot gas anti-icing of aircraft engine inlet part
CN114139393B (en) Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer
CN112678189B (en) Improved icing sensor installation position determining method
CN104354867A (en) Design method of big supercooling water droplet icing detector and detector
Cao et al. Numerical simulation of rime ice accretions on an aerofoil using an Eulerian method
Ayan et al. In-flight ice accretion simulation in mixed-phase conditions
CN114398844B (en) Steady-state anti-icing simulation method based on continuous water film flow
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
Liu et al. Icing performance of stratospheric airship in ascending process
Chang et al. Three-dimensional modelling and simulation of the ice accretion process on aircraft wings
Grenestedt et al. Dynamic soaring in hurricanes
WO2019186151A1 (en) Methods and apparatus for simulating liquid collection on aerodynamic components
Müller et al. UAV Icing: 3D Simulations of Propeller Icing Effects and Anti-Icing Heat Loads
Bendarkar et al. Rapid Assessment of Power Requirements and Optimization of Thermal Ice Protection Systems
Serkan et al. Parallel computing applied to three-dimensional droplet trajectory simulation in Lagrangian approach
Hann Numerical Simulation of In-Flight Icing of Unmanned Aerial Vehicles
Dong et al. Calculation and analysis of water film flow characteristics on anti-icing airfoil surface
Grift et al. Computational method for ice crystal trajectories in a turbofan compressor
Özgen et al. Ice accretion simulations on airfoils
Bai et al. Analysis of the influence of different droplet size distribution on wing impact characteristics and ice shape
Khalil et al. Flight Simulation and Drag Prediction for a Pitching-Accelerating Hypersonic Reentry Vehicle

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