Aero-engine blade sand dust erosion damage evolution prediction method and system

Abstract

The present application relates to the technical field of numerical simulation and structural damage prediction of aero-engine, and particularly relates to an aero-engine blade sand dust erosion damage evolution prediction method and system. The aero-engine blade sand dust erosion damage evolution prediction method comprises the following steps: obtaining a current airflow field F i , an initial erosion result R and a current erosion result E i of an initial complete blade; determining a current time step Δt i by combining wall thickness calculation and erosion depth evaluation according to a fast prediction method; sequentially performing blade geometry Boolean cutting and blade profile extraction and parameterized reconstruction, and outputting an erosion geometry grid M i+1 ; mapping the current airflow field F i to M i+1 by using a Fluent data interpolation function, taking the airflow field as an initial airflow field for next simulation, obtaining a new airflow field and a new erosion result; when a total blade erosion time reaches a preset value, verifying the current time step Δt i for irrelevance, and outputting a full-range damage evolution result.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation and structural damage prediction technology for aero-engines, specifically to a method and system for predicting the evolution of sand and dust erosion damage on aero-engine blades. Background Technology

[0002] When aero-engines operate in complex environments such as deserts and plateaus, a large number of sand and dust particles carried in the air enter the compressor. These high-speed sand and dust particles continuously impact the leading edge and surface area of ​​the compressor blades, causing wear, erosion, and changes in the blade material's geometry. These changes in blade morphology further affect the local flow field structure and particle trajectory, resulting in a distinct coupling and evolutionary characteristic in the subsequent damage process.

[0003] In existing engineering projects, numerical simulation methods are commonly used to analyze blade damage under sand-swallowing conditions. A common approach is based on coupled fluid dynamics and particle dynamics calculations, which calculate particle trajectories and collision behaviors to obtain the erosion or wear distribution on the blade surface. However, these existing methods generally suffer from the following shortcomings: First, the blade geometry remains unchanged over a long period of time. In most numerical simulations, the blade geometry model remains in its initial state throughout the calculation process, and the damage results are only statistically analyzed as an additional quantity. This fails to reflect the impact of morphological changes on subsequent flow and particle behavior, resulting in a significant decrease in long-term prediction accuracy.

[0004] Second, there is a lack of efficient methods for determining the geometric iteration cycle. Existing methods typically use a fixed time step for geometric updates, which cannot guarantee that each update corresponds to a meaningful morphological change, nor can they adaptively adjust the update frequency, resulting in a large amount of invalid calculations or the omission of key damage states.

[0005] Third, the simulation is re-initialized from zero after each geometry update. Whenever the blade geometry changes, a new round of simulation needs to start iteratively from the zero initial field, which has high convergence cost and low computational efficiency, severely restricting the engineering feasibility of multi-iteration, long-term damage evolution simulation.

[0006] Fourth, the ability to predict long-term damage evolution is insufficient. If calculations are always based on the initial geometry, the prediction results will deviate significantly from the actual service conditions; and existing incremental methods are limited in practicality due to the lack of the three key improvements mentioned above.

[0007] Therefore, there is an urgent need for a method to predict the evolution of blade sand and dust erosion damage that can automatically determine the geometric update time step, efficiently generate post-erosion meshes, and accelerate convergence using historical airflow fields at an acceptable computational cost. Summary of the Invention

[0008] The purpose of this invention is to provide a method for predicting the evolution of sand and dust erosion damage on aero-engine blades, in order to solve the problems of difficulty in determining the geometric update time step, low reconstruction efficiency, and the need to initialize from zero for each simulation in the existing numerical simulation methods.

[0009] To address the aforementioned problems, this invention proposes a method for predicting the evolution of sand and dust erosion damage on aero-engine blades. The technical solution employed is as follows: A method for predicting the evolution of sand and dust erosion damage on aero-engine blades includes the following steps: Step S1: Construct a computational domain mesh using the complete blade's three-dimensional geometric model, and perform CFD discrete-phase particle simulation using Fluent to obtain the current airflow field F of each stage of the initial complete blade. i Given the initial erosion rate R, each blade stage synchronously and sequentially performs erosion rate space mapping to obtain the current erosion result E. i , where i=0; Step S2: Based on the current erosion results E of each blade level i The wall thickness is calculated and the erosion depth is assessed, and the boundary contour of the erosion penetration area is extracted. A fast prediction method is used to determine the time period Δt required when the erosion volume of any first-stage blade reaches a preset damage threshold. i This serves as the current propulsion time step and simultaneously triggers a full-stage blade morphology update. Step S3: Using the boundary contour of the erosion penetration area of ​​each blade stage and the current propulsion time step Δt i As input, each stage of the blade synchronously and sequentially performs blade geometric Boolean cutting, blade profile extraction, and parametric reconstruction, outputting the post-erosion geometric mesh M of each stage of the blade. i+1 ; Step S4: Calculate the current airflow field F for each stage of the blades. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ; Step S5: Determine whether the total erosion time of the blade has reached the preset value, wherein: If the target is not reached, let i = i + 1 and return to steps S2-S5; If the condition is met, the entire damage evolution result will be output.

[0010] Furthermore, in step S1, the erosion rate spatial mapping specifically includes: The STL format three-dimensional geometric model of the blade is aligned with the coordinate units of each erosion point containing spatial coordinates and corresponding initial erosion rates, and spatial alignment is achieved through rotation transformation. The initial erosion rate corresponding to each erosion point is mapped to the mesh vertices of the blade's STL format three-dimensional geometric model to obtain a continuous erosion distribution field.

[0011] Furthermore, the spatial alignment specifically includes: Calculate the principal distribution angles of the mesh vertex coordinate set and the coordinate set of each erosion point in the Z-axis direction for the three-dimensional geometric model of the blade in STL format. Use the difference between the two principal distribution angles as the rotation angle, and determine the optimal alignment angle by taking the average nearest neighbor distance between the mesh vertex coordinate set and the coordinate set of each erosion point after rotation. The formula for calculating the optimal alignment angle is:

[0012] in: For the optimal alignment angle, This is the set of mesh vertex coordinates for the STL format 3D geometric model of the blade. For a single grid vertex in the set of grid vertex coordinates, (This is the 3D coordinate vector of the vertex.) This is the set of coordinates for each erosion point. R(θ) is the three-dimensional coordinate vector of a single erosion point in the set of coordinates of each erosion point, R(θ) is the rotation matrix about the Z-axis, and θ is the rotation angle. The initial erosion rate mapping for each erosion point is obtained using kernel density estimation and then smoothed using Laplace, specifically including: Each erosion point is projected onto the nearest neighbor point on the surface of the blade's STL format 3D geometric model, and erosion points whose projection distance is greater than a preset multiple of the projection distance are removed as outliers to obtain the projected erosion points. For each mesh vertex, a support domain is determined with the mesh vertex as the center and the preset kernel bandwidth radius as the support radius. Within this support domain, the same-side normal condition is applied to filter the mesh vertices and the projected erosion points, resulting in each filtered erosion point and each mesh vertex. The initial erosion rates corresponding to each erosion point after screening are weighted and summed according to the kernel function as the numerator, and the areas of the dual units of each grid vertex after screening are weighted and summed according to the same kernel function as the denominator, to obtain the erosion rate after mapping the grid vertex. The erosion rate of each mesh vertex after mapping is smoothed using Laplacian smoothing to eliminate noise; The formula for the kernel density estimation method is as follows:

[0013] in, For grid vertices Mapped erosion rate For grid vertices The collection of erosion points after center screening. For the erosion points after screening The initial erosion rate at the location, For the erosion points after screening The projected 3D coordinate vector For the erosion points after screening The position of the projected 3D coordinate vector and the mesh vertex The Euclidean distance between them For grid vertices The set of all grid vertices after filtering by the center. for The filtered grid vertices For the filtered grid vertices The area of ​​the dual unit. For the filtered grid vertices With grid vertices The Euclidean distance between them For kernel functions; The pair The calculation formula is:

[0014] in, For the filtered grid vertices The area of ​​the dual unit. To match the filtered grid vertices The set of adjacent triangular facets, for A single triangular facet in the middle, Triangular facet The area; The formula for Laplace smoothing is:

[0015] in, For grid vertices Erosion rate after smoothing treatment For grid vertices Mapped erosion rate For grid vertices a certain adjacent grid vertex Mapped erosion rate Let x be the set of adjacent grid vertices. It is a constant.

[0016] Furthermore, in step S2, the calculation of wall thickness and assessment of erosion depth, and the extraction of the boundary contour line of the erosion penetration area, specifically include: Based on the geometric characteristics of the thin-walled blade structure, the wall thickness of each grid vertex is automatically calculated by anti-normal nearest neighbor matching; the erosion depth is calculated by combining the continuous erosion distribution field, and the boundary contour line of the erosion penetration area is extracted. The wall thickness calculation specifically includes: Pairing point search: For a target vertex q, search for a pairing point Q within its K nearest neighbors that simultaneously satisfies the inverse normal angle condition and the distance projection condition. Wall thickness assignment: If a pairing point Q is found, the Euclidean distance between the target vertex q and the pairing point Q is used as the wall thickness of the target vertex q, and thus the wall thickness field is obtained; if a pairing point Q is not found, the wall thickness field is obtained by interpolating using the known wall thickness in the neighborhood through inverse distance weighted interpolation. Wall thickness field smoothing: Laplace smoothing is applied to the wall thickness field; The condition for the anti-normal included angle is: nq · nQ<-α Where: nq and nQ are the area-weighted normal vectors at the target vertex q and the paired point Q, respectively, and α is a constant; The distance projection condition is: |( q - Q) · nq| / ‖q - Q‖>β Where: β is a constant, ‖q - Q‖ is the Euclidean distance between the target vertex q and the paired point Q, and |( q - Q) ·nq| is the projection length of the position difference vector between the target vertex q and the paired point Q onto the normal vector direction.

[0017] Furthermore, the erosion depth assessment is based on the Finnie erosion model, combined with a continuous erosion distribution field, to obtain the erosion depth of each grid vertex; The erosion depth of each grid vertex is: h(x) = E(x) · T / ρ Where: h(x) is the erosion depth at grid vertex x, E(x) is the erosion rate after smoothing grid vertex x, T is the running time, and ρ is the blade material density; The extracted boundary contour line of the erosion penetration area includes: When the erosion penetration ratio is ≥1.0, it is judged as complete erosion penetration; When the erosion penetration ratio is less than 1.0, it is judged as incomplete erosion penetration; At the interface between the fully eroded and partially eroded areas, the contour lines of the erosion boundary are extracted and smoothed using a contour line tracing algorithm to obtain the boundary contour lines of the eroded area. The erosion penetration ratio is: r(x) = h(x) / t(x) Where r(x) is the erosion penetration ratio, h(x) is the erosion depth at mesh vertex x, and t(x) is the original material thickness at mesh vertex x.

[0018] Furthermore, the fast prediction method specifically includes: Based on the current erosion results of each blade level E i Each grid vertex and its corresponding erosion rate, combined with the blade material density ρ, according to h(x)=E i (x)·Δt / ρ, calculates the linear extrapolation relationship of the erosion volume at each grid vertex in each stage of the blade over time, where h(x) is the erosion depth at grid vertex x, Δt is the erosion time step, and E i (x) represents the current erosion rate at grid vertex x; The time interval required for the maximum erosion volume of any stage blade to first reach 2% of the current blade volume is defined as the current propulsion time step Δt. i .

[0019] Furthermore, the blade geometric Boolean cutting specifically includes: For each erosion-penetrated boundary contour line, calculate the average surface normal direction of the nearest neighbor point. Projecting the contour onto a vertical plane constructs a closed curve and generates a surface, along... Bidirectional stretching by a preset distance generates a cut entity; Perform Boolean subtraction operations on the blade solid model and each cut solid model in sequence, and output the damaged blade solid model.

[0020] Furthermore, the blade profile extraction and parameterized reconstruction specifically include: Based on the percentage of leaf height, each streamline of the damaged blade solid model is fitted as a B-spline curve and a surface of revolution is generated. The damaged blade solid model is then sectioned to obtain a two-dimensional contour point set of each streamline. Connectivity analysis is performed on the two-dimensional contour point set based on the K nearest neighbors. A closed curve is constructed by searching for Eulerian circuits, and adaptive point density adjustment is performed. Calculate the three-dimensional curvature κ at each point on the closed curve. i The maximum point of three-dimensional curvature is searched in the extreme region of the Z coordinate as the dividing point of the front and rear edges, and the closed curve is divided into two open curves, which are automatically labeled as pressure surface profile and suction surface profile according to geometric criteria. After the pressure surface profile and suction surface profile are processed by Gaussian smoothing or Savitzky-Golay filtering, all blade height streamlines are sorted by blade height percentage and assembled and output in geomTurbo format; The three-dimensional curvature κ i The formula is: κ i = 2‖(P i+1 - P i ) × (P i - P i-1 )‖ / (‖P i+1 - P i || · ||P i - P i-1 ‖ · ‖P i+1 - P i-1 ‖) Among them, P i Let P be the three-dimensional coordinates of the i-th point on the closed curve. i+1 Let P be the three-dimensional coordinates of the (i+1)th point on the closed curve. i-1 Let be the three-dimensional coordinates of the (i-1)th point on the closed curve.

[0021] Furthermore, the current airflow field F of each stage of the blades is... i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ,include: Complete the post-erosion geometry M in Fluent i+1 Importing, defining regions, and setting boundary conditions; The velocity, pressure, temperature, and turbulent flow rates after convergence of the CFD discrete-phase particle simulation are exported as nodal data and spatially interpolated onto the post-erosion geometric mesh M. i+1 This serves as the initial airflow field for the next CFD discrete phase particle simulation. Post-erosion geometric mesh M i+1 The interpolated airflow field was used as the initial airflow field for CFD discrete phase particle simulation to obtain the new airflow field F for each stage of the blades. i+1 and new erosion results E i+1 .

[0022] The present invention also provides a system for performing the above-described method for predicting the evolution of sand and dust erosion damage on aero-engine blades, comprising: The initial simulation module is used to construct a computational domain mesh using a complete three-dimensional geometric model of the blade, and to perform CFD discrete-phase particle simulation using Fluent to obtain the current airflow field F of each stage of the blade in the initial complete blade. i Given the initial erosion rate R, each blade stage synchronously and sequentially performs erosion rate space mapping to obtain the current erosion result E. i , where i=0; A rapid prediction module is used to predict erosion results based on the current erosion status of each blade level. i The wall thickness is calculated and the erosion depth is assessed, and the boundary contour of the erosion penetration area is extracted. A fast prediction method is used to determine the time period Δt required when the erosion volume of any first-stage blade reaches a preset damage threshold. i This serves as the current propulsion time step and simultaneously triggers a full-stage blade morphology update. The mesh reconstruction module is used to define the boundary contours of the erosion penetration areas of each blade stage and the current propulsion time step Δt. i As input, each stage of the blade synchronously and sequentially performs blade geometric Boolean cutting, blade profile extraction, and parametric reconstruction, outputting the post-erosion geometric mesh M of each stage of the blade. i+1 ; The airflow field interpolation module is used to interpolate the current airflow field F of each stage of the blades. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ; The iterative control module is used to determine whether the total erosion time of the blade has reached a preset value, wherein: If the target is not reached, let i = i + 1 and return to steps S2-S5; If the condition is met, the entire damage evolution result will be output.

[0023] Compared with the prior art, this application has the following beneficial effects: (i) Introducing an automatic time step determination mechanism for damage threshold determination into blade erosion damage evolution simulation: The threshold for triggering geometric update is set when the erosion volume of any stage of the blade reaches a preset damage threshold. The current propulsion time step Δt is calculated by fast linear extrapolation of the current erosion rate. i Compared with the fixed step size method, this mechanism automatically extends the update interval to reduce invalid calculations during the slow erosion stage and automatically shortens the interval to capture key morphological changes during the rapid erosion stage, adaptively adjusting the morphological update accuracy of the entire blade stage; under the same computational load, more intermediate morphological states of the entire blade stage can be resolved, reducing the damage processes missed due to excessively large time steps.

[0024] (ii) Automatic reconstruction code from simulated erosion results to a computable geometric mesh is provided, resolving the engineering bottleneck of geometric updates: using each erosion point exported from Fluent as the sole input, the code automatically completes erosion rate spatial mapping, wall thickness calculation and erosion depth assessment, blade geometric Boolean cutting, and blade profile extraction and parametric reconstruction sequentially. The final output can be directly imported into a standard profile file from Numeca AutoGrid. The above process requires no manual modeling intervention. Compared to methods relying on manual reconstruction of damage geometry, this method supports automatic and continuous multi-round iterative simulations.

[0025] (iii) Using the previous converged airflow field as the initial airflow field for the next Fluent simulation enables knowledge transfer and rapid convergence across iterative airflow fields: Since the blade morphology changes little and the airflow fields are highly similar between adjacent iterations, using the previous converged airflow field as the initial field can reduce the number of iterations required for simulation convergence by more than 50%, providing efficiency assurance for long-term, multi-cycle full-life evolution simulation.

[0026] (iv) Three major innovations synergistically constitute a closed-loop, efficient evolution prediction framework: a rapid prediction method ensures that each update corresponds to a meaningful morphological change; automatic reconstruction ensures the automation and efficiency of geometric reconstruction; and airflow field interpolation ensures rapid convergence of each CFD discrete phase particle simulation. The above three steps are sequentially linked, forming an iterative prediction process that combines erosion rate extrapolation, automatic geometric reconstruction, and airflow field transfer. Attached Figure Description

[0027] Figure 1 A schematic flowchart of the method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to the present invention; Figure 2 A flowchart illustrating steps S1, S2, and S3 in the method for predicting the evolution of sand and dust erosion damage to aero-engine blades according to the present invention. Figure 3 The erosion distribution cloud map drawn based on simulation results in the prediction method for sand and dust erosion damage evolution of aero-engine blades in this invention; Figure 4 A schematic diagram of the grid generated in the method for predicting the evolution of sand and dust erosion damage to aero-engine blades in this invention. Detailed Implementation

[0028] The embodiments of this application are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain this application, and should not be construed as limiting this application.

[0029] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0030] The following describes, with reference to the accompanying drawings, a method and system for predicting the evolution of sand and dust erosion damage on aero-engine blades according to embodiments of this application.

[0031] The following is combined with Figure 1 and Figure 2 This application provides a detailed description of the method for predicting the evolution of sand and dust erosion damage on aero-engine blades.

[0032] Step S1: Construct a computational domain mesh using the complete blade's three-dimensional geometric model, and perform computational fluid dynamics (CFD) discrete-phase particle simulation using Fluent to obtain the current airflow field F of each stage of the initial complete blade. i Given the initial erosion rate R, each blade stage synchronously and sequentially performs erosion rate space mapping to obtain the current erosion result E. i , where i=0.

[0033] Specifically, the spatial mapping of erosion rate includes: unifying the coordinate units of the three-dimensional geometric model of the blade in stereolithography (STL) format with each erosion point containing spatial coordinates and corresponding initial erosion rates, and achieving spatial alignment through rotation transformation; mapping the initial erosion rate corresponding to each erosion point to the mesh vertices of the blade STL format three-dimensional geometric model to obtain a continuous erosion distribution field.

[0034] In one specific embodiment, spatial alignment specifically includes: Calculate the principal distribution angles of the mesh vertex coordinate set and the coordinate set of each erosion point in the Z-axis direction for the three-dimensional geometric model of the blade in STL format. Use the difference between the two principal distribution angles as the rotation angle, and determine the optimal alignment angle by taking the average nearest neighbor distance between the mesh vertex coordinate set and the coordinate set of each erosion point after rotation. The formula for calculating the optimal alignment angle is:

[0035] in: For the optimal alignment angle, This is the set of mesh vertex coordinates for the STL format 3D geometric model of the blade. For a single grid vertex in the set of grid vertex coordinates, (This is the 3D coordinate vector of the vertex.) This is the set of coordinates for each erosion point. R(θ) is the three-dimensional coordinate vector of a single erosion point in the set of coordinates of each erosion point, R(θ) is the rotation matrix about the Z-axis, and θ is the rotation angle. The initial erosion rate mapping for each erosion point is obtained using kernel density estimation and then smoothed using Laplace, specifically including: Calculate the area of ​​the dual element of each mesh vertex in the STL format three-dimensional geometric model of the blade. The area of ​​the dual element is equal to one-third of the sum of the areas of all triangular faces adjacent to the mesh vertex, and is used as a measure of the discrete area represented by the mesh vertex. Each erosion point is projected onto the nearest neighbor point on the surface of the blade's STL format 3D geometric model, and erosion points whose projection distance is greater than a preset multiple of the median projection distance are removed as outliers to obtain the projected erosion points. For each mesh vertex, a support domain is determined with the vertex as the center and a preset kernel bandwidth radius as the support radius. Within this support domain, a same-side normal condition is applied to filter mesh vertices and projected erosion points. That is, only erosion points and mesh vertices whose dot product of their unit outward normal vector and the unit outward normal vector of the mesh vertex is not less than a preset constant are retained. This is to avoid the pressure surface and suction surface interfering with each other when the blade is a thin-walled structure. The filtered erosion points and mesh vertices are obtained. The initial erosion rates corresponding to the filtered erosion points are weighted and summed according to the kernel function as the numerator, and the areas of the dual units of the filtered mesh vertices are weighted and summed according to the same kernel function as the denominator. The two are divided to obtain the erosion rate after mapping the mesh vertex. The erosion rate of mesh vertices that do not contain projected erosion points within the support domain is zero. The erosion rate of each mesh vertex after mapping is smoothed using Laplacian smoothing to eliminate noise; The formula for kernel density estimation is:

[0036] in, For grid vertices Mapped erosion rate For grid vertices The set of erosion points after center filtering, i.e., the set of grid vertices. The distance between them is no greater than the kernel bandwidth radius. The set of erosion points after projection that satisfy the same-side normal condition. For the erosion points after screening The initial erosion rate at the location, For the erosion points after screening The projected 3D coordinate vector For the erosion points after screening The position of the projected 3D coordinate vector and the mesh vertex The Euclidean distance between them For grid vertices The set of all grid vertices after filtering by the center, i.e., the set of all grid vertices. The distance between them is no greater than the kernel bandwidth radius. The set of mesh vertices that satisfy the condition of being normal to the same side. for The filtered grid vertices in the middle, For the filtered grid vertices The area of ​​the dual unit. For the filtered grid vertices With grid vertices The Euclidean distance between them This is a kernel function whose independent variable is the Euclidean distance between two points, and the independent variable is no greater than the kernel bandwidth radius. When the value is non-negative, and the independent variable is greater than the kernel bandwidth radius. The value is zero; the denominator is the weighted coverage area of ​​the kernel function on the blade surface on the same side of the grid vertex, i.e., the discrete form of the area integral of the kernel function along the blade surface on that side. It automatically decreases at the blade's leading edge, trailing edge, and tip because the support domain is truncated by the blade edge, thus achieving boundary correction; kernel bandwidth radius. The value is the larger of the particle size of the sand and dust and the average sampling interval of each erosion point, up to 20 times that larger value. The calculation formula is:

[0037] in, For the filtered grid vertices The area of ​​the dual unit, i.e., the area of ​​the filtered mesh vertices. The discrete area measure it represents To match the filtered grid vertices The set of adjacent triangular facets, for A single triangular facet in the middle, Triangular facet The area of ​​the blade; the sum of the areas of the dual elements of all filtered mesh vertices is always equal to the total surface area of ​​the blade, they do not overlap and completely cover each other; the same-side normal condition is:

[0038] in, , These are the grid vertices. and filtered grid vertices The unit outward normal vector at that location, The constant is used for the projected erosion points. The mean of the unit outward normal vectors of each vertex of the triangle on which its projection lies. Substituting the same condition, that is Erosion points and mesh vertices projected that do not meet this condition are not included. and .

[0039] The formula for the kernel function is:

[0040] in, For kernel function, Let be the Euclidean distance between the two points. For the kernel bandwidth radius; when hour, Take zero. The kernel function can also take the top-hat kernel. or Gaussian kernel .

[0041] The formula for Laplace smoothing is:

[0042] in, For grid vertices Erosion rate after smoothing treatment For grid vertices Mapped erosion rate For grid vertices a certain adjacent grid vertex Mapped erosion rate For grid vertices The set of adjacent grid vertices, It is a constant, with a value of 0.5.

[0043] A combined kernel density estimation and Laplace smoothing erosion rate mapping method effectively improves the continuity, conservation, and robustness of the erosion distribution field. Addressing the spatial mismatch between the discrete point cloud data derived from simulation and the mesh vertices of the blade's STL-format 3D geometric model, the statistical range is defined by the kernel bandwidth radius of each mesh vertex. A same-side normal condition is applied. The numerator is the kernel-weighted sum of the erosion rates of all erosion points within this support domain after same-side filtering, while the denominator is the kernel-weighted sum of the areas of the dual elements of all mesh vertices within this support domain after same-side filtering. The denominator represents the weighted coverage area of ​​the kernel function on the same-side blade surface. Its function is to convert the kernel-weighted mass erosion rate into the erosion rate per unit area. This denominator is automatically reduced at the blade's leading edge, trailing edge, and tip because the support domain is truncated by the blade edge, achieving boundary correction and preventing the erosion rate at edge vertices from being underestimated. The same-side normal condition ensures that, when the blade is a thin-walled structure, the pressure and suction surfaces are statistically independent and do not intersect. Since the kernel bandwidth radius is the preset physical length rather than the grid scale, when the grid is locally refined, the number of grid vertices participating in the summation in the denominator increases proportionally to the area, and the area of ​​the dual unit of each grid vertex decreases proportionally. The sum converges to a certain value of the area integral of the kernel function along the same side of the leaf surface. The numerator is the kernel weighted sum of physical quantities, so the obtained erosion rate does not change with the density of the grid.

[0044] Step S2: Based on the current erosion results E of each blade level i The wall thickness is calculated and the erosion depth is assessed, and the boundary contour of the erosion penetration area is extracted. A fast prediction method is used to determine the time period Δt required when the erosion volume of any first-stage blade reaches a preset damage threshold. i This step uses the current propagation time step as the starting point and simultaneously triggers a full-stage blade morphology update. This step automatically determines the geometric update cycle by using "determining that the erosion volume of any stage of the blade reaches a preset damage threshold" as the threshold, balancing prediction accuracy and computational efficiency, and avoiding invalid iterations or information loss caused by a fixed step size. Here, the preset damage threshold can be adjusted according to requirements.

[0045] Specifically, the wall thickness calculation and erosion depth assessment include: automatically calculating the wall thickness of each grid vertex based on the geometric characteristics of the thin-walled structure of the blade by anti-normal nearest neighbor matching; calculating the erosion depth by combining the continuous erosion distribution field; and extracting the boundary contour line of the erosion penetration area.

[0046] Wall thickness calculation specifically includes: Pairing point search: For a target vertex q, search for a pairing point Q within its K nearest neighbors that simultaneously satisfies the normal angle condition and the distance projection condition. Wall thickness assignment: If a pairing point Q is found, the Euclidean distance between the target vertex q and the pairing point Q is used as the wall thickness of the target vertex q, and thus the wall thickness field is obtained; if a pairing point Q is not found, the wall thickness field is obtained by interpolating using the known wall thickness in the neighborhood through inverse distance weighted interpolation. Wall thickness field smoothing: Laplace smoothing is applied to the wall thickness field; The condition for the anti-normal angle is: nq · nQ<-α Where: nq and nQ are the area-weighted normal vectors at the target vertex q and the paired point Q, respectively, and α is a constant with a value of 0.3; The distance projection condition is: |(q - Q) · nq| / ‖q - Q‖>β Where: β is a constant with a value of 0.3, ‖q - Q‖ is the Euclidean distance between the target vertex q and the paired point Q, and |(q-Q)·nq| is the projection length of the position difference vector between the target vertex q and the paired point Q onto the normal vector direction.

[0047] In one specific embodiment, the erosion depth assessment is based on the Finnie erosion model, combined with a continuous erosion distribution field, to obtain the erosion depth of each grid vertex; the erosion depth of each grid vertex is: h(x) = E(x) · T / ρ Where: h(x) is the erosion depth at grid vertex x, E(x) is the erosion rate after smoothing grid vertex x, T is the running time (s), and ρ is the blade material density; Extracting the boundary contour of the eroded area includes: When the erosion penetration ratio is ≥1.0, it is judged as complete erosion penetration; When the erosion penetration ratio is less than 1.0, it is judged as incomplete erosion penetration; At the interface between the fully eroded and partially eroded areas, the contour lines of the erosion boundary are extracted and smoothed using a contour line tracing algorithm to obtain the boundary contour lines of the eroded area. The erosion penetration ratio is: r(x) = h(x) / t(x) Where r(x) is the erosion penetration ratio, h(x) is the erosion depth at mesh vertex x, and t(x) is the original material thickness at mesh vertex x. That is, smooth tracing of the erosion boundary contour is achieved through contour boundary extraction.

[0048] In one specific embodiment, the rapid prediction method specifically includes: based on the current erosion results E of each level of the blades. iEach grid vertex and its corresponding erosion rate, combined with the blade material density ρ, according to h(x)=E i (x)·Δt / ρ, calculates the linear extrapolation relationship of the erosion volume at each grid vertex in each stage of the blade over time, where h(x) is the erosion depth at grid vertex x, Δt is the erosion time step, and E i (x) represents the current erosion rate at grid vertex x; the time interval required for the maximum erosion volume of any first-stage blade to first reach 2% of the current blade volume is determined as the current propagation time step Δt. i .

[0049] Specifically, in the rapid prediction method, the damage threshold that triggers geometric update is defined as the first time the cumulative erosion volume of any level of blade reaches 2% of the current blade volume. This is achieved within a time step Δt. i Within the grid, the erosion depth of each grid vertex is determined by h(x) = E. i (x)·Δt i / ρ increases linearly with time. Multiplying the erosion depth of each grid vertex by the area of ​​its dual element and summing the results over all grid vertices yields the blade's erosion depth at the propulsion time step Δt. i Cumulative erosion volume within:

[0050] Among them, V i (Δt i ) represents the propulsion time step Δt of this stage of blades. i The cumulative erosion volume within the grid, h(x) is the erosion depth at grid vertex x, a(x) is the area of ​​the dual cell at grid vertex x, and E i (x) represents the current erosion rate at grid vertex x, ρ is the blade material density, and the summation applies to all grid vertices of this stage of the blade. Let the cumulative erosion volume reach 2% of the current blade volume, i.e., V i (Δt i ) = 0.02 · V i Due to the cumulative erosion volume affecting Δt i It is strictly linear, and the current propagation time step Δt can be directly obtained by inverse solving. i :

[0051] Among them, V i Let E be the current blade volume, ρ be the blade material density, and E be the weight of the blade. i(x) is the current erosion rate at grid vertex x, and a(x) is the area of ​​the dual cell at grid vertex x. The above formula is used to calculate the erosion rate for each level of blade, and the minimum value is taken as the current advance time step, which is the time required for any level of blade to reach the 2% damage threshold first. The physical meaning of this method is as follows: select the most dangerous point with the largest ratio of the current erosion rate of the mesh vertex to the wall thickness of the mesh vertex, and use the time when the current erosion rate is extrapolated to 2% damage as the current advancement time step to ensure that each geometric update corresponds to a substantial morphological change, while avoiding excessively frequent geometric updates that lead to computational waste.

[0052] Step S3: Using the boundary contour of the erosion penetration area of ​​each blade stage and the current propulsion time step Δt i As input, each stage of the blade synchronously and sequentially performs blade geometric Boolean cutting, blade profile extraction, and parametric reconstruction, outputting the post-erosion geometric mesh M of each stage of the blade. i+1 .like Figure 2 and Figure 3 As shown, this step sequentially executes steps S310 to S320, resulting in the post-erosion geometric mesh M. i+1 This is used by Fluent to generate the computational mesh. This step only performs geometric reconstruction and does not involve solving the airflow field.

[0053] Step S310: Blade geometric Boolean cutting includes: converting the Initial Graphics Exchange Specification (IGES) format blade parametric model from Shell topology to Solid topology, constructing a cutting entity based on the boundary contour line of each erosion penetration area and performing Boolean subtraction operation to generate a damaged blade solid model.

[0054] Specifically, this includes: calculating the average surface normal direction of the nearest neighbor point for each boundary contour line of the erosion penetration area. Projecting the contour onto a vertical plane constructs a closed curve and generates a surface, along... The blade solid model is generated by bidirectional stretching a preset distance; Boolean subtraction operations are performed on the blade solid model and each cut solid model in sequence to output the damaged blade solid model.

[0055] Step S320: Blade profile extraction and parametric reconstruction includes: sectioning and contouring the damaged blade according to the percentage of blade height, automatically dividing the pressure surface and suction surface, and outputting a standard format profile file for the next mesh generation.

[0056] Specifically, this includes: fitting B-spline curves to the streamlines of each blade height in the damaged blade solid model according to the percentage of blade height and generating a surface of revolution; performing a sectioning operation on the damaged blade solid model to obtain a two-dimensional contour point set of each blade height streamline; performing connectivity analysis on the two-dimensional contour point set based on the K nearest neighbors, constructing a closed curve through Eulerian circuit search, and performing adaptive point density adjustment; and calculating the three-dimensional curvature κ of each point on the closed curve. i The maximum curvature point in the extreme region of the Z coordinate is searched as the leading and trailing edge split point, dividing the closed curve into two open curves, which are automatically labeled as pressure surface profiles and suction surface profiles according to geometric criteria. After the pressure surface profiles and suction surface profiles are processed by Gaussian smoothing or Savitzky-Golay filtering, all blade height streamlines are sorted according to the percentage of blade height and assembled and output in geomTurbo format. Here, the closed curve is the closed curve of the K-Nearest Neighbors (KNN) algorithm graph. Three-dimensional curvature κ i The formula is: κ i = 2‖(P i+1 - P i ) × (P i - P i-1 )‖ / (‖P i+1 - P i ‖ · ‖P i - P i-1 ‖ · ‖P i+1 - P i-1 ‖) Among them, P i Let P be the three-dimensional coordinates of the i-th point on the closed curve. i+1 Let P be the three-dimensional coordinates of the (i+1)th point on the closed curve. i-1 Let be the three-dimensional coordinates of the (i-1)th point on the closed curve.

[0057] Step S4: Calculate the current airflow field F for each stage of the blades. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 This step aims to utilize the convergent airflow field from the previous simulation to provide a high-quality initial field for the next simulation, thereby significantly reducing the number of iterations required for convergence.

[0058] Specifically, such as Figure 4 As shown, step S4 includes: completing the post-erosion geometry mesh M in Fluent. i+1Importing, defining the region, and setting boundary conditions; calling Fluent's "Interpolate Data" function to export the velocity field, pressure field, temperature field, and turbulent flow rate after the convergence of the CFD discrete phase particle simulation as nodal data, and mapping it to the post-erosion geometric mesh M by spatial interpolation. i+1 This serves as the initial airflow field for the next Fluent simulation, enabling rapid convergence; the geometric mesh M after erosion... i+1 The interpolated airflow field was used as the initial airflow field for Fluent simulation to obtain the new airflow field F for each stage of the blades. i+1 and new erosion results E i+1 Here, spatial interpolation uses inverse distance weighting by default.

[0059] The core principle of this airflow field initiation method is that, because this invention ensures that the cumulative erosion volume change corresponding to each geometric update of the blades at each stage is 2%, the change in blade morphology between adjacent iterations is small, and the airflow field structure is highly similar. Therefore, the convergent airflow field of the previous step is very close to the final convergent field of the next step. Using it as the initial field can significantly reduce the number of iterations required from the initial residual to convergence. Engineering tests show that, compared to zero initialization, this method can reduce the number of iterations required for a single simulation convergence by more than 50%.

[0060] Step S5: Determine whether the total erosion time of the blade has reached the preset value, wherein: If the target is not reached, let i = i + 1 and return to steps S2-S5; If the condition is met, the entire damage evolution result will be output.

[0061] Specifically, the total erosion duration is the total simulation time, and the preset can be set as needed.

[0062] If the total erosion time of the leaf does not reach the preset value, let i = i + 1, return to step S2 to continue the loop iteration, and use the new erosion result E i+1 Determine the next propagation time step, and so on. If the total erosion time of the blade reaches the preset value, summarize the damage evolution results of each round (blade morphology at each time point, airflow field at each stage, and profile file), and output the full process results of blade sand and dust erosion damage evolution, which will be used for subsequent remaining life assessment, maintenance strategy formulation, and anti-erosion optimization design.

[0063] In another specific embodiment, the method for predicting the evolution of sand and dust erosion damage on aero-engine blades further includes: in step S5, if the total erosion duration of the blade reaches a preset damage threshold and the corresponding erosion morphology result is obtained, then the irrelevance of the current propulsion time step is verified, specifically including: Based on a preset damage threshold, steps S2-S4 are repeated with a smaller damage threshold to obtain the erosion morphology results corresponding to the smaller damage threshold, where: If the error in the compressor aerodynamic performance calculated by the erosion morphology result corresponding to the preset damage threshold and the erosion morphology result corresponding to the smaller damage threshold is less than 5%, then the erosion morphology result corresponding to the preset damage threshold is convergent and reliable, the full damage evolution result is output and the process is terminated. If the error in the compressor aerodynamic performance calculated from the erosion morphology result corresponding to the preset damage threshold is greater than or equal to 5% compared with the erosion morphology result corresponding to a smaller damage threshold, then the damage threshold is further reduced and steps S2-S4 are repeated until the erosion morphology result corresponding to the preset damage threshold converges reliably. The entire damage evolution result is then output and the process terminates.

[0064] For example, based on the time required for the maximum erosion volume of any stage blade to first reach 2% of the current blade volume, the current propulsion time step Δt is determined. i The results show that when verifying the independence of propulsion time step, the threshold for smaller damage can be set when the maximum erosion volume of any stage blade first reaches 1.5% of the current blade volume, and the above operation can be performed.

[0065] This application also provides a system for performing the above-described method for predicting the evolution of sand and dust erosion damage on aero-engine blades, comprising: The initial simulation module is used to construct a computational domain mesh using a complete three-dimensional geometric model of the blade. Fluent is then used for CFD discrete-phase particle simulation to obtain the current airflow field Fi and initial erosion rate R of each stage of the initial complete blade. Subsequently, each stage of the blade synchronously and sequentially performs erosion rate space mapping to obtain the current erosion result E. i , where i=0; A rapid prediction module is used to predict erosion results based on the current erosion status of each blade level. i The wall thickness is calculated and the erosion depth is assessed, and the boundary contour of the erosion penetration area is extracted. A fast prediction method is used to determine the time period Δt required when the erosion volume of any first-stage blade reaches a preset damage threshold. i This serves as the current propulsion time step and simultaneously triggers a full-stage blade morphology update. The mesh reconstruction module is used to define the boundary contours of the erosion penetration areas of each blade stage and the current propulsion time step Δt. i As input, each stage of the blade synchronously and sequentially performs blade geometric Boolean cutting, blade profile extraction, and parametric reconstruction, outputting the post-erosion geometric mesh M of each stage of the blade. i+1 ; The airflow field interpolation module is used to interpolate the current airflow field F of each stage of the blades. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ; The iterative control module is used to determine whether the total erosion time of the blade has reached a preset value, wherein: If the target is not reached, let i = i + 1 and return to steps S2-S5; If the condition is met, the entire damage evolution result will be output.

[0066] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. The scope of patent protection of the present invention shall be determined by the claims. Similarly, any equivalent structural changes made based on the description and drawings of the present invention shall also be included within the scope of protection of the present invention.

Claims

1. A method for predicting the evolution of sand and dust erosion damage on aero-engine blades, characterized in that, Includes the following steps: Step S1: Construct a computational domain mesh using the complete blade's three-dimensional geometric model, and perform CFD discrete-phase particle simulation using Fluent to obtain the current airflow field F of each stage of the initial complete blade. i Given the initial erosion rate R, each blade stage synchronously and sequentially performs erosion rate space mapping to obtain the current erosion result E. i , where i=0; Step S2: Based on the current erosion results E of each blade level i The wall thickness is calculated and the erosion depth is assessed, and the boundary contour of the erosion penetration area is extracted. A fast prediction method is used to determine the time period Δt required when the erosion volume of any first-stage blade reaches a preset damage threshold. i This serves as the current propulsion time step and simultaneously triggers a full-stage blade morphology update. Step S3: Using the boundary contour of the erosion penetration area of ​​each blade stage and the current propulsion time step Δt i As input, each stage of the blade synchronously and sequentially performs blade geometric Boolean cutting, blade profile extraction, and parametric reconstruction, outputting the post-erosion geometric mesh M of each stage of the blade. i+1 ; Step S4: Calculate the current airflow field F for each stage of the blades. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ; Step S5: Determine whether the total erosion time of the blade has reached the preset value, wherein: If the target is not reached, let i = i + 1 and return to steps S2-S5; If the condition is met, the entire damage evolution result will be output.

2. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 1, characterized in that, In step S1, the erosion rate spatial mapping specifically includes: The STL format three-dimensional geometric model of the blade is aligned with the coordinate units of each erosion point containing spatial coordinates and corresponding initial erosion rates, and spatial alignment is achieved through rotation transformation. The initial erosion rate corresponding to each erosion point is mapped to the mesh vertices of the blade's STL format three-dimensional geometric model to obtain a continuous erosion distribution field.

3. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 2, characterized in that, The spatial alignment specifically includes: Calculate the principal distribution angles of the mesh vertex coordinate set and the coordinate set of each erosion point in the Z-axis direction for the three-dimensional geometric model of the blade in STL format. Use the difference between the two principal distribution angles as the rotation angle, and determine the optimal alignment angle by taking the average nearest neighbor distance between the mesh vertex coordinate set and the coordinate set of each erosion point after rotation. The formula for calculating the optimal alignment angle is: in: For the optimal alignment angle, This is the set of mesh vertex coordinates for the STL format 3D geometric model of the blade. For a single grid vertex in the set of grid vertex coordinates, (This is the 3D coordinate vector of the vertex.) This is the set of coordinates for each erosion point. R(θ) is the three-dimensional coordinate vector of a single erosion point in the set of coordinates of each erosion point, R(θ) is the rotation matrix about the Z-axis, and θ is the rotation angle. The initial erosion rate mapping corresponding to each erosion point adopts the kernel density estimation method and is smoothed by Laplacian. Specifically, each erosion point is projected to the nearest neighbor point on the surface of the blade STL format three-dimensional geometric model, and erosion points with a projection distance greater than a preset multiple of the projection distance are removed as outliers to obtain the projected erosion points. For each mesh vertex, a support domain is determined with the mesh vertex as the center and the preset kernel bandwidth radius as the support radius. Within this support domain, the same-side normal condition is applied to filter the mesh vertices and the projected erosion points, resulting in each filtered erosion point and each mesh vertex. The initial erosion rates corresponding to each erosion point after screening are weighted and summed according to the kernel function as the numerator, and the areas of the dual units of each grid vertex after screening are weighted and summed according to the same kernel function as the denominator, to obtain the erosion rate after mapping the grid vertex. The erosion rate of each mapped mesh vertex is smoothed using Laplacian smoothing to eliminate noise; the formula for the kernel density estimation method is: in, For grid vertices Mapped erosion rate For grid vertices The collection of erosion points after center screening. For the erosion points after screening The initial erosion rate at the location, For the erosion points after screening The projected 3D coordinate vector For the erosion points after screening The position of the projected 3D coordinate vector and the mesh vertex The Euclidean distance between them For grid vertices The set of all grid vertices after filtering by the center. for The filtered grid vertices For the filtered grid vertices The area of ​​the dual unit. For the filtered grid vertices With grid vertices The Euclidean distance between them For kernel functions; The The calculation formula is: in, For the filtered grid vertices The area of ​​the dual unit. To match the filtered grid vertices The set of adjacent triangular facets, for A single triangular facet in the middle, Triangular facet The area; The formula for Laplace smoothing is: in, For grid vertices Erosion rate after smoothing treatment For grid vertices Mapped erosion rate For grid vertices a certain adjacent grid vertex Mapped erosion rate Let x be the set of adjacent grid vertices. It is a constant.

4. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 3, characterized in that, In step S2, the calculation of wall thickness and assessment of erosion depth, and the extraction of the boundary contour line of the erosion penetration area, specifically include: Based on the geometric characteristics of the thin-walled blade structure, the wall thickness of each grid vertex is automatically calculated by anti-normal nearest neighbor matching; the erosion depth is calculated by combining the continuous erosion distribution field, and the boundary contour line of the erosion penetration area is extracted. The wall thickness calculation specifically includes: Pairing point search: For a target vertex q, search for a pairing point Q within its K nearest neighbors that simultaneously satisfies the inverse normal angle condition and the distance projection condition. Wall thickness assignment: If a pairing point Q is found, the Euclidean distance between the target vertex q and the pairing point Q is used as the wall thickness of the target vertex q, and thus the wall thickness field is obtained; if a pairing point Q is not found, the wall thickness field is obtained by interpolating using the known wall thickness in the neighborhood through inverse distance weighted interpolation. Wall thickness field smoothing: Laplace smoothing is applied to the wall thickness field; The condition for the anti-normal included angle is: nq · nQ<-α Where: nq and nQ are the area-weighted normal vectors at the target vertex q and the paired point Q, respectively, and α is a constant; The distance projection condition is: |( q - Q) · nq| / ‖q - Q‖>β Where: β is a constant, ‖q - Q‖ is the Euclidean distance between the target vertex q and the paired point Q, and |( q - Q) · nq| is the projection length of the position difference vector between the target vertex q and the paired point Q in the direction of the normal vector.

5. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 4, characterized in that, The erosion depth assessment is based on the Finnie erosion model, combined with a continuous erosion distribution field, to obtain the erosion depth of each grid vertex; The erosion depth of each grid vertex is: h(x) = E(x) · T / ρ Where: h(x) is the erosion depth at grid vertex x, E(x) is the erosion rate after smoothing grid vertex x, T is the running time, and ρ is the blade material density; The extracted boundary contour line of the erosion penetration area includes: When the erosion penetration ratio is ≥1.0, it is judged as complete erosion penetration; When the erosion penetration ratio is less than 1.0, it is judged as incomplete erosion penetration; At the interface between the fully eroded and partially eroded areas, the contour lines of the erosion boundary are extracted and smoothed using a contour line tracing algorithm to obtain the boundary contour lines of the eroded area. The erosion penetration ratio is: r(x) = h(x) / t(x) Where r(x) is the erosion penetration ratio, h(x) is the erosion depth at mesh vertex x, and t(x) is the original material thickness at mesh vertex x.

6. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 5, characterized in that, The rapid prediction method specifically includes: Based on the current erosion results of each blade level E i Each grid vertex and its corresponding erosion rate, combined with the blade material density ρ, according to h(x)=E i (x)·Δt / ρ, calculates the linear extrapolation relationship of the erosion volume at each grid vertex in each stage of the blade over time, where h(x) is the erosion depth at grid vertex x, Δt is the erosion time step, and E i (x) represents the current erosion rate at grid vertex x; The time interval required for the maximum erosion volume of any stage blade to first reach 2% of the current blade volume is defined as the current propulsion time step Δt. i .

7. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 1, characterized in that, The blade geometric Boolean cutting specifically includes: For each boundary contour line of the erosion penetration area, calculate the average surface normal direction of the nearest neighbor point. Projecting the contour onto a vertical plane constructs a closed curve and generates a surface, along... Bidirectional stretching by a preset distance generates a cut entity; Perform Boolean subtraction operations on the blade solid model and each cut solid model in sequence, and output the damaged blade solid model.

8. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 7, characterized in that, The blade profile extraction and parameter reconstruction specifically include: Based on the percentage of leaf height, each streamline of the damaged blade solid model is fitted as a B-spline curve and a surface of revolution is generated. The damaged blade solid model is then sectioned to obtain a two-dimensional contour point set of each streamline. Connectivity analysis is performed on the two-dimensional contour point set based on the K nearest neighbors. A closed curve is constructed by searching for Eulerian circuits, and adaptive point density adjustment is performed. Calculate the three-dimensional curvature κ at each point on the closed curve. i The maximum point of three-dimensional curvature is searched in the extreme region of the Z coordinate as the dividing point of the front and rear edges, and the closed curve is divided into two open curves, which are automatically labeled as pressure surface profile and suction surface profile according to geometric criteria. After the pressure surface profile and suction surface profile are processed by Gaussian smoothing or Savitzky-Golay filtering, all blade height streamlines are sorted by blade height percentage and assembled and output in geomTurbo format; The three-dimensional curvature κ i The formula is: κ i = 2‖(P i+1 - P i ) × (P i - P i-1 )‖ / (‖P i+1 - P i ‖ · ‖P i - P i-1 ‖ · ‖P i+1 - P i-1 ‖) Among them, P i Let P be the three-dimensional coordinates of the i-th point on the closed curve. i+1 Let P be the three-dimensional coordinates of the (i+1)th point on the closed curve. i-1 Let be the three-dimensional coordinates of the (i-1)th point on the closed curve.

9. The method for predicting the evolution of sand and dust erosion damage on aero-engine blades according to claim 1, characterized in that, The current airflow field F of each stage of the blades is described. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ,include: Complete the post-erosion geometry M in Fluent i+1 Importing, defining regions, and setting boundary conditions; The velocity, pressure, temperature, and turbulent flow rates after convergence of the CFD discrete-phase particle simulation are exported as nodal data and spatially interpolated onto the post-erosion geometric mesh M. i+1 This serves as the initial airflow field for the next CFD discrete phase particle simulation. Post-erosion geometric mesh M i+1 The interpolated airflow field was used as the initial airflow field for CFD discrete phase particle simulation to obtain the new airflow field F for each stage of the blades. i+1 and new erosion results E i+1 .

10. A system for performing the method for predicting the evolution of sand and dust erosion damage on aero-engine blades as described in any one of claims 1-9, characterized in that, include: The initial simulation module is used to construct a computational domain mesh using a complete three-dimensional geometric model of the blade. Fluent is then used for CFD discrete-phase particle simulation to obtain the current airflow field Fi and initial erosion rate R of each stage of the initial complete blade. Subsequently, each stage of the blade synchronously and sequentially performs erosion rate space mapping to obtain the current erosion result E. i , where i=0; A rapid prediction module is used to predict erosion results based on the current erosion status of each blade level. i The wall thickness is calculated and the erosion depth is assessed, and the boundary contour of the erosion penetration area is extracted. A fast prediction method is used to determine the time period Δt required when the erosion volume of any first-stage blade reaches a preset damage threshold. i This serves as the current propulsion time step and simultaneously triggers a full-stage blade morphology update. The mesh reconstruction module is used to define the boundary contours of the erosion penetration areas of each blade stage and the current propulsion time step Δt. i As input, each stage of the blade synchronously and sequentially performs blade geometric Boolean cutting, blade profile extraction, and parametric reconstruction, outputting the post-erosion geometric mesh M of each stage of the blade. i+1 ; The airflow field interpolation module is used to interpolate the current airflow field F of each stage of the blades. i Mapped to the post-erosion geometric mesh M using Fluent data interpolation. i+1 As the initial airflow field for the next CFD discrete phase particle simulation, the new airflow field F for each stage of the blades is obtained. i+1 and new erosion results E i+1 ; The iterative control module is used to determine whether the total erosion time of the blade has reached a preset value, wherein: If the target is not reached, let i = i + 1 and return to steps S2-S5; If the condition is met, the entire damage evolution result will be output.

Images