CN114139393A - Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer - Google Patents

Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer Download PDF

Info

Publication number
CN114139393A
CN114139393A CN202111478862.1A CN202111478862A CN114139393A CN 114139393 A CN114139393 A CN 114139393A CN 202111478862 A CN202111478862 A CN 202111478862A CN 114139393 A CN114139393 A CN 114139393A
Authority
CN
China
Prior art keywords
component
water
water film
icing
flow
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
CN202111478862.1A
Other languages
Chinese (zh)
Other versions
CN114139393B (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.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
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 Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202111478862.1A priority Critical patent/CN114139393B/en
Publication of CN114139393A publication Critical patent/CN114139393A/en
Application granted granted Critical
Publication of CN114139393B publication Critical patent/CN114139393B/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
    • 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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • 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/06—Power analysis or power 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/08—Thermal analysis or thermal optimisation
    • 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
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Fluid Mechanics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

The invention discloses a three-dimensional anti-icing numerical simulation method for electric heating of a component considering water film flow heat transfer, which comprises the following steps of: carrying out three-dimensional anti-icing two-phase flow modeling on the part; meshing the anti-icing two-phase flow calculation domain of the part by using commercial software; solving a Reynolds average N-S equation to obtain an air flow field; solving a supercooled water drop motion equation and calculating a local water collection coefficient of the surface of the component; solving a continuous equation, a momentum equation and an energy equation of the water film flow of the anti-icing surface of the component, and calculating the water film flow, heat transfer and phase change of the surface of the component; the method adopts a numerical simulation method, accurately and efficiently simulates the three-dimensional anti-icing process of the component, and lays an important foundation for the engineering analysis of three-dimensional electric heating anti-icing.

Description

Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer
Technical Field
The invention relates to the technical field of electric heating anti-icing, in particular to a component electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer.
Background
In the aeronautical field, icing is one of the most serious risks that may be encountered. Icing of the wings or engines may cause loss of maneuverability and controllability of the aircraft, seriously threatens flight safety, and the resulting flight accidents often occur. The research reports provided by the NASA scientific and technical research project in the united states show that the percentage of icing-related events in a total accident of up to 16.2% is the highest among the flight accidents with casualties.
The electric heating anti-icing method is the most widely used anti-icing method for modern airplanes, and has high efficiency, controllability and safety. For the research of electric heating anti-icing of the component, the research is mostly developed aiming at a two-dimensional model in China, the research on three-dimensional water film flowing on the surface of the component in the electric heating process is less, and a momentum equation of the water film flowing established by three-dimensional anti-icing software FENSAP-ICE developed abroad is only a simple integral equation considering the influence of airflow drag force in a show water model. At home and abroad, a calculation method for simulating an electric heating anti-icing process by taking an energy item obtained by calculating the flowing heat transfer and phase change of a water film as a heat source item of an anti-icing surface is not available.
Disclosure of Invention
Based on the problems, the invention provides a component electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer, which realizes accurate and efficient simulation of the component anti-icing process; and an important foundation is laid for engineering analysis of three-dimensional electric heating anti-icing.
In order to solve the above problems, the present invention provides a method for simulating an electric heating three-dimensional anti-icing numerical value of a component considering water film flow heat transfer, comprising the steps of:
s1: and establishing a calculation domain of the anti-icing two-phase flow calculation of the part according to the geometrical parameters of the part.
S2: the computational domain is meshed using business software.
S3: acquiring flow field working conditions (the speed, the temperature, the liquid water content, the diameter of the supercooled water drops, the turbulence intensity and the like of incoming air and supercooled water drops) and electric heating power, and calculating an air flow field through a Reynolds average Navier-Stokes equation (RANS) to obtain a flow field solution, wherein the calculation formula is as follows:
Figure BDA0003394680840000021
Figure BDA0003394680840000022
Figure BDA0003394680840000023
Figure BDA0003394680840000024
where ρ isaIs the density of air; t is time, uiIs the velocity component; i and j respectively take 1, 2 and 3 to represent three coordinates of x, y and z; p is the pressure on the fluid microcell; mu.saIs the kinetic viscosity of air; k is the turbulence energy; mu.stIs the turbulent viscosity; e represents the total energy; k is a radical ofeffIs the effective thermal conductivity; t is the temperature; (τ)ij)effIs the bias stress tensor; shIndicating an internal heat source.
Calculating the heat conduction in the solid of the component to obtain the temperature of the electric heating surface of the component under the heat exchange of air flow, wherein the solid heat conduction calculation formula is as follows:
Figure BDA0003394680840000025
wherein h represents the enthalpy of the solid; lambda [ alpha ]sRepresents the thermal conductivity of the solid;
Figure BDA0003394680840000026
is the heat transferred due to the rotational or translational motion of the solid.
S4: based on the calculation result of S3, the movement speed of the supercooled water drops is calculated according to the flow field working conditions and the geometric shapes of the components
Figure BDA0003394680840000027
The calculation formula is as follows:
Figure BDA0003394680840000028
Figure BDA0003394680840000029
Figure BDA00033946808400000210
Figure BDA00033946808400000211
where ρ iswIs the density of the supercooled water droplets; alpha is alphawIs the volume fraction of supercooled water droplets; u. ofa1,ua2,ua3Is the velocity of air
Figure BDA0003394680840000031
Components at three coordinates; u. ofd1,ud2,ud3As the velocity of water droplets
Figure BDA0003394680840000032
Components at three coordinates; g1，g2，g3Is acceleration of gravity
Figure BDA0003394680840000033
Components at three coordinates;
Figure BDA0003394680840000034
is the momentum exchange coefficient between the air and the supercooled water droplets; d is the diameter of the water droplet;
Figure BDA0003394680840000035
is the drag coefficient;
Figure BDA0003394680840000036
is the drag force coefficient;
Figure BDA0003394680840000037
is the relative reynolds number.
The local water collection coefficient beta of the surface of the part can be calculated according to the calculation result of the supercooled water drop track and is as follows:
Figure BDA0003394680840000038
wherein u isd,normalIs the normal velocity of the water droplets impacting the wall surface; u. ofd,∞Is the incoming flow velocity of the supercooled water droplets; LWC is the liquid water content.
S5: according to the calculation result of S4, calculating the water film flow, heat transfer and phase change of the surface of the component, and obtaining the water film thickness H of the surface of the component by solving the continuous equation, momentum equation and energy equation of the water film flow of the anti-icing surface of the componentwFlow velocity of
Figure BDA0003394680840000039
Height of ice formation HiEnergy of super-cooled water drops impacting on component surface
Figure BDA00033946808400000310
Energy taken away by surface water film evaporation
Figure BDA00033946808400000311
And heat released by surface water film icing
Figure BDA00033946808400000312
The continuous equation calculation formula of the current water film flow is as follows:
Figure BDA00033946808400000313
Figure BDA00033946808400000314
Figure BDA00033946808400000315
wherein n represents a time step;
Figure BDA00033946808400000316
the sum of the mass of the water films flowing out of each unit surface of the component surface infinitesimal control body is represented; u shape∞Is the incoming flow velocity; h represents the surface convection heat transfer coefficient; rhoaIs the density of air; c. Cp,aRepresents the specific heat capacity at constant pressure of air; sc is the Schmidt number; mwIs the molecular weight of water; r represents a molar gas constant; p is a radical ofwall，TwallRespectively controlling the saturated vapor pressure and temperature of the surface of the body; p is a radical of∞，T∞Respectively saturated vapour pressure and temperature outside the boundary layer.
The momentum equation calculation formula of the water film flow is as follows:
Figure BDA0003394680840000041
wherein
Figure BDA0003394680840000042
Is the pressure gradient of the water film;
Figure BDA0003394680840000043
is the water film flow velocity vector; the positive direction of the z coordinate axis is the outward normal direction vertical to the wall surface;
Figure BDA0003394680840000044
boundary conditions of the upper and lower surfaces of the water film;
Figure BDA0003394680840000045
is the component surface shear force vector.
According to the momentum equation of water film flow and the boundary conditions of upper and lower surfaces of water film, the speed of water film flow
Figure BDA0003394680840000046
The calculation formula of (2) is as follows:
Figure BDA0003394680840000047
is as follows;
according to the energy equation of water film flow
Figure BDA0003394680840000048
The calculation formula of (2) is as follows:
Figure BDA0003394680840000049
Figure BDA00033946808400000410
Figure BDA00033946808400000411
Figure BDA00033946808400000412
wherein
Figure BDA00033946808400000413
Heat provided for electrical heating; l isfLatent heat of freezing; l iseIs the steaming of waterGenerating latent heat; t iswIs the temperature of the water film; t isdIs the temperature of the supercooled water droplets impinging on the surface of the component; u. ofdIs the velocity of the supercooled water droplets; t isaIs the temperature of the incoming air.
Height H of ice formation on the surface of the componentiThe calculation formula of (2) is as follows:
Figure BDA0003394680840000051
wherein
Figure BDA0003394680840000052
The original ice thickness on the surface of the part, initially,
Figure BDA0003394680840000053
is 0.
Thickness H of water film on surface of componentwThe calculation formula of (2) is as follows:
Figure BDA0003394680840000054
wherein
Figure BDA0003394680840000055
The original thickness of the water film on the surface of the part is obtained, initially,
Figure BDA0003394680840000056
is 0.
S6: according to the calculation result of S5, the energy items obtained by the calculation of S5 are processed
Figure BDA0003394680840000057
And as a heat source item of the anti-icing surface of the component, performing coupling calculation of an external air flow field, internal electric heating and solid heat conduction of the component to obtain the temperature distribution of the anti-icing surface of the component.
Further, when the energy change caused by the impact of supercooled water drops and the evaporation and phase change of water films on the surface of the component is not considered, there is no change
Figure BDA0003394680840000058
And
Figure BDA0003394680840000059
the three items, the coupling calculation of the external air, the internal electric heating and the solid heat conduction of the component is directly carried out by commercial software, and the temperature of the surface of the component under the air flow heat exchange in the step S3 is obtained. Energy term obtained by calculating flow heat transfer and phase change of water film
Figure BDA00033946808400000510
As a heat source item of the anti-icing surface of the component, the coupling calculation of the external air, the internal electric heating and the solid heat conduction of the component is performed again, and the temperature distribution of the anti-icing surface of the component in step S6, which takes into account the impact of supercooled water drops and the evaporation and phase change of the surface water film, is obtained.
Compared with the prior art, the beneficial results of the invention are as follows: the complicated electric heating anti-icing process is divided into four modules to be solved: calculating an external air flow field of the component, calculating flow and impact characteristics of supercooled water drops, calculating heat transfer and phase change of water film flow of an anti-icing surface, and calculating the temperature of the anti-icing surface; the momentum equation of water film flow is not only a simple integral equation considering the influence of airflow drag force, but is established based on an N-S equation to realize three-dimensional flow simulation of the water film on the surface of the part; the method can accurately and efficiently simulate the three-dimensional anti-icing process of the component, and lays an important foundation for the engineering analysis of the three-dimensional electric heating anti-icing of the component.
Drawings
FIG. 1 is a flow chart of a method for simulating an electric heating three-dimensional anti-icing numerical value of a component considering water film flow heat transfer in the embodiment.
Fig. 2 is a schematic view of the profile of a NACA0012 airfoil.
FIG. 3 is a schematic illustration of seven electric heaters distributed on the leading edge of an NACA0012 airfoil.
FIG. 4 is a graph comparing calculated and experimentally measured results for the temperature of the anti-icing surface of an airfoil of NACA 0012.
Detailed Description
In order to make the objects, calculation steps and advantages of the present invention easily understood by those skilled in the art, the present invention will be further described in detail with reference to the following examples and the accompanying drawings, which are used for illustrating the present invention and are not to be construed as limiting the present invention.
Example (b):
the invention takes a three-dimensional component as a research object, simulates the electric heating anti-icing process of the component by considering the flowing, heat transfer and phase change of a water film on the surface of the component, and comprises the following specific implementation steps:
s1: and establishing a calculation domain of the anti-icing two-phase flow calculation of the part according to the geometrical parameters of the part.
S2: the computational domain is meshed using business software.
S3: acquiring flow field working conditions (the speed, the temperature, the liquid water content, the diameter of the supercooled water drops, the turbulence intensity and the like of incoming air and supercooled water drops) and electric heating power, and calculating an air flow field through a Reynolds average Navier-Stokes equation (RANS) to obtain a flow field solution, wherein the calculation formula is as follows:
Figure BDA0003394680840000061
Figure BDA0003394680840000062
Figure BDA0003394680840000063
Figure BDA0003394680840000064
where ρ isaIs the density of air; t is time, uiIs the velocity component; i and j respectively take 1, 2 and 3 to represent three coordinates of x, y and z; p is the pressure on the fluid microcell; mu.saIs the kinetic viscosity of air; k is the turbulence energy; mu.stIs the turbulent viscosity; e represents the total energy; k is a radical ofeffIs the effective thermal conductivity; t is the temperature; (τ)ij)effIs the bias stress tensor; shIndicating an internal heat source.
Calculating the heat conduction in the solid of the component to obtain the temperature of the electric heating surface of the component under the heat exchange of air flow, wherein the solid heat conduction calculation formula is as follows:
Figure BDA0003394680840000071
wherein h represents the enthalpy of the solid; lambda [ alpha ]sRepresents the thermal conductivity of the solid;
Figure BDA0003394680840000072
is the heat transferred due to the rotational or translational motion of the solid.
S4: based on the calculation result of S3, the movement speed of the supercooled water drops is calculated according to the flow field working conditions and the geometric shapes of the components
Figure BDA0003394680840000073
The calculation formula is as follows:
Figure BDA0003394680840000074
Figure BDA0003394680840000075
Figure BDA0003394680840000076
Figure BDA0003394680840000077
where ρ iswIs the density of the supercooled water droplets; alpha is alphawIs the volume fraction of supercooled water droplets; u. ofa1,ua2,ua3Is the velocity of air
Figure BDA0003394680840000078
Components at three coordinates; u. ofd1,ud2,ud3As the velocity of water droplets
Figure BDA0003394680840000079
Components at three coordinates; g1，g2，g3Is acceleration of gravity
Figure BDA00033946808400000710
Components at three coordinates;
Figure BDA00033946808400000711
is the momentum exchange coefficient between the air and the supercooled water droplets; d is the diameter of the water droplet;
Figure BDA00033946808400000712
is the drag coefficient;
Figure BDA00033946808400000713
is the drag force coefficient;
Figure BDA00033946808400000714
is the relative reynolds number.
The local water collection coefficient beta of the surface of the part can be calculated according to the calculation result of the supercooled water drop track
Figure BDA00033946808400000715
Wherein u isd,normalIs the normal velocity of the water droplets impacting the wall surface; u. ofd,∞Is the incoming flow velocity of the supercooled water droplets; LWC is the liquid water content.
S5: according to the calculation result of S4, calculating the water film flow, heat transfer and phase change of the surface of the component, and obtaining the water film thickness H of the surface of the component by solving the continuous equation, momentum equation and energy equation of the water film flow of the anti-icing surface of the componentwFlow velocity of
Figure BDA0003394680840000081
Height of ice formation HiEnergy of super-cooled water drops impacting on component surface
Figure BDA0003394680840000082
Energy taken away by surface water film evaporation
Figure BDA0003394680840000083
And heat released by surface water film icing
Figure BDA0003394680840000084
The continuous equation calculation formula of the current water film flow is as follows:
Figure BDA0003394680840000085
Figure BDA0003394680840000086
Figure BDA0003394680840000087
wherein n represents a time step;
Figure BDA0003394680840000088
the sum of the mass of the water films flowing out of each unit surface of the component surface infinitesimal control body is represented; u shape∞Is the incoming flow velocity; h represents the surface convection heat transfer coefficient; rhoaIs the density of air; c. Cp,aRepresents the specific heat capacity at constant pressure of air; sc is the Schmidt number; mwIs the molecular weight of water; r represents a molar gas constant; p is a radical ofwall，TwallRespectively controlling the saturated vapor pressure and temperature of the surface of the body; p is a radical of∞，T∞Respectively saturated vapour pressure and temperature outside the boundary layer.
The momentum equation calculation formula of the water film flow is as follows:
Figure BDA0003394680840000089
wherein
Figure BDA00033946808400000810
Is the pressure gradient of the water film;
Figure BDA00033946808400000811
is the water film flow velocity vector; the positive direction of the z coordinate axis is the outward normal direction vertical to the wall surface;
Figure BDA00033946808400000812
boundary conditions of the upper and lower surfaces of the water film;
Figure BDA00033946808400000813
is the component surface shear force vector.
According to the momentum equation of water film flow and the boundary conditions of upper and lower surfaces of water film, the speed of water film flow
Figure BDA00033946808400000814
The calculation formula of (2) is as follows:
Figure BDA0003394680840000091
according to the energy equation of water film flow
Figure BDA0003394680840000092
The calculation formula of (2) is as follows:
Figure BDA0003394680840000093
Figure BDA0003394680840000094
Figure BDA0003394680840000095
Figure BDA0003394680840000096
wherein
Figure BDA0003394680840000097
Heat provided for electrical heating; l isfLatent heat of freezing; l iseIs the latent heat of evaporation of water; t iswIs the temperature of the water film; t isdIs the temperature of the supercooled water droplets impinging on the surface of the component; u. ofdIs the velocity of the supercooled water droplets; t isaIs the temperature of the incoming air.
Height H of ice formation on the surface of the componentiThe calculation formula of (2) is as follows:
Figure BDA0003394680840000098
wherein
Figure BDA0003394680840000099
The original ice thickness on the surface of the part, initially,
Figure BDA00033946808400000910
is 0.
Thickness H of water film on surface of componentwThe calculation formula of (2) is as follows:
Figure BDA00033946808400000911
wherein
Figure BDA00033946808400000912
The original thickness of the water film on the surface of the part is obtained, initially,
Figure BDA00033946808400000913
is 0.
S6: according to the calculation result of S5, the energy items obtained by the calculation of S5 are processed
Figure BDA00033946808400000914
And as a heat source item of the anti-icing surface of the component, performing coupling calculation of an external air flow field, internal electric heating and solid heat conduction of the component to obtain the temperature distribution of the anti-icing surface of the component.
Further, when the energy change caused by the impact of supercooled water drops and the evaporation and phase change of water films on the surface of the component is not considered, there is no change
Figure BDA00033946808400000915
And
Figure BDA00033946808400000916
the three items, the coupling calculation of the external air, the internal electric heating and the solid heat conduction of the component is directly carried out by commercial software, and the temperature of the surface of the component under the air flow heat exchange in the step S3 is obtained. Energy term obtained by calculating flow heat transfer and phase change of water film
Figure BDA0003394680840000101
As a heat source item of the anti-icing surface of the component, the coupling calculation of the external air, the internal electric heating and the solid heat conduction of the component is performed again, and the temperature distribution of the anti-icing surface of the component in step S6, which takes into account the impact of supercooled water drops and the evaporation and phase change of the surface water film, is obtained.
The invention will be verified and its feasibility analysed by experimental data as follows:
taking an NACA0012 airfoil (shown in figure 1) with the chord length of 0.9144m as an example, seven independent electric heaters (shown in figure 2) are distributed on the leading edge of the airfoil, and the heating power of the heater A is 4.805kW/m2The heating power of the heater B is 4.03kW/m2The heating power of the heater C is 2.945kW/m2The heating power of the heater D is 4.805kW/m2The heating power of the heater E is 3.41kW/m2The heating power of the heater F is 2.635kW/m2The heating power of the heater G is 2.325kW/m2The velocity of the incoming air and supercooled water droplets was 44.7m/s, the temperature was-6.67 ℃ and the liquid water content was 0.78g/m3The method has no attack angle, and is used for carrying out numerical simulation analysis on the electric heating three-dimensional anti-icing process of the wing and simultaneously carrying out comparison verification on the temperature of the anti-icing surface of the wing and experimental data, and the method is as follows.
Establishing a calculation domain of wing anti-icing two-phase flow calculation according to geometrical parameters of wings, meshing the calculation domain by adopting commercial software, firstly carrying out coupling calculation of wing external air flow heat transfer, internal electric heating and solid heat conduction to obtain the heat flow distribution of an air flow field outside the wings and an anti-icing surface, then calculating a motion equation of super-cooled water drops to obtain the motion speed of the super-cooled water drops, calculating the local water collection coefficient distribution of the wing surface, then calculating the flow, heat transfer and phase change of a wing anti-icing surface water film to obtain energy terms brought by water drop impact and water film heat transfer, parameters such as the thickness and flow speed of the water film, finally adding the energy terms to the wing anti-icing surface, coupling external air heat transfer, internal electric heating and solid heat conduction calculation again to obtain the anti-icing surface temperature distribution.
And finally, verifying the calculation result of the wing anti-icing surface temperature by comparing the experiment measurement result, wherein FIG. 4 is a comparison graph of the calculation result of the NACA0012 wing section anti-icing surface temperature and the experiment measurement result. As shown in fig. 4, the temperature error between the calculated result and the experimental data is within 3 ℃, which proves the feasibility of the invention.
The above is an embodiment of the present invention. The embodiments and specific parameters in the embodiments are only for the purpose of clearly illustrating the verification process of the invention and are not intended to limit the scope of the invention, which is defined by the claims, and all equivalent structural changes made by using the contents of the specification and the drawings of the present invention should be covered by the scope of the present invention.

Claims (2)

1. The method for simulating the three-dimensional anti-icing numerical value of the electric heating component by considering the flowing heat transfer of the water film is characterized by comprising the following steps of:
s1: establishing a calculation domain of the anti-icing two-phase flow calculation of the component according to the geometric parameters of the component;
s2: carrying out grid division on the calculation domain by adopting commercial software;
s3: acquiring flow field working conditions (the speed, the temperature, the liquid water content, the diameter of the supercooled water drops, the turbulence intensity and the like of incoming air and supercooled water drops) and electric heating power, and calculating an air flow field through a Reynolds average Navier-Stokes equation (RANS) to obtain a flow field solution, wherein the calculation formula is as follows:
Figure FDA0003394680830000011
Figure FDA0003394680830000012
Figure FDA0003394680830000013
Figure FDA0003394680830000014
where ρ isaIs the density of air; t is time, uiIs the velocity component; i and j respectively take 1, 2 and 3 to represent three coordinates of x, y and z; p is the pressure on the fluid microcell; mu.saIs the kinetic viscosity of air; k is the turbulence energy; mu.stIs the turbulent viscosity; e represents the total energy; k is a radical ofeffIs the effective thermal conductivity; t is the temperature; (τ)ij)effIs the bias stress tensor; shRepresents an internal heat source;
calculating the heat conduction in the solid of the component to obtain the temperature of the electric heating surface of the component under the heat exchange of air flow, wherein the solid heat conduction calculation formula is as follows:
Figure FDA0003394680830000015
wherein h represents the enthalpy of the solid; lambda [ alpha ]sRepresents the thermal conductivity of the solid;
Figure FDA0003394680830000016
heat transferred due to the rotational or translational movement of the solid;
s4: based on the calculation result of S3, the movement speed of the supercooled water drops is calculated according to the flow field working conditions and the geometric shapes of the components
Figure FDA0003394680830000017
The calculation formula is as follows:
Figure FDA0003394680830000021
Figure FDA0003394680830000022
Figure FDA0003394680830000023
Figure FDA0003394680830000024
where ρ iswIs the density of the supercooled water droplets; alpha is alphawIs the volume fraction of supercooled water droplets; u. ofa1,ua2,ua3Is the velocity of air
Figure FDA0003394680830000025
Components at three coordinates; u. ofd1,ud2,ud3As the velocity of water droplets
Figure FDA0003394680830000026
Components at three coordinates; g1，g2，g3Is acceleration of gravity
Figure FDA0003394680830000027
Components at three coordinates;
Figure FDA0003394680830000028
is the momentum exchange coefficient between the air and the supercooled water droplets; d is the diameter of the water droplet;
Figure FDA0003394680830000029
is the drag coefficient;
Figure FDA00033946808300000210
is the drag force coefficient;
Figure FDA00033946808300000211
is the relative Reynolds number;
the local water collection coefficient beta of the surface of the part can be calculated according to the calculation result of the supercooled water drop track
Figure FDA00033946808300000212
Wherein u isd,normalIs the normal velocity of the water droplets impacting the wall surface; u. ofd,∞Is the incoming flow velocity of the supercooled water droplets; LWC is liquid water content;
s5: according to the calculation result of S4, calculating the water film flow, heat transfer and phase change of the surface of the component, and obtaining the water film thickness H of the surface of the component by solving the continuous equation, momentum equation and energy equation of the water film flow of the anti-icing surface of the componentwFlow velocity of
Figure FDA00033946808300000213
Height of ice formation HiImpingement of super-cooled water droplets on component surfacesBrought energy
Figure FDA00033946808300000214
Energy taken away by surface water film evaporation
Figure FDA00033946808300000215
And heat released by surface water film icing
Figure FDA00033946808300000216
The continuous equation calculation formula of the current water film flow is as follows:
Figure FDA0003394680830000031
Figure FDA0003394680830000032
Figure FDA0003394680830000033
wherein n represents a time step;
Figure FDA0003394680830000034
the sum of the mass of the water films flowing out of each unit surface of the component surface infinitesimal control body is represented; u shape∞Is the incoming flow velocity; h represents the surface convection heat transfer coefficient; rhoaIs the density of air; c. Cp,aRepresents the specific heat capacity at constant pressure of air; sc is the Schmidt number; mwIs the molecular weight of water; r represents a molar gas constant; p is a radical ofwall，TwallRespectively controlling the saturated vapor pressure and temperature of the surface of the body; p is a radical of∞，T∞Respectively the saturated steam pressure and the temperature outside the boundary layer;
the momentum equation calculation formula of the water film flow is as follows:
Figure FDA0003394680830000035
wherein
Figure FDA0003394680830000036
Is the pressure gradient of the water film;
Figure FDA0003394680830000037
is the water film flow velocity vector; the positive direction of the z coordinate axis is the outward normal direction vertical to the wall surface;
Figure FDA0003394680830000038
boundary conditions of the upper and lower surfaces of the water film;
Figure FDA0003394680830000039
is the component surface shear force vector;
according to the momentum equation of water film flow and the boundary conditions of upper and lower surfaces of water film, the speed of water film flow
Figure FDA00033946808300000310
The calculation formula of (2) is as follows:
Figure FDA00033946808300000311
according to the energy equation of water film flow
Figure FDA00033946808300000312
The calculation formula of (2) is as follows:
Figure FDA00033946808300000313
Figure FDA00033946808300000314
Figure FDA0003394680830000041
Figure FDA0003394680830000042
wherein
Figure FDA0003394680830000043
Heat provided for electrical heating; l isfLatent heat of freezing; l iseIs the latent heat of evaporation of water; t iswIs the temperature of the water film; t isdIs the temperature of the supercooled water droplets impinging on the surface of the component; u. ofdIs the velocity of the supercooled water droplets; t isaIs the temperature of the incoming air;
height H of ice formation on the surface of the componentiThe calculation formula of (2) is as follows:
Figure FDA0003394680830000044
wherein
Figure FDA0003394680830000045
The original ice thickness on the surface of the part, initially,
Figure FDA0003394680830000046
is 0;
thickness H of water film on surface of componentwThe calculation formula of (2) is as follows:
Figure FDA0003394680830000047
wherein
Figure FDA0003394680830000048
The original thickness of the water film on the surface of the part is obtained, initially,
Figure FDA0003394680830000049
is 0;
s6: according to the calculation result of S5, the energy items obtained by the calculation of S5 are processed
Figure FDA00033946808300000410
And as a heat source item of the anti-icing surface of the component, performing coupling calculation of an external air flow field, internal electric heating and solid heat conduction of the component to obtain the temperature distribution of the anti-icing surface of the component.
2. The electric heating three-dimensional anti-icing numerical simulation method of the component considering the water film flow heat transfer according to claim 1, is characterized in that: when the energy change caused by the impact of supercooled water drops and the evaporation and phase change of water films on the surfaces of components is not considered, the energy change does not exist
Figure FDA00033946808300000411
And
Figure FDA00033946808300000412
the coupling calculation of external air, internal electric heating and solid heat conduction of the component is directly carried out through commercial software, and the temperature of the surface of the component under the condition of air flow heat exchange is obtained; then calculating the energy term obtained by the flow heat transfer and phase change of the water film
Figure FDA00033946808300000413
And as a heat source item of the anti-icing surface of the component, performing coupling calculation of external air, internal electric heating and solid heat conduction of the component again to obtain the temperature distribution of the anti-icing surface of the component, which takes the impact of supercooled water drops and surface water film evaporation and phase change into consideration.
CN202111478862.1A 2021-12-06 2021-12-06 Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer Active CN114139393B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111478862.1A CN114139393B (en) 2021-12-06 2021-12-06 Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111478862.1A CN114139393B (en) 2021-12-06 2021-12-06 Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer

Publications (2)

Publication Number Publication Date
CN114139393A true CN114139393A (en) 2022-03-04
CN114139393B CN114139393B (en) 2023-04-07

Family

ID=80384398

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111478862.1A Active CN114139393B (en) 2021-12-06 2021-12-06 Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer

Country Status (1)

Country Link
CN (1) CN114139393B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115587506A (en) * 2022-12-09 2023-01-10 四川大学 Design method of electric heating ice preventing and removing system
CN115688453A (en) * 2022-11-08 2023-02-03 东莞理工学院 Performance calculation method of cross-flow membrane module for flue gas heat and humidity recovery
CN115817822A (en) * 2023-02-09 2023-03-21 中国空气动力研究与发展中心低速空气动力研究所 Heat load distribution design method of electric heating anti-icing system

Citations (7)

* 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
WO2013078628A1 (en) * 2011-11-30 2013-06-06 天津空中代码工程应用软件开发有限公司 Flight icing numerical simulation method of helicopter rotor wing
CN103218507A (en) * 2012-12-13 2013-07-24 中国电力科学研究院 Two-dimensional numerical simulation method for icing process of power transmission line
CN104298886A (en) * 2014-10-20 2015-01-21 上海电机学院 Icing 3-D numerical simulation method of aeroengine rotating part
JP2015215772A (en) * 2014-05-12 2015-12-03 株式会社東芝 Heat transfer simulation device and heat transfer simulation method
CN109751204A (en) * 2019-02-18 2019-05-14 中国空气动力研究与发展中心低速空气动力研究所 A kind of wind energy conversion system icing method for numerical simulation
CN110759458A (en) * 2019-09-19 2020-02-07 江苏新宜中澳环境技术有限公司 Optimization design method for reaction distribution in ozone reactor

Patent Citations (7)

* 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
WO2013078628A1 (en) * 2011-11-30 2013-06-06 天津空中代码工程应用软件开发有限公司 Flight icing numerical simulation method of helicopter rotor wing
CN103218507A (en) * 2012-12-13 2013-07-24 中国电力科学研究院 Two-dimensional numerical simulation method for icing process of power transmission line
JP2015215772A (en) * 2014-05-12 2015-12-03 株式会社東芝 Heat transfer simulation device and heat transfer simulation method
CN104298886A (en) * 2014-10-20 2015-01-21 上海电机学院 Icing 3-D numerical simulation method of aeroengine rotating part
CN109751204A (en) * 2019-02-18 2019-05-14 中国空气动力研究与发展中心低速空气动力研究所 A kind of wind energy conversion system icing method for numerical simulation
CN110759458A (en) * 2019-09-19 2020-02-07 江苏新宜中澳环境技术有限公司 Optimization design method for reaction distribution in ozone reactor

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
申晓斌;郭琦;林贵平;穆作栋;王博;: "三维旋转部件复杂表面结冰数值模拟" *
苏长明;胡娅萍;曹广州;吉洪湖;: "考虑水膜蒸发的三维明冰积冰数值研究" *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115688453A (en) * 2022-11-08 2023-02-03 东莞理工学院 Performance calculation method of cross-flow membrane module for flue gas heat and humidity recovery
CN115587506A (en) * 2022-12-09 2023-01-10 四川大学 Design method of electric heating ice preventing and removing system
CN115587506B (en) * 2022-12-09 2023-03-10 四川大学 Design method of electric heating ice preventing and removing system
CN115817822A (en) * 2023-02-09 2023-03-21 中国空气动力研究与发展中心低速空气动力研究所 Heat load distribution design method of electric heating anti-icing system

Also Published As

Publication number Publication date
CN114139393B (en) 2023-04-07

Similar Documents

Publication Publication Date Title
CN114139393B (en) Part electric heating three-dimensional anti-icing numerical simulation method considering water film flow heat transfer
Beaugendre et al. FENSAP-ICE's three-dimensional in-flight ice accretion module: ICE3D
Zhu et al. 3D ice accretion simulation for complex configuration basing on improved messinger model
Son et al. Ice accretion on helicopter fuselage considering rotor-wake effects
CN114169077B (en) Strong-coupling three-dimensional numerical simulation method for hot gas anti-icing of aircraft engine inlet part
CN114896906B (en) Ice accretion simulation method considering heat conduction in ice layer and solid wall surface
Cao et al. Extension to the Myers model for calculation of three-dimensional glaze icing
CN114398844B (en) Steady-state anti-icing simulation method based on continuous water film flow
Cao et al. Numerical simulation of rime ice accretions on an aerofoil using an Eulerian method
Pena et al. Development of a three-dimensional icing simulation code in the NSMB flow solver
Dong et al. Heat transfer and temperature analysis of an aeroengine strut under icing conditions
Chang et al. Three-dimensional modelling and simulation of the ice accretion process on aircraft wings
Das et al. Ice shape prediction for turbofan rotating blades
Cao et al. Insight into rime ice accretion on an aircraft wing and corresponding effects on aerodynamic performance
Sang et al. Numerical simulation of icing effect and ice accretion on three-dimensional configurations
Dong et al. Numerical Simulation of Hot Air Anti-icing Charateristics of an Aero-engine Strut
Yirtici et al. Wind turbine performance losses due to the ice accretion on the turbine blades
Cao et al. Numerical simulation of ice accretion in supercooled large droplet conditions
Wang et al. Performance calculation of electrothermal anti-icing system on three-dimensional surface
Lei et al. Numerical investigation of the electrothermal de-icing process of a rotor blade
Wang et al. Anti-icing simulation in wet air of a piccolo system using FENSAP-ICE
Wada et al. Evaluation of aerodynamic performance of airfoil using the e-mps method after icing
Kinzel et al. A CFD Approach for Predicting 3D Ice Accretion on Aircraft
Zayni et al. RECENT DEVELOPMENTS OF THE MULTI-PHYSICS SOLVER CHAMPS-ICE
Liang et al. Thermal Feature Analysis of a New Hot-Air Anti-Icing Structure

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