A numerical method for three-dimensional ice shape simulation based on the synchronous solution of the rotating and stationary domains
The three-dimensional ice shape simulation method, which solves the problem of simulating the coupled icing process of rotating and stationary components by synchronously solving the rotation-stationary domain, solves the problem and achieves high-precision ice shape prediction, supporting the anti-icing and de-icing design of complex mechanical equipment.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- NORTHWESTERN POLYTECHNICAL UNIV
- Filing Date
- 2026-02-03
- Publication Date
- 2026-06-23
AI Technical Summary
Existing technologies cannot accurately simulate the three-dimensional icing process of rotating and stationary components in a real coupled flow field, resulting in blind spots in the anti-icing and de-icing design and safety assessment of rotating machinery.
A three-dimensional ice shape simulation numerical calculation method with simultaneous solution in the rotating and stationary domains was adopted. The control equations of water droplet motion in the rotating and stationary coordinate systems were established respectively. The centrifugal force and Coriolis force of the rotation effect were considered. Combined with the air flow field and water droplet impact characteristics, a three-dimensional icing thermodynamic model was constructed, and the ice shape growth was updated by dynamic mesh technology.
It realizes synchronous integrated numerical simulation of rotating and stationary components under real coupled flow fields, improves the accuracy and physical fidelity of icing prediction, and provides a reliable analysis tool for the anti-icing and de-icing design of complex mechanical equipment.
Smart Images
Figure CN122263704A_ABST
Abstract
Description
Technical Field
[0001] This application belongs to the field of icing simulation technology for rotating machinery, specifically involving a numerical calculation method for three-dimensional ice shape simulation by simultaneous solution of rotational and stationary domains. Background Technology
[0002] Icing poses a serious threat to rotating machinery that relies on aerodynamic performance, a challenge particularly prominent in the aviation and wind power sectors. In aviation, aircraft icing is a key factor endangering flight safety: when an aircraft passes through clouds containing supercooled water droplets, the droplets impact the aircraft surface and freeze, significantly altering the aerodynamic shape, leading to decreased lift, increased drag, and severely deteriorated handling and stability. For transport aircraft using turboprop engines, the icing risk is even more complex—both the rotating propeller and the stationary air intake are exposed to the icing environment, and their icing processes are coupled. While the propeller itself ices, the resulting slipstream significantly alters the flow field and water droplet impact characteristics in the air intake region, greatly increasing safety risks and the difficulty of protection.
[0003] Similarly, in the wind power sector, wind turbines, often deployed in high-altitude, low-temperature, and high-humidity environments to capture greater wind energy, are highly susceptible to blade icing. Blade icing not only accelerates material fatigue and shortens service life but also alters the blade's aerodynamic shape, leading to significant losses in annual power generation. Furthermore, the shedding and falling of ice blocks poses a potential threat to equipment and the safety of on-site personnel. These two areas collectively highlight the icing risks faced by rotating machinery in complex climatic environments and the severe challenges they pose to equipment performance, operational safety, and economic efficiency.
[0004] Currently, icing prediction and protection design for rotating machinery mainly rely on two technical methods: icing wind tunnel testing and numerical simulation. While icing wind tunnel testing can provide realistic ice shape data, it has inherent limitations for complex configurations containing large rotating components. Limited by wind tunnel size, dynamic conditions, and testing technology, it is difficult to conduct full-condition icing tests on full-size propeller-intake integrated models, making it impossible to obtain ice growth data of rotating and stationary components under actual coupled flow fields, resulting in experimental blind spots in protection design.
[0005] In numerical simulation, most existing mainstream icing simulation software and methods are developed for stationary components, and their theoretical basis is based on the model of water droplet motion and ice growth in a stationary coordinate system. For the icing problem of rotating components, existing numerical methods usually adopt simplification strategies, such as applying equivalent incoming flow conditions to the rotating component in the stationary coordinate system, or completely separating the rotating component from the stationary component for independent simulation. These methods fail to accurately characterize the dynamic aerodynamic-water droplet coupling effect between the rotating and stationary domains, especially the influence of the unsteady flow field induced by the rotating component on the water droplet impact characteristics on the surface of the stationary component, and also fail to simulate the effects of rotational centrifugal force and Coriolis force on the flow and freezing process of liquid water film on the component surface. Due to the lack of synchronous and integrated simulation capabilities for the icing process of rotating and stationary components under real physical coupling conditions, existing technologies are unable to accurately predict the three-dimensional ice growth of such complex systems.
[0006] Therefore, existing technologies cannot accurately and synchronously simulate the three-dimensional icing process under conditions where rotating and stationary components coexist and are coupled with each other, which has become a technical bottleneck restricting the refined anti-icing and de-icing design and safety assessment of high-end rotating machinery equipment. Summary of the Invention
[0007] To address the shortcomings of existing technologies in simulating integrated icing of rotating and stationary components under coupled conditions, a three-dimensional numerical calculation method for simulating ice shape is provided, which simultaneously solves the rotational and stationary domains. By establishing and simultaneously solving the control equations for water droplet motion in the rotating and stationary coordinate systems, and comprehensively considering the influence of centrifugal force, Coriolis force, and air shear force generated by the rotational effect on the water droplet impact and overflow process, synchronous and high-precision prediction of ice shape growth in the rotating and stationary domains is achieved. This provides a reliable analytical tool for the anti-icing and de-icing design of complex mechanical equipment containing moving parts.
[0008] To achieve the above technical objectives, this application specifically adopts the following technical solution: In one aspect of this application, a numerical calculation method for three-dimensional ice shape simulation by simultaneous solution of rotational and stationary domains is provided, comprising the following steps: S1. Establish a computational domain for the geometric model containing rotating and stationary components and perform mesh generation. Import the mesh into numerical calculation software and solve the airflow field. Specifically, the Reynolds-averaged Navier-Stokes equations are solved for the computational domain containing the rotating components using the multiple reference frame method, while the Reynolds-averaged Navier-Stokes equations are solved directly for the computational domain containing the stationary components. S2. Based on the Euler method, water droplet phase control equations applicable to both stationary and rotating components are established, and solved using a user-defined scalar UDS, simultaneously obtaining the water droplet impact characteristics on the surfaces of the rotating and stationary components; among which, For a stationary component, the governing equation for the water droplet phase is:
[0009]
[0010] For the rotating component, the governing equation for the water droplet phase is:
[0011]
[0012] in, The volume fraction of the water droplet; The density of the water droplet; The velocity of the water droplet; The relative velocity of the water droplets; Air speed; Relative air velocity; The air-water droplet exchange coefficient; It is the acceleration due to gravity; For Hamiltonian operators; It is the rotational angular velocity; It is a position vector; The numerical diffusion coefficient; S3. Based on the obtained water droplet impact characteristics, and combining the laws of mass conservation and energy conservation, a three-dimensional icing thermodynamic model is constructed. The driving force for the surface water film flow is determined for both rotating and stationary components, and the overflow water direction and flow rate of each surface grid cell are calculated. The three-dimensional icing thermodynamic model is solved to obtain the icing mass of each surface grid cell within each icing time step. With icing height ; S4. For each mesh node, traverse and mark all surface mesh elements sharing that node; calculate the node icing height based on the area-weighted average method. Unit vector in the direction of ice growth :
[0013]
[0014] in, The first The icing height and area of a shared grid cell. The unit vector of the outer normal of the mesh cell. This represents the total number of grid cells sharing this node; S5. Based on the icing height of the node Unit vector in the direction of ice growth Calculate the displacement of the mesh node coordinates and update the node coordinates. :
[0015] in, These are the original coordinates of the node.
[0016] S6. Based on the calculated new node coordinates, update the surface mesh using the dynamic mesh function of the numerical calculation software to generate the current icing geometry. S7. Determine whether the cumulative calculation time has reached the preset freezing time: If it has, then determine the current geometric shape as the final ice shape; if it has not, then use the updated mesh as the new input and return to step S1 for the next round of iterative calculation until the cumulative calculation time reaches the preset freezing time.
[0017] In one implementation, in step S3, the mass conservation equation is:
[0018] in, The mass flow rate resulting from the impact of water droplets; This represents the total mass flow rate of overflow water flowing into the current grid cell; The total mass flow rate of the overflow water flowing out of the current grid cell; The mass flow rate carried away by water evaporation; This represents the icing mass flow rate of the current grid cell; The energy conservation equation is as follows:
[0019] in, This represents the total energy flowing in; This represents the energy generated by the impact on the water droplet; The energy released when water droplets freeze; To prevent the flow of hot and cold air; The heat carried away by water evaporation; This represents the total energy flowing out; This refers to the heat that is carried away by convection and heat exchange.
[0020] In one implementation, in step S3, a freezing coefficient is introduced. Describe the freezing state: ; By the freezing coefficient Solving the equations simultaneously with the mass conservation equation and the energy conservation equation yields the icing mass of each grid cell. and overflow water flow rate.
[0021] In one implementation, step S3 involves determining the driving force for the flow of water film on the surfaces of the rotating and stationary components, respectively. For stationary components, the water film flow velocity equal to air speed ; For rotating components, the water film flow rate Calculated by the following formula:
[0022] in, For water film thickness, The dynamic viscosity of water, The air shear stress experienced by the water film. This is the density of water.
[0023] In one implementation, the water film flow velocity is calculated. With the outward normal vector of the grid cell edge dot product To determine the direction of overflow water: If This indicates that water flows out from that side; if If the side is 0, it means that no water flows out from that side or water flows in.
[0024] In one implementation, in step S3, the icing height of each surface grid cell is... Calculated using the following formula:
[0025] in, This represents the icing mass flow rate of the grid cell. For the freezing time step, The density of ice, This represents the area of the grid cell.
[0026] In another aspect of this application, a numerical calculation system for three-dimensional ice shape simulation with simultaneous solution of rotational and stationary domains is provided, comprising: The airflow field calculation module performs mesh generation and airflow field solution for the model calculation domain containing rotating and stationary components. The solution process includes: solving the Reynolds-averaged Navier-Stokes equations for the rotating domain using the multiple reference frame method, and directly solving the Reynolds-averaged Navier-Stokes equations for the stationary domain. The water droplet impact characteristic calculation module, based on the Euler method, establishes and solves different water droplet phase control equations for the rotating and stationary components respectively. The water droplet phase control equations include control equations for the stationary component that neglect centrifugal force and Coriolis force, and control equations for the rotating component that include centrifugal force and Coriolis force terms. The module also outputs the water droplet impact characteristics on the surfaces of the rotating and stationary components simultaneously. The icing thermodynamics solution module constructs and solves a three-dimensional icing thermodynamics model based on the water droplet impact characteristics. The model includes mass conservation equations and energy conservation equations. Different water film flow driving models are used for rotating and stationary parts to calculate the overflow water direction and flow rate, and the icing mass and icing height of each surface grid unit within the time step are obtained. The node ice shape calculation module traverses all grid nodes, finds all surface grid cells that share the node for each node, and calculates the node ice height and ice growth direction unit vector for the node by weighted averaging based on the area, ice height and outward normal unit vector of each grid cell. The dynamic update and judgment module calculates the displacement of the grid nodes and updates their coordinates to generate a new icing surface geometry based on the output of the node ice shape calculation module. It also determines whether the cumulative calculation time has reached the preset total time, and if not, controls the updated grid to be re-inputted into the air flow field calculation module to start a new round of iterative calculation until the preset total time is reached, and outputs the final ice shape.
[0027] In another aspect of this application, a computer device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the numerical calculation method for simultaneous solution of the rotation-stationary domain in three-dimensional ice shape simulation.
[0028] In another aspect of this application, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the numerical calculation method for simultaneous solution of the rotation-stationary domain in three-dimensional ice shape simulation.
[0029] The beneficial effects of this application are as follows: 1. A synchronous and integrated numerical simulation of the icing process of rotating and stationary components under realistic coupled flow field conditions was achieved. By establishing and simultaneously solving the control equations of water droplet motion in the rotating and stationary coordinate systems, the physical distortion caused by existing methods that separate the two or use equivalent simplifications is overcome, and the influence of the slip flow of the rotating component on the impact characteristics of water droplets on the surface of the stationary component is accurately reflected.
[0030] 2. The role of rotation in the entire freezing process is more accurately characterized. During the droplet motion stage, centrifugal force and Coriolis force terms are introduced into the governing equations; during the ice growth stage, the overflow water treatment model considers the driving effect of centrifugal force on the flow of water film on the surface of rotating parts, thus improving the physical fidelity of the simulation of freezing of rotating machinery.
[0031] 3. A three-dimensional numerical prediction framework for icing, applicable to scenarios of coexistence of rotation and stillness, was established. This method integrates airflow field solving, water droplet impact calculation, a three-dimensional thermodynamic model considering overflow, and an ice shape growth and update module based on dynamic mesh, forming a complete and iterative numerical simulation process. This provides an effective analytical tool for the anti-icing design of complex systems such as propeller-inlet and rotor-fuselage.
[0032] 4. By introducing a node displacement algorithm with mesh area weighted average, a stable and reasonable mapping from surface icing amount to geometric shape update is achieved, ensuring the reliability of mesh deformation and geometric smoothness when performing multi-step ice shape growth calculations on complex three-dimensional surfaces.
[0033] 5. By comparing with typical experimental cases of rotating and stationary components, the ice shape calculated by the method of this application is in good agreement with the experimental data, which verifies its prediction accuracy and engineering applicability, and provides a reliable numerical approach to solve the long-standing problem of predicting icing of rotating machinery coupling. Attached Figure Description
[0034] Figure 1 This is a schematic diagram of the numerical calculation method for three-dimensional ice shape simulation using simultaneous rotation-stationary domain solution in this application. Figure 2 This is a schematic diagram of the ice surface boundary movement method in the numerical calculation method for three-dimensional ice shape simulation of this application; Figure 3 The geometric models (a) of the rotating component and the stationary component (b) used for numerical calculation verification in the embodiments of this application are shown. Figure 4 This is the ice shape verification result of Exp1 in the embodiments of this application; Figure 5 This is the ice shape verification result of Exp2 in the embodiments of this application. Detailed Implementation
[0035] 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.
[0036] To address the challenge of handling the coupling effects between rotating and stationary components in existing numerical methods for icing, this application proposes a novel approach for simultaneous solution within a unified framework. For both the non-inertial frame of reference of the rotating component and the inertial frame of reference of the stationary component, separate governing equations for droplet motion, incorporating rotational effects such as Coriolis force and centrifugal force, are established. These equations are then calculated simultaneously during the numerical solution process, directly yielding the droplet impact characteristics under coupled flow fields. Based on this, a unified three-dimensional icing thermodynamic model is constructed, with differentiated modeling of the overflow driving forces on the rotating and stationary surfaces, respectively incorporating the dominant effects of centrifugal force and air shear force. By converting the calculated local icing amount into the normal displacement of the mesh nodes and using dynamic meshing technology to iteratively update the computational domain geometry and flow field, simultaneous and dynamic prediction of ice growth in both the rotating and stationary domains is achieved.
[0037] In one specific implementation, refer to Figure 1 As shown, this application presents a numerical calculation method for three-dimensional ice shape simulation using simultaneous rotation-stationary domain solution, which specifically includes the following steps: S1. Establish a computational domain for the geometric model containing rotating and stationary components and perform mesh generation. Import the mesh into numerical calculation software and solve the airflow field. Specifically, the Reynolds-averaged Navier-Stokes equations are solved using the multiple reference frame method for the computational domain containing the rotating components, while the Reynolds-averaged Navier-Stokes equations are solved directly for the computational domain containing the stationary components.
[0038] Specifically, based on the physical geometry of the rotating and stationary components, a corresponding three-dimensional computational domain is established and meshed to generate a computational mesh file containing volume elements.
[0039] Import the computational mesh file into the selected numerical computation software (FLUENT). In the software, the computational domain is divided into different regions based on the motion state of the component: the region surrounding the rotating component is defined as the rotational computational domain, and the remaining regions are defined as the stationary computational domain.
[0040] In some embodiments, a multi-reference-frame approach is employed to handle the rotational motion of the rotating computational domain. A reference coordinate system that rotates synchronously with the rotating component is applied to this region, and the governing equations of the fluid are solved in this rotating reference frame. For the stationary computational domain, the governing equations are solved directly in an absolutely stationary reference frame. The governing equations for both parts employ the Reynolds-averaged Navier-Stokes equations to describe the airflow under turbulent conditions.
[0041] To close the Reynolds-averaged equations and simulate turbulence effects, the standard equations are selected. Turbulence Model. This model combines the advantages of different turbulence models in the near-wall and mainstream regions. In the mesh region near the wall, the standard wall function method is used to handle the boundary layer flow, ensuring computational accuracy while reducing the requirements for near-wall mesh resolution.
[0042] After completing the above model and algorithm settings, the solution calculation is performed in the numerical calculation software to obtain a stable airflow field distribution in the computational domain, including detailed information on parameters such as velocity and pressure.
[0043] S2. Based on the airflow field, calculate the water droplet impact characteristics of the rotating and stationary parts, including the following steps: Step S21: This invention uses the Euler method to calculate the water droplet impact characteristics. The Euler method treats the water droplet phase as a continuous phase. Since the particle concentration of water droplets is very small, the influence of the water droplet phase on the air phase is considered negligible. Therefore, based on the stable calculation of the air field, the governing equations of the water droplet phase are solved to obtain the volume integral number and velocity distribution of the water droplets.
[0044] The water droplet phase control equations were solved using the Fluent UDS module, where UDS stands for User-Defined Scalars. The water droplet phase control equations are shown below: For a stationary component, the phase continuity equation and momentum equation for the water droplet are as follows:
[0045]
[0046] For the rotating component, the phase continuity equation and momentum equation for the water droplet are as follows:
[0047]
[0048] In the formula, The volume fraction of the water droplet; The density of the water droplet; The velocity of the water droplet; Air speed; It is the acceleration due to gravity; For Hamiltonian operators; The air-water droplet exchange coefficient; Relative air velocity; It is the rotational angular velocity; To control the volume vector; Let be the relative velocity of the water droplets. According to the velocity composition theorem:
[0049]
[0050] In the formula, The absolute velocity of the water droplet; For the speed of the impact; The air-water droplet exchange coefficient is calculated using the following formula:
[0051] In the formula, Indicates aerodynamic viscosity; The diameter of the water droplet; The density of the water droplet; The resistance function is calculated using the following formula:
[0052]
[0053] In the formula, This is the water droplet resistance coefficient; The relative Reynolds number can be calculated using the following formula:
[0054] In the formula, air density; It is the aerodynamic viscosity.
[0055] To avoid calculation divergence caused by abnormal local water droplet volume integrals during the calculation process, a diffusion term needs to be added to the water droplet phase continuity equation, as shown below: Rotating components:
[0056] Stationary components:
[0057] in, is the numerical diffusion coefficient.
[0058] S3. Based on the obtained water droplet impact characteristics, and combining the laws of mass conservation and energy conservation, a three-dimensional icing thermodynamic model is constructed. The driving force for the surface water film flow is determined for both rotating and stationary components, and the overflow water direction and flow rate of each surface grid cell are calculated. The three-dimensional icing thermodynamic model is solved to obtain the icing mass of each surface grid cell within each icing time step. With icing height .
[0059] The icing thermodynamic model is established based on the conservation of mass and energy of surface grid cells within each icing time step. The icing thermodynamic model follows the basic idea of the Messier model and specifically considers the mass and energy exchange caused by surface water film flow under three-dimensional conditions.
[0060] The core of solving the freezing thermodynamic model lies in determining the flow direction and flow rate of unfrozen liquid water (overflow) within each surface grid cell. The driving force for the surface water film flow is determined separately for rotating and stationary components, considering their different mechanical environments. For stationary components, the main driving force for water film flow is the shear stress exerted by air on the water film surface. For rotating components, the driving force for water film flow is composed of air shear stress and volume forces generated by the rotational effect, including centrifugal force and Coriolis force.
[0061] Based on the above driving force analysis, the governing equation for the overflow velocity of the surface water film is established.
[0062] For the surface of rotating components, the water film flow velocity Solve using the following coupling equations:
[0063] in, For water film thickness, The dynamic viscosity of water, The air shear stress experienced by the water film. This is the density of water.
[0064] Air shear stress The calculation formula is: , Let be the air velocity. The water film flow velocity is determined by both inertial acceleration and shear stress. The inertial acceleration is related to this velocity itself and there is a coupling relationship. In the solution process, it needs to be handled by a coupled iterative solution method.
[0065] For a stationary component, the water film velocity is equal to the air velocity near the wall, i.e. .
[0066] After determining the water film flow velocity, the velocity vector and the outward normal vector of each edge of the grid cell are calculated. dot product Determine the direction of overflow. If This indicates that water flows out from the i-th edge of the current grid cell; if If , it means that there is no overflow water flowing out of the current edge.
[0067] After clarifying the overflow inflow and outflow relationship of each grid cell, mass conservation equations and energy conservation equations are established for each surface grid cell.
[0068] The mass conservation equation is expressed as:
[0069] in, The mass flow rate resulting from the impact of water droplets; This represents the total mass flow rate of overflow water flowing into the current grid cell; The total mass flow rate of the overflow water flowing out of the current grid cell; The mass flow rate carried away by water evaporation; This represents the icing mass flow rate of the current grid cell. All units are kg / s.
[0070] To characterize the freezing state, a freezing coefficient is introduced. It is defined as the ratio of the actual freezing rate to the actual total inflow rate within the grid cell: ; Depend on As can be seen from the definition, The range of values is .when When =1, all water within the grid cell freezes without overflow; this is frost / ice. At this time, some of the water in the grid cell freezes into ice, and some overflows into adjacent grid cells through the outflow surface; the ice formed at this time is open ice. When the value is 0, it means that the water flowing into the grid cell is not frozen, and only overflow occurs.
[0071] For the mesh elements of the icy surface, the mass conservation equation is:
[0072] In the formula, This represents the total energy flowing in; This represents the energy brought by the impact of the water droplet. The latent heat released when water droplets freeze; To prevent icing heat flow, the wall surface is generally insulated during icing calculations. ; The heat carried away by water evaporation; This represents the total energy flowing out; The heat carried away by convective heat transfer is expressed in W.
[0073] Solving the above mass and energy conservation equations simultaneously yields the icing mass of each surface grid cell in the current time step. Combining the area A of the grid cell and the icing time step... And the density of ice The icing height of the grid cell within this time step can then be calculated. : .
[0074] This yields the icing mass distribution and geometric growth height of all mesh cells on the entire icing surface at this time step.
[0075] S4, Reference Figure 2 As shown, for each mesh node, all surface mesh elements sharing that node are traversed and marked; the node icing height is calculated based on the area-weighted average method. Unit vector in the direction of ice growth .
[0076] Since ice growth is a continuous physical process, while numerical models use discrete grid cells for description, each grid node is typically shared by multiple adjacent surface grid cells. To obtain a smooth and reasonable ice growth at a node, the ice growth information of all surface grid cells sharing that node needs to be fused.
[0077] The processing iterates through all surface mesh nodes. For each mesh node to be computed, all surface mesh elements sharing this node are identified and marked. Let the number of marked elements for a given node be... .
[0078] The icing height at the nodes was determined using an area-weighted average method. The icing height of each marked unit was... Its grid cell area Multiply the products of all marked units, sum the products, and then divide this sum by the sum of the areas of all marked units. This calculation is expressed by the following formula:
[0079] In the formula, This represents the calculated icing height of the node. This method shows that larger mesh cells contribute more to the icing height of a node than smaller cells.
[0080] The direction of ice growth is also determined based on the principle of area-weighted averaging, assuming that the ice grows along the surface normal direction. Each marked surface mesh cell has its own unit outward normal vector. The unit vector of the ice growth direction of the node. Calculate using the following steps: .
[0081] S5. Based on the icing height of the node Unit vector in the direction of ice growth Calculate the displacement of the grid node coordinates and update the node coordinates.
[0082] For each surface mesh node, its new coordinates after icing are obtained by calculating its positional displacement. The displacement is calculated based on the icing height of that node. Its unit vector in the direction of ice growth The components of this unit vector in three-dimensional space are represented as follows: .
[0083] The node coordinates are updated according to the following relationship: the node's new coordinates after ice growth. From its original coordinates This is obtained by adding the displacement vector along the growth direction. The magnitude of this displacement vector is the nodal icing height. The direction is determined by the unit vector. definition.
[0084] Specifically, the calculation formulas for each component of the new coordinate system are as follows:
[0085] Through the above calculations, all surface mesh nodes will move along their normal direction according to the local icing amount, thereby updating the geometric description of the part's shape after icing.
[0086] S6. Based on the calculated new node coordinates, update the surface mesh using the dynamic mesh function of the numerical calculation software to generate the geometry of the iced surface at the current moment.
[0087] The dynamic mesh feature allows the positions of mesh nodes to be changed during calculation according to user-defined rules, thereby simulating boundary movement or deformation. To simulate surface deformation caused by ice growth, specific user-defined functions need to be written within this feature framework.
[0088] Specifically, this is accomplished by defining the macro function DEFINE_GRID_MOTION. Within this macro function, the program reads the new coordinate values corresponding to each surface mesh node, calculated in step S5. Write program logic to update the currently stored mesh node coordinates in the software with these newly calculated values. This function controls the dynamic mesh module to move the coordinates of boundary nodes defined as icy walls only; internal mesh nodes are automatically and smoothly adjusted by the software's dynamic mesh algorithm based on boundary changes.
[0089] Once the custom function is written and configured in the software, the dynamic mesh update process is initiated. The software executes the function, applies the displacements of all nodes, and reconnects the mesh topology to generate a new computational mesh that adapts to the updated boundary shape. The geometry corresponding to this new mesh is the shape of the component after icing at the end of the current icing time step.
[0090] S7. Determine whether the cumulative calculation time has reached the preset freezing time: If it has, then determine the current geometric shape as the final ice shape; if it has not, then use the updated mesh as the new input and return to step S1 for the next round of iterative calculation until the cumulative calculation time reaches the preset freezing time.
[0091] Step S7 is used to control the time progression and iteration termination of the icing simulation. This step compares the current cumulative calculation time with the preset total icing time.
[0092] If the result indicates that the cumulative calculation time has reached or exceeded the preset ice summation time, then the geometric shape generated in step S6 is determined as the final ice shape, and the entire numerical simulation process ends.
[0093] If the result indicates that the cumulative calculation time has not yet reached the preset summarization time, iterative calculation for the next time step needs to continue. At this point, the surface mesh information containing the new geometry generated after the dynamic mesh update in step S6 is output as a boundary mesh file.
[0094] The boundary mesh file is imported into the mesh generation software. In the mesh generation software, a new volume mesh suitable for subsequent flow field calculations is generated based on this deformed boundary geometry. After the new volume mesh is generated, it is loaded back into the main numerical computation software.
[0095] Subsequently, the simulation process returns to step S1, using the newly loaded volume mesh as the computational domain, and begins the next round of iterative calculations. Specifically, the airflow field, water droplet impact characteristics, icing thermodynamic processes, nodal displacements, and mesh updates are recalculated. The decision logic in step S7 will be executed again at the new computation time point.
[0096] The cycle of calculation, judgment, mesh update, and recalculation will continue until the cumulative calculation time meets the preset ice-summing time condition. The final output is the icing geometry of the component at the end of the entire time history.
[0097] Example To verify the effectiveness and accuracy of the proposed numerical calculation method for simultaneous solution of rotating and stationary domains in three-dimensional ice shape simulation, a comparative verification was conducted using numerical calculations and wind tunnel experiments. Given the lack of icing data in currently available experiments showing the rotating and stationary components operating under realistic coupling conditions, this embodiment adopts a strategy of verifying the rotating and stationary components separately.
[0098] The geometric model used for verification is as follows: Figure 3As shown, the enlarged portions of the figure are airfoil NACA0018 (a) and airfoil NACA0012 (b). First, the numerical calculation method of this invention is used to simulate the icing test conditions of the rotating components listed in Table 1. Then, the exact same test conditions are reproduced in the icing wind tunnel to obtain the experimental ice shape of the rotating component. Figure 4 The results show a comparison between the ice shapes obtained from numerical calculations and those obtained from wind tunnel experiments.
[0099] Table 1 Icing Conditions of Rotating Components
[0100] Furthermore, to verify the applicability of the method to stationary components, the airfoil icing test conditions listed in Table 2 were selected, and the icing calculations were performed using the numerical method of this invention. Icing wind tunnel experiments were conducted under the same incoming flow conditions to obtain the experimental ice shapes. Figure 5 The results show a comparison between the calculated ice shape and the experimental ice shape of the stationary airfoil.
[0101] Table 2 Icing Conditions for Stationary Components
[0102] Through the comparative verifications conducted on rotating and stationary components respectively, the computational accuracy of the method of this invention can be evaluated when handling rotational effects and icing problems in stationary flow. The verification results show that this method can effectively predict the three-dimensional ice growth on the surface of components in different motion states, laying a reliable numerical foundation for simulating complex coupled icing scenarios where rotating and stationary components coexist.
[0103] 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 numerical calculation method for three-dimensional ice shape simulation using simultaneous rotation-stationary domain solution, characterized in that, Includes the following steps: S1. Establish a computational domain for the geometric model containing rotating and stationary components and perform mesh generation. Import the mesh into numerical calculation software and solve the airflow field. Specifically, the Reynolds-averaged Navier-Stokes equations are solved for the computational domain containing the rotating components using the multiple reference frame method, while the Reynolds-averaged Navier-Stokes equations are solved directly for the computational domain containing the stationary components. S2. Based on the Euler method, water droplet phase control equations applicable to both stationary and rotating components are established, and solved using a user-defined scalar UDS, simultaneously obtaining the water droplet impact characteristics on the surfaces of the rotating and stationary components; among which, For a stationary component, the governing equation for the water droplet phase is: For the rotating component, the governing equation for the water droplet phase is: in, The volume fraction of the water droplet; The density of the water droplet; The velocity of the water droplet; The relative velocity of the water droplets; Air speed; Relative air velocity; The air-water droplet exchange coefficient; It is the acceleration due to gravity; For Hamiltonian operators; It is the rotational angular velocity; It is a position vector; The numerical diffusion coefficient; S3. Based on the obtained water droplet impact characteristics, and combining the laws of mass conservation and energy conservation, a three-dimensional icing thermodynamic model is constructed. The driving force for the surface water film flow is determined for both rotating and stationary components, and the overflow water direction and flow rate of each surface grid cell are calculated. The three-dimensional icing thermodynamic model is solved to obtain the icing mass of each surface grid cell within each icing time step. With icing height ; S4. For each mesh node, traverse and mark all surface mesh elements sharing that node; calculate the node icing height based on the area-weighted average method. Unit vector in the direction of ice growth : in, The first The icing height and area of a shared grid cell. The unit vector of the outer normal of the mesh cell. This represents the total number of grid cells sharing this node; S5. Based on the icing height of the node Unit vector in the direction of ice growth Calculate the displacement of the mesh node coordinates and update the node coordinates. : in, These are the original coordinates of the node. S6. Based on the calculated new node coordinates, update the surface mesh using the dynamic mesh function of the numerical calculation software to generate the current icing geometry. S7. Determine whether the cumulative calculation time has reached the preset freezing time: If it has, then determine the current geometric shape as the final ice shape; if it has not, then use the updated mesh as the new input and return to step S1 for the next round of iterative calculation until the cumulative calculation time reaches the preset freezing time.
2. The numerical calculation method for three-dimensional ice shape simulation according to claim 1, characterized in that, In step S3, the mass conservation equation is: in, The mass flow rate resulting from the impact of water droplets; This represents the total mass flow rate of overflow water flowing into the current grid cell; The total mass flow rate of the overflow water flowing out of the current grid cell; The mass flow rate carried away by water evaporation; This represents the icing mass flow rate of the current grid cell; The energy conservation equation is as follows: in, This represents the total energy flowing in; This represents the energy generated by the impact on the water droplet; The energy released when water droplets freeze; To prevent the flow of hot and cold air; The heat carried away by water evaporation; This represents the total energy flowing out; This refers to the heat that is carried away by convection and heat exchange.
3. The numerical calculation method for three-dimensional ice shape simulation according to claim 2, characterized in that, In step S3, a freezing coefficient is introduced. Describe the freezing state: ； By the freezing coefficient Solving the equations simultaneously with the mass conservation equation and the energy conservation equation yields the icing mass of each grid cell. and overflow water flow rate.
4. The numerical calculation method for three-dimensional ice shape simulation according to claim 1, characterized in that, In step S3, the driving force for the flow of water film on the surface of the rotating and stationary components is determined separately, specifically as follows: For stationary components, the water film flow velocity equal to air speed ; For rotating components, the water film flow rate Calculated by the following formula: in, For water film thickness, The dynamic viscosity of water, The air shear stress experienced by the water film. This is the density of water.
5. The numerical calculation method for three-dimensional ice shape simulation according to claim 4, characterized in that, By calculating the water film flow velocity With the outward normal vector of the grid cell edge dot product To determine the direction of overflow water: If This indicates that water flows out from that side; if If the side is 0, it means that no water flows out from that side or water flows in.
6. The numerical calculation method for three-dimensional ice shape simulation according to claim 1, characterized in that, In step S3, the icing height of each surface mesh element is... Calculated using the following formula: in, This represents the icing mass flow rate of the grid cell. For the freezing time step, The density of ice, This represents the area of the grid cell.
7. A numerical calculation system for three-dimensional ice shape simulation with simultaneous solution of rotational and stationary domains, characterized in that, include: The airflow field calculation module performs mesh generation and airflow field solution for the model calculation domain containing rotating and stationary components. The solution process includes: solving the Reynolds-averaged Navier-Stokes equations for the rotating domain using the multiple reference frame method, and directly solving the Reynolds-averaged Navier-Stokes equations for the stationary domain. The water droplet impact characteristic calculation module, based on the Euler method, establishes and solves different water droplet phase control equations for the rotating and stationary components respectively. The water droplet phase control equations include control equations for the stationary component that neglect centrifugal force and Coriolis force, and control equations for the rotating component that include centrifugal force and Coriolis force terms. The module also outputs the water droplet impact characteristics on the surfaces of the rotating and stationary components simultaneously. The icing thermodynamics solution module constructs and solves a three-dimensional icing thermodynamics model based on the water droplet impact characteristics. The model includes mass conservation equations and energy conservation equations. Different water film flow driving models are used for rotating and stationary parts to calculate the overflow water direction and flow rate, and the icing mass and icing height of each surface grid unit within the time step are obtained. The node ice shape calculation module traverses all grid nodes, finds all surface grid cells that share the node for each node, and calculates the node ice height and ice growth direction unit vector for the node by weighted averaging based on the area, ice height and outward normal unit vector of each grid cell. The dynamic update and judgment module calculates the displacement of the grid nodes and updates their coordinates to generate a new icing surface geometry based on the output of the node ice shape calculation module. It also determines whether the cumulative calculation time has reached the preset total time, and if not, controls the updated grid to be re-inputted into the air flow field calculation module to start a new round of iterative calculation until the preset total time is reached, and outputs the final ice shape.
8. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the three-dimensional ice shape simulation numerical calculation method according to any one of claims 1-6.
9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the numerical calculation method for three-dimensional ice shape simulation as described in any one of claims 1-6.