An ice crystal melting and icing whole-process simulation calculation method in an aero-engine compressor passage

Abstract

The application belongs to the technical field of aero-engine icing simulation, and discloses a kind of ice crystal melting icing whole process simulation calculation method in aero-engine compressor passage, first, the fluid domain model of compressor is established and initial flow field is obtained;Subsequently, the movement, impact and phase change process of ice crystal are simulated by combining discrete phase model and user-defined function;Further, the ice accumulation is calculated based on the icing thermodynamic model, and the dynamic mesh technology is driven to update the grid node position of blade surface in real time, and the icing topography is reconstructed;Finally, the reconstructed ice model is calculated for flow field, and the aerodynamic performance loss caused by icing is analyzed.The application realizes the whole process dynamic coupling simulation from ice crystal transport, ice growth to aerodynamic influence evaluation, and provides accurate and effective numerical analysis means for compressor ice crystal icing risk assessment.

Description

Technical Field

[0001] This application belongs to the field of aero-engine icing simulation technology, specifically involving a simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine. Background Technology

[0002] The main factors contributing to icing vary depending on the cloud altitude. When flying below 5000 m, supercooled water in the clouds is the key factor leading to icing. However, at altitudes above 7000 m, water vapor exists primarily as solid ice crystals. Ice crystal formation is often closely linked to strong convective weather and hot, humid environments. Under the influence of updrafts, the altitude of water vapor increases dramatically, and the temperature drops sharply, leading to the formation of ice crystals in large quantities. Typically, the median mass diameter (MMD) of ice crystals in convective clouds ranges from 20 μm to 500 μm, and the ice crystal concentration can reach 9 g / m³. 3 When an aircraft flies through clouds containing ice crystals, the ice crystals are drawn into the engine by the airflow, thus inducing ice formation.

[0003] Ice crystal formation poses a serious potential hazard to compressors. When ice detaches from the stator blades or leading-edge components of the low-pressure compressor, it can be carried downstream of the engine, causing physical damage to the blades and potentially leading to compressor surge, stall, and engine shutdown. Simultaneously, icing obstructs airflow, reduces engine efficiency, and causes imbalances in internal engine components, thus placing higher demands on flight control and structural strength. Furthermore, ice crystals can interfere with sensor equipment, such as pitot tubes and angle-of-attack sensors, resulting in erroneous flight data, misleading pilot decisions, and potentially causing malfunctions in autopilot systems and other critical flight assistance systems, increasing the difficulty of manual operation. In extreme cases, ice crystals entering the engine can affect combustion efficiency in the combustion chamber, leading to uneven airflow and fuel mixing, and causing thrust loss. Therefore, the FAA first included engine icing caused by ice crystals in its airworthiness certification review in 2015. However, due to the location and characteristics of ice crystal formation, icing within the compressor cannot be detected immediately, which significantly hinders in-depth research on ice crystal formation.

[0004] Currently, research on ice crystal icing within compressor channels mainly involves two approaches: First, constructing a single-stage axial-flow compressor model and associated facilities to simulate the environmental characteristics of an engine compressor. However, this method is experimentally challenging, and the techniques for preparing ice crystals of different sizes, controlling ambient temperature, and collecting ice crystal icing data are still immature. Furthermore, building an experimental system requires significant financial and time investment. Second, establishing a single-stage axial-flow compressor model and using CFD simulation software to obtain the icing patterns of ice crystals within the compressor channels.

[0005] Current research methods for this problem mainly fall into two categories: experimental approaches and numerical simulation approaches. Experimental approaches involve building test benches that simulate the compressor environment to directly observe ice crystal behavior and icing consequences. However, this method faces technical challenges such as difficulty in accurately controlling ice crystal size and concentration, high costs of high-altitude environment simulation, and incomplete data acquisition. Furthermore, it is time-consuming and costly, making it difficult to cover all operating conditions. Numerical simulation approaches, especially coupled computational methods based on computational fluid dynamics (CFD), have become an important research tool in this field due to their flexibility and cost advantages. General-purpose commercial CFD software provides the basic framework for flow field solving and particle tracking, but its native models typically do not include complete descriptions of complex physical processes such as ice crystal phase transitions, impact adhesion, and ice growth. With the continuous upgrading of simulation technology and the significant improvement in computing resources and capabilities, the priority of this research approach is constantly increasing. However, since existing CFD simulation software cannot directly calculate ice crystal formation, users need to perform secondary development on existing open-source CFD simulation software in conjunction with basic equations to accurately simulate the movement and phase change process of ice crystals. Furthermore, existing simulation technology can only obtain ice thickness data on compressor blades and cannot update the morphology of compressor blades in real time during the icing process, let alone conduct aerodynamic analysis after compressor blades icing.

[0006] Therefore, neither experiments nor simulations can simulate the dynamic icing process of ice crystals in the compressor and its subsequent impact on aerodynamic performance in a real-time, complete, and low-cost manner. Summary of the Invention

[0007] To address the problem of difficulty in reconstructing the icing morphology of compressor blades in existing technologies, this application provides a simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine. By coupling the icing thermodynamic model with dynamic mesh technology, the blade icing morphology is reconstructed and updated in real time during the simulation process, thereby realizing the simulation from ice crystal transport, phase change, ice growth to post-icing aerodynamic performance evaluation, providing a technical basis for accurately analyzing the compressor performance loss caused by icing.

[0008] To achieve the above technical objectives, this application specifically adopts the following technical solution: One aspect of this application provides a simulation calculation method for the entire process of ice crystal melting and icing in the compressor passage of an aero-engine, comprising the following steps: S1. Based on the compressor's blade profile, casing profile, and hub profile data, establish a compressor blade geometric model containing at least one stage of rotor blades and one stage of stator blades, and generate a fluid domain model containing a rotating fluid domain and a stationary fluid domain based on the compressor blade geometric model; wherein the rotating fluid domain and the stationary fluid domain correspond to the regions where the rotor blades and stator blades are located, respectively. S2. Mesh the fluid domain model and import it into computational fluid dynamics software to perform steady-state flow field calculations and obtain the initial flow field distribution in the compressor channel; S3. Based on the initial flow field, enable the Discrete Phase Model (DPM), load the user-defined function (UDF) to simulate the melting, impact and rebound behavior of ice crystal particles, simulate the motion and phase change process of ice crystal particles after they are sucked into the compressor channel, and obtain the distribution and state of ice crystal particles on the blade surface. S4. Based on the distribution and state of ice crystal particles on the blade surface, an icing thermodynamic model is established using the loaded UDF. The surface mass, momentum and energy conservation equations are solved to calculate the icing rate and icing amount on the blade surface. Combining the dynamic mesh method, the position of the mesh nodes on the blade surface is updated according to the icing amount to reconstruct the compressor blade morphology model after icing. S5. Mesh the reconstructed compressor blade morphology model and perform steady-state flow field calculations to obtain the flow field distribution in the compressor channel after icing, and then analyze the aerodynamic losses caused by icing.

[0009] In one implementation, step S1 includes: importing the compressor blade profile, casing profile, and hub profile data into the blade modeling module to generate a compressor blade geometric model; importing the compressor blade geometric model into the geometry processing module to generate the fluid domain model, wherein the rotating fluid domain and the stationary fluid domain are connected through an interface.

[0010] In one implementation, in step S2, the fluid domain model is divided into unstructured tetrahedral meshes, and the leading edge region of the blade is locally meshed.

[0011] In one implementation, in step S2, flow field calculations are performed based on multiple grid schemes with different densities. Grid independence is verified by comparing the convergence of key flow field parameters, and the grid used for subsequent calculations is selected based on the verification results.

[0012] In one implementation, in step S2, the boundary conditions for the steady-state flow field calculation include: setting the inlet as a mass flow rate inlet, setting the outlet as a free flow outlet, setting the walls on both sides of the rotor and stator fluid domain as rotational periodic boundaries, setting the interface between them as Interface, and setting the remaining walls as no-slip boundaries.

[0013] In one embodiment, the inlet mass flow rate is set to 0.561 kg / s, the rotor speed is set to 1800 rad / s, and the fluid medium is air.

[0014] In one implementation, in step S3, the UDF used to simulate the behavior of ice crystal particles includes at least: an impact model written based on the DEFINE_DPM_EROSION macro, a particle state update model written based on the DEFINE_DPM_SCALAR_UPDATE macro, a particle motion law model written based on the DEFINE_DPM_LAW macro, and a source term model written based on the DEFINE_DPM_SOURCE macro.

[0015] In one implementation, in step S4, the dynamic mesh method is applied to the fluid domain where the stationary blade is located. The motion law of the mesh nodes is controlled by UDF. The UDF drives the corresponding mesh nodes to move according to the amount of ice formation of the mesh cells within the current time step, thereby realizing the real-time update of the blade morphology.

[0016] In one implementation, in steps S2 and S5, the steady-state flow field calculation employs the multiple reference frame (MRF) method to handle the rotational motion of the rotor region.

[0017] The beneficial effects of this application are as follows: This application, building upon existing methods for calculating ice crystal formation in aero-engine compressor passages, solves the challenge of aligning the amount of ice formation within a grid cell with changes in grid cell size by developing and loading UDF-related programs. This enables real-time updates of the compressor blade morphology after icing, effectively addressing the difficulty in reconstructing the morphology of iced compressor blades. It enriches and improves the research on ice crystal formation in aero-engine compressor passages and establishes a comprehensive calculation method for ice crystal formation in these passages. This application provides crucial technical support for subsequent hazard assessment and risk prediction of ice crystal formation in compressor passages. Attached Figure Description

[0018] Figure 1 This is a flowchart illustrating the calculation method for the entire process of ice crystal melting and icing in the compressor passage of an aero-engine according to an embodiment of this application. Figure 2 This is a single-level stage37 model diagram provided in the embodiments of this application; Figure 3 This is a mesh diagram of a single-level stage37 model provided in the embodiments of this application; wherein, a is the stage37 model mesh, and b is the local refinement of the stage37 model mesh; Figure 4 This is a mesh independence analysis diagram of an embodiment of this application; Figure 5 This is a pressure distribution diagram of a single-stage stage37 model according to an embodiment of this application; Figure 6This is a 30s icing distribution diagram of the single-stage stage37 model in an embodiment of this application; where a is the icing distribution of stage37; and b is the leaf morphology after 20s icing. Figure 7 This is a pressure comparison diagram of the single-stage stage37 model before and after freezing in an embodiment of this application; where a is the pressure distribution before freezing; and b is the pressure distribution after freezing for 20 seconds. Detailed Implementation

[0019] The technical solution of this application will be clearly and completely described below with reference to specific embodiments. However, those skilled in the art will understand that the embodiments described below are only some embodiments of this application, not all embodiments, and are only used to illustrate this application, and should not be regarded as limiting the scope of this application. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0020] To address the limitations of existing numerical simulation methods that rely on fixed geometry and fail to reflect the dynamic processes of icing and real-time flow field feedback, this application establishes a coupled computational flow mechanism within the framework of computational fluid dynamics, integrating ice crystal dynamics, phase change thermodynamics, and geometric evolution. By introducing user-defined functions (UDFs), the impact, adhesion, and freezing processes of ice crystals are accurately described, and the resulting local ice accumulation is calculated. Simultaneously, the calculated ice accumulation in the mesh cells is converted into boundary conditions driving the displacement of mesh nodes. Dynamic meshing technology is used to allow the blade surface mesh to deform in real-time and accordingly as ice grows. This transforms the traditional "one-way" icing prediction into a dynamic simulation of the interaction between "geometry-flow field-ice crystals." Therefore, the altered aerodynamic shape after icing can be directly used for subsequent flow field analysis, thus obtaining a complete physical picture from initial icing to performance evaluation within a coherent numerical experiment.

[0021] In one specific embodiment, this application provides a simulation calculation method for the entire process of ice crystal melting and icing in the compressor passage of an aero-engine, comprising the following steps: S1. Based on the compressor's blade profile, casing profile, and hub profile data, establish a compressor blade geometric model that includes at least one stage of rotor blades and one stage of stator blades.

[0022] In some embodiments, geometric modeling is performed using a specialized computer-aided engineering (CAE) software platform. First, the compressor's airfoil, casing, and hub profile data are imported into a dedicated blade modeling module. Based on the input profile data, the blade modeling module automatically generates a three-dimensional compressor blade geometric model. This model includes at least one rotor blade and one stator blade, constituting a complete "stage."

[0023] Subsequently, the blade geometry model is imported into the geometry processing module. Within this module, the fluid computational domain, or fluid domain model, is generated around the blade entity. This process divides the continuous fluid space into two parts: a rotating fluid domain surrounding the rotor blades and a stationary fluid domain surrounding the stator blades. To simulate the relative motion and flow field transfer between the rotor and stator in actual flow, an "interface" is set at the adjacent interface of these two fluid domains in the software. This interface allows data to be transferred across regions during computation while maintaining the geometric independence of the two fluid domains, thus enabling the definition of different motion properties.

[0024] In some embodiments, the publicly disclosed single-stage Stage 37 compressor is selected as the object. During implementation, three data files—the airfoil profile, upper casing profile, and hub profile corresponding to the Stage 37 blades—are imported into the Bladegen module of the ANSYS Workbench platform. This module generates a three-dimensional geometric model of the Stage 37, containing one rotor blade and one stator blade. During modeling, specific geometric parameters can be set, such as the axial distance from the leading edge of the rotor blade to the inlet section, the axial distance from the trailing edge of the stator blade to the outlet section, and the tip clearance.

[0025] S2. Generate a fluid domain model containing a rotating fluid domain and a stationary fluid domain based on the compressor blade geometric model; wherein the rotating fluid domain and the stationary fluid domain correspond to the regions where the rotor blade and stator blade are located, respectively.

[0026] Generating a fluid domain model that includes both rotating and stationary fluid domains involves: using the geometry processing capabilities of computer-aided engineering (CAE) software to create a continuous cavity entity that completely encloses all blade entities (including rotor blades and stator blades); this cavity serves as the computational region for fluid flow. Subsequently, based on the motion properties of the blades themselves, this continuous fluid domain is logically and geometrically divided.

[0027] Specifically, in the software, the fluid region surrounding the rotor blades is defined as a separate entity, namely the rotating fluid domain. The rotating fluid domain simulates the flow conditions that rotate with the rotor. Correspondingly, the fluid region surrounding the stator blades is defined as another separate entity, namely the stationary fluid domain. The flow within the stationary fluid domain is considered to be stationary relative to the ground coordinate system.

[0028] To ensure the correct transfer of flow parameters (such as pressure and velocity) between the two domains in the calculation, an interface is defined on their shared, virtual geometric interface. This interface is not a solid wall, but a computational boundary that allows for information exchange. The resulting fluid domain model geometrically represents the compressor passage and physically defines the boundaries between the rotating and stationary regions, laying the foundation for subsequent settings of different computational models (such as the multiple reference frame method).

[0029] In some embodiments, after importing the compressor blade geometry model, including the stage37 moving and stationary blades, into the Geometry module of ANSYS Workbench, the fill command is used to generate a fluid domain entity that encloses the blades. Subsequently, the split tool is used to divide this fluid domain entity into two independent volumes, using the axial clearance plane between the moving and stationary blades as the boundary. One volume is named the rotating fluid domain and is specified to contain the moving blade, while the other volume is named the stationary fluid domain and is specified to contain the stationary blade. The contact surface between the two is labeled "Interface" in the software, serving as the interaction surface for subsequent flow field calculations.

[0030] S3. Mesh the fluid domain model and import it into computational fluid dynamics software to perform steady-state flow field calculations and obtain the initial flow field distribution within the compressor channel.

[0031] In some embodiments, unstructured tetrahedral meshes are used to partition the fluid domain model. Unstructured meshes are better suited to complex geometries than structured meshes. Specifically, the fluid domain model is imported into dedicated mesh generation software (such as ANSYS Fluent Meshing), and the software algorithm automatically generates an unstructured tetrahedral mesh that fills the entire fluid domain.

[0032] In some embodiments, considering the large airflow velocity gradient and drastic flow changes in the leading edge region of the blade, local mesh refinement is implemented in this region to accurately capture the flow field details. Specifically, in the mesh generation software, a geometric region near the leading edge of the blade is defined, and the target mesh size in this region is specified to be smaller than that in other regions. The software then generates a mesh that transitions from dense to sparse.

[0033] In some embodiments, grid independence verification is performed to ensure that the final flow field calculation results are not excessively affected by the number and density of grids.

[0034] Specifically, based on the same fluid geometry model, multiple mesh schemes with different global size settings are generated, each with a different total number of elements (density). These meshes are then imported into the flow field solver, and steady-state flow field calculations are performed under identical physical models and boundary conditions. After the calculations are completed, key flow field parameters obtained from each mesh scheme are compared and analyzed, such as the total pressure ratio at the compressor inlet and outlet, isentropic efficiency, or pressure distribution at specific locations on the blade surface. When the mesh is further refined (increasing the total number of elements), if the change in these key parameters is less than a pre-set threshold (e.g., 1%), the calculation results are considered to be largely unaffected by the number of elements, achieving "mesh independence." At this point, one mesh scheme that achieves a balance between computational accuracy and time consumption can be selected as the final mesh scheme for all subsequent calculations.

[0035] In some embodiments, the finalized mesh model is imported into a computational fluid dynamics (CFD) solver (such as ANSYS Fluent) for steady-state flow field calculations. Steady-state calculations assume that the physical quantities at each point in the flow field do not change over time, making them suitable for obtaining stable operating point flow fields.

[0036] Specifically, steady-state flow field calculations involve solvers and models, multi-reference frame processing, boundary conditions, and material properties.

[0037] Solver and Model: The energy equation is enabled to calculate the temperature field; the influence of gravity is considered; an appropriate turbulence model is selected based on the flow characteristics inside the compressor, such as the SST k-omega model, which is relatively accurate in predicting the separation of flow against the pressure gradient; and the corresponding wall function is selected to handle near-wall flow.

[0038] Multi-reference frame processing: To simulate the rotation of the rotor blades, the multi-reference frame (MRF) method is employed. The rotating fluid domain is set as the rotating reference frame and assigned a rotational speed about its axis; the stationary fluid domain remains in the stationary reference frame. Flow field data is exchanged and matched between the two domains through a predefined interface.

[0039] Boundary conditions include: Inlet: Set as mass flow inlet, specifying the mass flow rate flowing into the flow field.

[0040] Outlet: Set to free flow outlet (Pressure Outlet), given static pressure or using the default settings, allowing the flow field to develop freely.

[0041] Walls: The two circumferential walls of the rotor and stator fluid domains are set as rotational periodic boundaries to simulate the real situation of multiple blades throughout the cycle. All solid walls, including rotor blades, stator blade surfaces, casing, and hub, are set as no-slip boundary conditions.

[0042] Interface: The interface type between the rotating domain and the stationary domain is set to Interface to ensure data transfer.

[0043] Material properties: Define the fluid domain medium as air and specify its physical properties such as density, specific heat capacity, thermal conductivity, and viscosity. These parameters can be set to vary with temperature or be constant (such as the physical properties of air at 15℃).

[0044] S4. Based on the initial flow field, enable the Discrete Phase Model (DPM), load the user-defined function (UDF) to simulate the melting, impact and rebound behavior of ice crystal particles, simulate the motion and phase change process of ice crystal particles after they are sucked into the compressor channel, and obtain the distribution and state of ice crystal particles on the blade surface.

[0045] Enable the Discrete Phase Model (DPM) in the flow field solver. DPM treats ice crystal particles as discrete particles, tracking their motion in Lagrangian coordinates, while considering the two-way coupling (momentum and energy exchange) between the particles and the continuous air phase. The calculation needs to be converted from steady-state to transient to simulate the transport process of particles over time.

[0046] In some embodiments, the initial injection conditions of ice crystal particles need to be defined in the DPM settings, including injection location, particle size and distribution, physical properties and concentration or mass flow rate.

[0047] Specifically, the injection location is typically set at the compressor inlet section; the particle size and distribution are the median diameter of the specified particles (e.g., 0.01 mm) and possible distribution patterns (e.g., Rosin-Rammler distribution); the physical properties are the particle's material density, specific heat capacity, latent heat, etc.; the initial state is usually set as solid ice crystals, and an initial temperature and initial velocity are assigned. The concentration or mass flow rate is the total mass of injected particles per unit time, reflecting the ice crystal concentration in the cloud.

[0048] However, the standard DPM model cannot directly describe the complex behaviors unique to ice crystals, such as partial melting, adhesion, rebound, or breakage that may occur upon impact with walls, and the mass and energy source terms caused by phase transitions. Therefore, the model's functionality is extended by writing and loading user-defined functions (UDFs).

[0049] In some embodiments, the UDF contains multiple sub-models, which are developed by calling the predefined macro function interface provided by the solver. Specifically: The impact model, written using the DEFINE_DPM_EROSION macro, defines the fate of ice crystal particles when they come into contact with the blade or flow channel wall. The UDF calculates a capture rate based on local conditions at the moment of particle impact (such as velocity, angle, wall temperature, and liquid water content), determining whether the particle completely adheres, partially adheres and bounces, or bounces completely.

[0050] The particle state update model is written using the DEFINE_DPM_SCALAR_UPDATE macro. This model updates a custom scalar for each tracked particle at each computation time step. These scalars can be used to record key state parameters of the particles, such as the current temperature, the mass fraction of melted liquid (melting ratio), or the total mass change.

[0051] The particle motion law model, written using the DEFINE_DPM_LAW macro, is used to define or modify the force laws acting on ice crystal particles. For example, it adjusts their drag coefficient under specific conditions (such as partial melting leading to shape changes), thereby more realistically affecting their motion trajectory.

[0052] The source term model, written using the DEFINE_DPM_SOURCE macro, calculates the mass, momentum, and energy source terms transferred to the continuous phase air by ice crystal particles during their motion due to phase transitions (melting or freezing). For example, the heat absorbed by melting ice crystals creates a cooling source term in the air; simultaneously, the decrease in particle mass is reflected in the mass source term.

[0053] After compiling the aforementioned UDF and mounting it into the DPM calculation, transient calculation is initiated. Within each time step, the solver simultaneously solves the equations of motion and state for both the continuous phase flow field and the discrete phase particles. This coupled calculation yields key results characterizing the transport and transformation processes of ice crystal particles, including the complete trajectory of particles from the inlet to the outlet or wall impact point, the quantity, mass, or energy distribution of particles impacting the stator and rotor blade surfaces, and the individual state parameters of each particle, such as survival, melting, and the specific melting ratio. These results collectively form the quantitative basis for analyzing ice initiation and growth.

[0054] S5. Based on the distribution and state of ice crystal particles on the blade surface, an icing thermodynamic model is established using a loaded UDF. The surface mass, momentum and energy conservation equations are solved to calculate the icing rate and icing amount on the blade surface. Combining the dynamic mesh method, the positions of the mesh nodes on the blade surface are updated according to the icing amount to reconstruct the morphology model of the compressor blade after icing.

[0055] In some embodiments, the ice growth calculation process is implemented by loading a specially written user-defined function (UDF) into the solver, which embeds an ice thermodynamic model. This model physically models the heat and mass transfer processes on the ice surface based on information about ice crystal particles impacting the blade surface (such as impact mass, location, and state) and local continuous phase flow field conditions (wall temperature, pressure, and shear force).

[0056] Specifically, the model code in the UDF solves a coupled system of surface mass, momentum, and energy conservation. During the calculation, the model considers multiple physical mechanisms: the direct freezing of trapped solid ice crystals, the overflow and redistribution of a possible liquid water film (formed by partially melted ice crystals) under surface shear forces, the phase change of the water film due to evaporation or secondary freezing, and the convection and heat conduction between the solid / liquid interface and the air.

[0057] By numerically solving these interrelated processes, UDF calculates the net icing rate (the increase in ice thickness per unit time and unit area) for each micro-element region (typically corresponding to a grid surface) on the blade surface within each time step. Integrating this rate over time yields the total ice volume (expressed as ice thickness or added mass) for that region during the current total computation time.

[0058] To reflect the feedback effect of geometric changes caused by icing on the subsequent flow field and ice crystal impact, the geometry of the computational model is changed in real time using dynamic mesh technology based on the calculated icing amount.

[0059] In some embodiments, dynamic mesh updates are primarily applied to the fluid domain containing stationary blades (such as stators). The implementation process is as follows: Driven by UDF, write another part of UDF code to read the amount of ice (such as the equivalent ice growth height) of each surface grid cell calculated by the aforementioned icing thermodynamic model within the current time step.

[0060] Node displacement mapping (UDF) translates the required ice growth height for each surface mesh cell into a displacement command that drives all nodes of that cell to move outward along the surface normal. The displacement is typically proportional to the amount of ice and ensures a smooth transition to avoid mesh distortion.

[0061] Dynamic mesh execution involves activating the dynamic mesh model in the solver and specifying the aforementioned UDF as a user-defined function controlling the movement of nodes in the mesh region (i.e., the stator blade surface and its adjacent volume mesh layers). At the end of each transient computation time step, the solver calls this UDF to move the relevant mesh nodes to new positions based on the latest icing data.

[0062] The morphology update and loop: after the mesh nodes are displaced, the geometry of the blade surface is updated, becoming thicker or changing shape to simulate ice growth. The updated mesh is then used for flow field solving and ice crystal tracking in the next time step, thus forming a dynamic closed-loop simulation of flow field-impact-icing-geometric update-new flow field.

[0063] S5. Mesh the reconstructed compressor blade morphology model and perform steady-state flow field calculations to obtain the flow field distribution in the compressor channel after icing, and then analyze the aerodynamic losses caused by icing.

[0064] Because the blade surface geometry changes due to icing, the original computational mesh is no longer applicable. A new mesh needs to be generated based on the exported post-icing blade geometry model.

[0065] In some embodiments, the geometry file containing the updated topography is imported into mesh generation software (such as ANSYS Fluent Meshing), and a new unstructured computational mesh is generated using a strategy similar to that in step S3. It is understood that, to ensure computational consistency and reduce errors introduced by the mesh itself, the type of the newly modeled mesh, global size control, and near-wall treatment principles should be consistent with the initial mesh scheme.

[0066] The newly generated mesh model is imported into the CFD solver to perform steady-state flow field calculations after icing. The settings for this calculation, including the turbulence model, wall functions, and fluid medium properties, are exactly the same as those used in step S3 to obtain the initial flow field, ensuring comparability of the results. Boundary conditions, such as inlet mass flow rate, outlet static pressure, and rotor speed, also need to remain unchanged. For calculations involving rotating components, the multiple reference frame (MRF) method is still used to handle the rotational motion of the rotor fluid domain. The only change is the geometry on which the calculation is based; in this case, it is a blade with an icing profile.

[0067] After the steady-state flow field calculations following icing have been completed and converged, the complete flow field distribution within the compressor passage under icing conditions can be obtained, including the pressure, velocity, and temperature fields. By comparing this flow field with the initial flow field obtained in step S3 before icing, the aerodynamic impact caused by icing can be quantitatively assessed. The focus of the analysis is typically on pressure loss.

[0068] In some embodiments, the changes in static pressure and total pressure at the same axial position or the same blade surface measurement point are compared before and after icing; the degree of decline in overall performance parameters of the compressor stage, such as total pressure rise, pressure ratio, or isentropic efficiency, is calculated and compared before and after icing; flow field cloud maps or vector maps are observed to analyze whether more severe flow separation, low-speed regions, or vortex structures are induced in key areas such as the blade suction surface and trailing edge after icing.

[0069] For example, by extracting and comparing the static pressure distribution curves along the flow channel centerline before and after icing, the additional pressure drop region and its magnitude caused by ice blockage or streamline disturbance can be clearly identified. Finally, the analysis results are presented in graphical form, clearly showing the specific loss patterns of compressor aerodynamic performance caused by icing.

[0070] Example Reference Figure 1As shown, a calculation method for the entire process of ice crystal melting and icing in the compressor passage of an aero-engine is presented, which focuses on solving the problems of difficulty in reconstructing the compressor blade morphology after icing and difficulty in assessing the aerodynamic losses after icing. The method includes the following steps: Step S1: Import the casing profile data, hub profile data, and two blade (moving blade and stationary blade) profile data corresponding to the stage37 model into the Bladegen module under the ANSYS Workbench platform to generate the stage37 model including the moving blade and stationary blade. The distance from the leading edge of the moving blade to the inlet is 100mm, the distance from the trailing edge of the stationary blade to the outlet is 50mm, and the blade tip clearance is 2mm. Step S2: Import the stage37 model containing the moving and stationary blades into the Geometry module of the ANSYS Workbench platform to generate a single-level stage37 model (containing one moving blade and one stationary blade) fluid domain. Then, import the fluid domain into Spaceclaim and name each boundary condition. In addition, the fluid domain containing the moving blade and the fluid domain containing the stationary blade need to be divided into two entities to facilitate subsequent meshing of the two parts separately. Figure 2 As shown; Step S3: The fluid domain models containing both moving and stationary blades are imported into Fluent Meshing for tetrahedral mesh generation. Local mesh refinement is applied to the leading edge of the blades to ensure Y+>1. Six models with mesh sizes of 394W, 550W, 634W, 780W, 906W, and 1018W are obtained for mesh independence analysis. Based on the results, a mesh with a size of 780W is ultimately selected. (Specific details are as follows...) Figure 3 and Figure 4 As shown; Step S4: Import the two meshed fluid domain models into ANSYS Fluent for solving the temperature and velocity fields. During the solution process, first click "Scale" to confirm the correct size of the stage37 model. Next, enable the energy equation and gravity options. Then, select the turbulence model as SST k-omega. Next, set the boundary conditions and computational domains. Finally, begin the steady-state flow field calculation for the stage37 compressor blade model and obtain the temperature and pressure field distribution patterns, as shown below. Figure 5 As shown.

[0071] The specific boundary conditions and computational domain settings are as follows: the walls on both sides of the moving blade fluid domain and the stationary blade fluid domain are set as rotational periodic boundaries; the interface between the moving blade fluid domain and the stationary blade fluid domain is set as Interface; the inlet is set as a mass flow inlet with a value of 0.561 kg / s; the outlet is set as a free flow outlet; the medium in the fluid domain is air at 15°C; the moving blade computational domain uses the multiple reference frame (MRF) method to simulate the blade rotation, with the rotational speed set to 1800 rad / s; all other wall boundaries are set as no-slip boundaries.

[0072] Step S5: After completing step S4, firstly, the steady-state solution used in S4 is changed to a transient solution, with a time step of 0.01s. Then, the discrete phase DPM model is opened, and the particle diameter is set to 0.01mm, the temperature to -10℃, the velocity to 167m / s, the particle shape to be spherical, and the particle material to be water, etc. In addition, an ice accumulation model, scalar, and source-related UDFs need to be loaded into the DPM, and UDFs related to the ice crystal motion model need to be loaded onto the moving blade. Finally, the calculation time step is set to 2000. After the calculation is completed, the ice crystal concentration distribution on the surface of the stator blade can be obtained. Step S6: After completing step S5, turn off the DPM model and turn on the dynamic mesh option. Select the stationary blade as the dynamic mesh region. The mesh motion is controlled by the loaded UDF. Finally, under the condition of a time step of 0.01s, calculate 2000 time steps to obtain the compressor blade morphology 20s after icing. Figure 6 As shown; Step S7: Export the compressor blade model after 20 seconds of icing, then import the exported .cas file into Fluent Meshing to re-mesh it. The size and number of mesh cells should remain basically the same as before. Finally, import the re-meshed compressor blade model into Fluent for calculation, and repeat step S4 to obtain the pressure and velocity distribution in the compressor passage after 20 seconds of icing, thereby obtaining the pressure loss analysis before and after icing. Figure 7 As shown.

[0073] Although the embodiments of this application have been described above in conjunction with the accompanying drawings, this application is not limited to the specific embodiments and application fields described above. The specific embodiments described above are merely illustrative and instructive, not restrictive. Those skilled in the art can make many other forms based on the guidance of this specification and without departing from the scope of protection of the claims of this application, and these are all within the scope of protection of this application.

Claims

1. A simulation calculation method for the entire process of ice crystal melting and icing in the compressor passage of an aero-engine, characterized in that, Includes the following steps: S1. Based on the compressor's blade profile, casing profile, and hub profile data, establish a compressor blade geometric model containing at least one stage of rotor blades and one stage of stator blades, and generate a fluid domain model containing a rotating fluid domain and a stationary fluid domain based on the compressor blade geometric model; wherein the rotating fluid domain and the stationary fluid domain correspond to the regions where the rotor blades and stator blades are located, respectively. S2. Mesh the fluid domain model and import it into computational fluid dynamics software to perform steady-state flow field calculations and obtain the initial flow field distribution in the compressor channel; S3. Based on the initial flow field, enable the Discrete Phase Model (DPM), load the user-defined function (UDF) to simulate the melting, impact and rebound behavior of ice crystal particles, simulate the motion and phase change process of ice crystal particles after they are sucked into the compressor channel, and obtain the distribution and state of ice crystal particles on the blade surface. S4. Based on the distribution and state of ice crystal particles on the blade surface, an icing thermodynamic model is established using the loaded UDF. The surface mass, momentum and energy conservation equations are solved to calculate the icing rate and icing amount on the blade surface. Combining the dynamic mesh method, the position of the mesh nodes on the blade surface is updated according to the icing amount to reconstruct the compressor blade morphology model after icing. S5. Mesh the reconstructed compressor blade morphology model and perform steady-state flow field calculations to obtain the flow field distribution in the compressor channel after icing, and then analyze the aerodynamic losses caused by icing.

2. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 1, characterized in that, Step S1 includes: importing the compressor blade profile, casing profile, and hub profile data into the blade modeling module to generate a compressor blade geometric model; importing the compressor blade geometric model into the geometry processing module to generate the fluid domain model, wherein the rotating fluid domain and the stationary fluid domain are connected through an interface.

3. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine as described in claim 1, characterized in that, In step S2, the fluid domain model is divided into unstructured tetrahedral meshes, and the leading edge region of the blade is locally meshed.

4. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 3, characterized in that, In step S2, flow field calculations are performed based on multiple grid schemes with different densities. Grid independence is verified by comparing the convergence of key flow field parameters, and the grid used for subsequent calculations is selected based on the verification results.

5. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 1, characterized in that, In step S2, the boundary conditions for the steady-state flow field calculation include: setting the inlet as a mass flow rate inlet, setting the outlet as a free flow outlet, setting the walls on both sides of the rotor and stator fluid domains as rotational periodic boundaries, setting the interface between them as Interface, and setting the remaining walls as no-slip boundaries.

6. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 5, characterized in that, The inlet mass flow rate is set to 0.561 kg / s, the rotor speed is set to 1800 rad / s, and the fluid medium is air.

7. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 1, characterized in that, In step S3, the UDF used to simulate the behavior of ice crystal particles includes at least: an impact model written based on the DEFINE_DPM_EROSION macro, a particle state update model written based on the DEFINE_DPM_SCALAR_UPDATE macro, a particle motion law model written based on the DEFINE_DPM_LAW macro, and a source term model written based on the DEFINE_DPM_SOURCE macro.

8. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 1, characterized in that, In step S4, the dynamic mesh method is applied to the fluid domain where the stationary blade is located. The motion law of the mesh nodes is controlled by UDF. The UDF drives the corresponding mesh nodes to move according to the amount of ice formation of the mesh cells within the current time step, thereby realizing the real-time update of the blade morphology.

9. The simulation calculation method for the entire process of ice crystal melting and icing in the compressor channel of an aero-engine according to claim 1, characterized in that, In steps S2 and S5, the steady-state flow field calculation uses the multiple reference frame (MRF) method to handle the rotational motion of the rotor region.

Images