Method and system for predicting crack propagation of a water turbine runner under transient hydraulic excitation

Abstract

The application discloses a kind of crack propagation prediction method and system of water turbine runner under transient hydraulic excitation, belong to the technical field of water conservancy and hydropower engineering, to solve the technical problem that existing prediction method cannot simulate crack propagation under transient hydraulic excitation.The application is by constructing finite element grid and level set field, encryption grid, identify enhanced node and solve runner time sequence displacement field, then calculate the stress intensity factor and energy release rate of crack front node, and calculate crack propagation direction and rate based on this, thereby update and correct level set field, to obtain the crack after expansion, to realize the accurate simulation of crack propagation under transient hydraulic excitation, provide reliable basis for runner crack monitoring, expansion trend prediction and fracture life assessment.Further, the application can track the crack propagation process from the current small condition to the penetration by multiple iterations, and accurately predict the fracture life of the runner according to the propagation law.

Description

Technical Field

[0001] This application relates to the field of water conservancy and hydropower engineering technology, and in particular to a method and system for predicting crack propagation in a turbine runner under transient hydraulic excitation. Background Technology

[0002] With the improvement of turbine efficiency and speed, and the thinning of runner blades, crack initiation and propagation caused by metal structural fatigue under long-term pulsating hydraulic excitation conditions are the main failure modes of turbine runners. In actual blade manufacturing, uncertainties exist in materials and processes, and the long intervals between periodic disassembly and inspection are costly. This makes it impossible for turbine maintenance personnel to predict and detect runner cracks in a timely manner, and to determine whether known micro-cracks will continue to propagate under operating conditions, or the operating time required for micro-cracks to evolve into through-cracks. Once a runner blade fractures, the turbine unit efficiency decreases, and shaft oscillation and resonance may occur due to mass asymmetry. In severe cases, the detached blade fragment, under high-speed rotation, may enter the volute, draft tube, and other components, causing chain reactions and resulting in serious power plant accidents and loss of life and property. Therefore, predicting crack propagation paths and fracture lives using numerical simulation methods for known cracks is of great significance for achieving predictive maintenance throughout the entire turbine lifecycle.

[0003] Numerical simulation of turbine runner blade fracture failure generally employs the pre-crack method. This method requires creating mesh elements for the crack surface and crack lead in a finite element mesh. Its drawback is the need for real-time updates of the cracked mesh during crack propagation. Therefore, when dealing with cracks in complex topologies such as turbine runners, it suffers from difficulties in mesh preprocessing and is time-consuming, making long-term safety assessments impossible.

[0004] This led to the development of the Extended Finite Element Method (EFEM). However, existing EFEM methods are rarely used in the field of turbine fracture fault assessment, and can only be used for simple evaluations. Due to the complex operating conditions of turbines in actual operation, existing EFEM methods cannot be directly applied to turbine fracture fault assessment. For example, patent document CN119167595A discloses a method for simulating turbine runner blade defects based on CFD and FEM. It uses computational fluid dynamics and dynamic mesh technology to obtain the distribution of multiple physical quantities in the internal flow field under faulty runner blade conditions, and obtains the natural frequency of the faulty runner through modal calculations. Combining the natural frequency and the hydraulic excitation frequency, it compares and judges the probability of resonance in the faulty runner. However, the patent solution can only analyze the impact of known runner wear and fracture faults on structural safety through numerical simulation. It cannot track the crack propagation process from micro-cracks to full penetration. The defect morphology needs to be preset during the calculation, which makes this solution unusable in predicting crack propagation technology, let alone calculating the crack propagation rate under transient hydraulic excitation and predicting fracture life. Due to the above-mentioned shortcomings of the existing technology, it is necessary to improve it. Summary of the Invention

[0005] To address the aforementioned problems, or one of them, the present invention aims to provide a method for predicting crack propagation in turbine runners under transient hydraulic excitation. This method constructs a finite element mesh and a horizontal concentrated field, refines the mesh, identifies reinforced nodes, and solves the time-series displacement field of the runner. It then calculates the stress intensity factor and energy release rate of the crack front nodes, and uses this as a basis to calculate the crack propagation direction and rate. The horizontal concentrated field is then updated and corrected to obtain the propagated crack, thereby achieving accurate simulation of crack propagation under transient hydraulic excitation and providing a reliable basis for runner crack monitoring, propagation trend prediction, and fracture life assessment.

[0006] To address the aforementioned problems or one of them, the second objective of this invention is to provide a crack propagation prediction system for turbine runners under transient hydraulic excitation. This system constructs a finite element mesh and a horizontal concentrated field, refines the mesh, identifies reinforced nodes, and solves the time-series displacement field of the runner. It then calculates the stress intensity factor and energy release rate of the crack front nodes, and uses these as a basis to calculate the crack propagation direction and rate. Based on this, the horizontal concentrated field is updated and corrected to obtain the propagated crack, thereby achieving accurate simulation of crack propagation under transient hydraulic excitation and providing a reliable basis for runner crack monitoring, propagation trend prediction, and fracture life assessment.

[0007] To achieve one of the above objectives, the first technical solution of the present invention is as follows: A method for predicting crack propagation in a turbine runner under transient hydraulic excitation includes the following steps: Step 1: Draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, including the turbine blades; draw the horizontal concentrated field based on the cracks on the turbine blades, including the normal horizontal concentrated field and the tangential horizontal concentrated field; Step 2: Determine the crack front based on the horizontal field, refine the mesh near the crack front to form the target mesh; screen the nodes of the target mesh, identify the reinforcing nodes, determine the type of reinforcing element where the reinforcing nodes are located, and solve the time-series displacement field of the turbine runner under transient hydraulic excitation; Step 3: Based on the time-series displacement field of the turbine runner, calculate the stress intensity factor and energy release rate of the crack front node, and use this as a basis to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation; Step 4: Based on the crack propagation direction and crack propagation rate obtained in Step 3, update and correct the horizontal field to obtain the propagated crack.

[0008] As a preferred technical measure, it also includes step five: predicting the fracture life of the impeller; Step four further includes: determining whether a through crack has formed based on the expanded crack; if a through crack is formed, proceed to step five; otherwise, return to step two for iteration.

[0009] As a preferred technical measure: In step two, the mesh near the crack tip is a tetrahedral mesh; the step of refining the mesh near the crack tip to form the target mesh is as follows: Step A1: Determine the target mesh size based on the crack length, and then determine the number of mesh refinements by combining the original mesh size near the crack leading edge with the target mesh size; Step A2: Use the distance control index field to filter out grid cells that meet the encryption requirements; Step A3: Process the mesh cells that meet the encryption requirements, set the midpoints of each of their edges as newly added mesh points, and divide the original tetrahedral mesh cells into four new mesh cells. Step A4: Determine whether the encryption count has been reached. If yes, complete the encryption; otherwise, return to step A2.

[0010] As a preferred technical measure: In step two, the method for filtering the nodes of the target mesh and identifying the enhanced nodes is as follows: Find the mesh elements in the target mesh that cross the zero isosurface at the crack front, and use the nodes of these mesh elements as candidate nodes; All candidate nodes are traversed to filter out crack-enhanced nodes and crack tip-enhanced nodes; the crack-enhanced nodes are used to characterize the displacement discontinuity of the crack surface; the crack tip-enhanced nodes are used to characterize the stress singularity of the crack tip. In step two, the method for determining the type of enhancement unit where the enhancement node resides is as follows: A mesh element containing at least one crack reinforcement node is defined as a crack reinforcement element; a mesh element containing at least one crack tip reinforcement node is defined as a crack tip reinforcement element; and a mesh element containing both is defined as a hybrid reinforcement element.

[0011] As a preferred technical measure: the method for solving the time-series displacement field of the runner under transient hydraulic excitation in step two is as follows: Based on the horizontal field, the reinforcement nodes near the crack tip are selected and the reinforcement elements in which they are located are determined; Based on the aforementioned enhancement unit, the unit stiffness matrix and mass matrix are calculated, and then assembled into the overall stiffness matrix and mass matrix of the impeller. Based on fluid dynamics simulation, data on pulsating fluid pressure loads under transient hydraulic excitation are obtained. Based on the overall stiffness matrix and mass matrix of the runner, as well as the pulsating fluid pressure load data, the acoustic-structure interaction dynamic equilibrium equation of the runner is solved to obtain the time-series displacement field of the runner.

[0012] As a preferred technical measure: In step three, the step of calculating the stress intensity factor and energy release rate of the crack front node based on the time-series displacement field of the rotor is as follows: Step B1: Obtain the location of the crack tip, take several nodes that are evenly distributed near the crack tip as the original calculation nodes, and establish a local coordinate system of the crack tip that matches the geometric direction of the crack tip. Step B2: Obtain the displacement field at each moment based on the time-series displacement field of the rotor, and extract the displacement results on both sides of the crack surface at each moment; Step B3: Transform the displacement results obtained in step B2 into the local coordinate system of the crack front, and then project them onto the original calculation nodes to form projected calculation nodes; Step B4: Based on the projection calculation nodes obtained in Step B3, draw a scatter plot and obtain the stress intensity factor of the crack front node at each time point through the scatter plot. Step B5: Based on the stress intensity factor obtained in step B4, the energy release rate of the crack front node at each time step is calculated using the Evan formula.

[0013] As a preferred technical measure: In step three, the step of calculating the crack propagation direction under transient hydraulic excitation is as follows: Step C1: Based on the stress intensity factor of the crack leading edge node at each moment under transient hydraulic excitation, the maximum circumferential stress criterion is used to calculate the deflection angle of each crack leading edge node at each moment; the deflection angle is located in the plane formed by the crack leading edge tangent and the crack surface normal. Step C2: Collect the deflection angle data of all crack leading edge nodes at each time step; Step C3: Calculate the average deflection angle data collected in step C2 over time. The average deflection angle obtained is the crack propagation direction. In step three, the method for calculating the crack propagation rate under transient hydraulic excitation is as follows: Calculate the corresponding equivalent stress intensity factor using the energy release rate of the crack front node at each time point; The equivalent stress intensity factor at each time point is statistically analyzed across the entire range, and the amplitude of the equivalent stress intensity factor corresponding to all load cycles is extracted. Obtain the crack propagation threshold value of the rotating wheel; The amplitude of the equivalent stress intensity factor corresponding to each load cycle is judged. When it is greater than the crack propagation threshold, the crack propagation rate is calculated by the Paris formula; otherwise, the assumed value is forcibly set as the crack propagation rate.

[0014] As a preferred technical measure: Step four: Based on the crack propagation direction and crack propagation rate obtained in step three, update and correct the horizontal field to obtain the propagated crack. Step D1: Based on the tangential level set, divide all meshes into cracked and uncracked regions; For the already cracked area, the propagation rate of the normal horizontal set is fixed to zero, so that the existing crack surface remains unchanged. For the uncracked zone, the crack propagation rate is decomposed into normal propagation rate and tangential propagation rate based on the crack propagation direction obtained in step three. Step D2: Based on the normal propagation rate and the tangential propagation rate, update the horizontal field after crack propagation in the uncracked zone to characterize the updated crack surface and crack front topology, so as to obtain the expanded crack. Step D3: Correct the horizontal field obtained after crack propagation in step D2, ensuring that the absolute value of the new horizontal gradient field is always equal to 1, and that the normal horizontal field and the tangential horizontal field on the crack surface are orthogonal.

[0015] As a preferred technical measure, step five, the method for predicting the fracture life of the impeller, includes: Extract the crack length, maximum equivalent stress intensity factor, and maximum equivalent stress intensity factor amplitude in each iteration step; Plot the first relationship curve between the maximum equivalent stress intensity factor and the crack length. Based on this first relationship curve and the fracture toughness of the runner material, determine whether the runner blade will undergo brittle fracture. If there is a possibility of brittle fracture, determine the critical crack length for brittle fracture and provide monitoring prompts. A second relationship curve is plotted between the amplitude of the maximum equivalent stress intensity factor and the crack length. Based on this second relationship curve and the fatigue crack propagation threshold value of the impeller material, it is determined whether fatigue fracture will occur due to fatigue crack propagation. When fatigue fracture is possible, the critical crack length at which propagation begins is determined, and the time from the critical crack length at which propagation begins to the formation of a through crack is calculated. The calculation result is the fracture life of the impeller.

[0016] To achieve one of the above objectives, the second technical solution of the present invention is as follows: A crack propagation prediction system for a turbine runner under transient hydraulic excitation includes the following modules: The mesh generation module is used to draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, the rotating components including the runner blades; and to draw a horizontal concentrated field based on the cracks on the runner blades, the horizontal concentrated field including a normal horizontal concentrated field and a tangential horizontal concentrated field. The temporal displacement field generation module is used to determine the crack front based on the horizontal concentrated field, refine the mesh near the crack front to form the target mesh, filter the nodes of the target mesh, identify the reinforcing nodes, determine the type of reinforcing element where the reinforcing nodes are located, and solve the temporal displacement field of the runner under transient hydraulic excitation. The crack propagation direction and rate generation module is used to calculate the stress intensity factor and energy release rate of the crack front node based on the time-series displacement field of the runner, and to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation based on these. Crack propagation module: Used to update and correct the horizontal field based on the crack propagation direction and crack propagation rate to obtain the propagated crack.

[0017] Compared with existing technical solutions, the present invention has the following beneficial effects: This invention constructs a finite element mesh and a horizontal concentrated field, and accurately characterizes the crack geometry through the horizontal concentrated field, eliminating the need for real-time updates of the cracked mesh during crack propagation. This solves the technical problems of difficult and time-consuming mesh preprocessing in existing technologies. By refining the mesh, the horizontal concentrated field's characterization of the crack geometry can be improved, enhancing the accuracy of stress intensity factor calculation and ensuring accurate crack propagation simulation. Furthermore, it avoids full-mesh refinement, balancing accuracy and efficiency. By identifying reinforced nodes and solving the time-series displacement field of the turbine runner, it accurately obtains the turbine's force and displacement data under transient hydraulic excitation, avoiding the limitations of traditional extended finite element methods due to load limitations. This method addresses calculation errors caused by inaccurate load acquisition, improving the fundamental accuracy of crack mechanical parameter calculations. It calculates the stress intensity factor and energy release rate at the crack front node based on the turbine's time-series displacement field, thereby determining the crack propagation direction and rate. This enables precise capture of crack propagation patterns under transient hydraulic excitation, updating and correcting the horizontal field to obtain the propagated crack. This allows for accurate simulation of crack propagation under transient hydraulic excitation, providing reliable geometric input for multi-round iterative prediction, aiding in the prediction of turbine crack propagation trends and the formulation of operation and maintenance plans, and providing a reliable basis for turbine crack fracture life assessment. Furthermore, through multiple iterations, this invention can track the expansion process of cracks from their current minor state to full penetration, and accurately predict the fracture life of the turbine runner based on its expansion pattern. This helps hydropower stations plan maintenance and turbine runner replacement in advance, effectively avoiding shutdown accidents caused by sudden turbine runner fracture, ensuring the safe and stable operation of the turbine, reducing operation and maintenance costs, and improving the overall operating efficiency of the hydropower station. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the specific embodiments of this application or the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0019] Figure 1 This is a flowchart of the method in Embodiment 2 of the present invention; Figure 2 This is a diagram showing the mesh refinement result of the final crack propagation simulation of the crack surface and the cracked impeller blade in Embodiment 3 of the present invention. Figure 3 This is a graph showing the variation of the maximum equivalent stress intensity factor with crack length in Embodiment 3 of the present invention. Figure 4 This is a graph showing the variation of the maximum equivalent stress intensity factor with crack length in Embodiment 3 of the present invention. Detailed Implementation

[0020] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings.

[0021] Example 1 This embodiment describes a method for predicting crack propagation in a turbine runner under transient hydraulic excitation, including the following steps: Step 1: Draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, including the turbine blades; draw the horizontal concentrated field based on the cracks on the turbine blades, including the normal horizontal concentrated field and the tangential horizontal concentrated field; Step 2: Determine the crack front based on the horizontal field, refine the mesh near the crack front to form the target mesh; screen the nodes of the target mesh, identify the reinforcing nodes, determine the type of reinforcing element where the reinforcing nodes are located, and solve the time-series displacement field of the turbine runner under transient hydraulic excitation; Step 3: Based on the time-series displacement field of the turbine runner, calculate the stress intensity factor and energy release rate of the crack front node, and use this as a basis to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation; Step 4: Based on the crack propagation direction and crack propagation rate obtained in Step 3, update and correct the horizontal field to obtain the propagated crack.

[0022] This embodiment constructs a finite element mesh and a horizontal concentrated field, refines the mesh, identifies enhanced nodes, and solves the time-series displacement field of the runner. Then, it calculates the stress intensity factor and energy release rate of the crack front node, and uses this as a basis to calculate the crack propagation direction and rate. The horizontal concentrated field is then updated and corrected to obtain the propagated crack, thereby achieving accurate simulation of crack propagation under transient hydraulic excitation. This provides a reliable basis for runner crack monitoring, propagation trend prediction, and fracture life assessment.

[0023] Example 2 like Figure 1 As shown in the figure, this embodiment describes a method for predicting crack propagation in a turbine runner under transient hydraulic excitation, including the following steps: Step 1: Draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, including the turbine blades; draw the horizontal concentrated field based on the cracks on the turbine blades, including the normal horizontal concentrated field and the tangential horizontal concentrated field; Step 2: Determine the crack front based on the horizontal field, refine the mesh near the crack front to form the target mesh; screen the nodes of the target mesh, identify the reinforcing nodes, determine the type of reinforcing element where the reinforcing nodes are located, and solve the time-series displacement field of the turbine runner under transient hydraulic excitation; Step 3: Based on the time-series displacement field of the turbine runner, calculate the stress intensity factor and energy release rate of the crack front node, and use this as a basis to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation; Step 4: Based on the crack propagation direction and crack propagation rate obtained in Step 3, update and correct the horizontal field to obtain the propagated crack; based on the propagated crack, determine whether a through crack has formed; if a through crack has formed, proceed to Step 5; otherwise, return to Step 2 for iteration. Step 5: Predict the fracture life of the impeller; In this embodiment, the rotating component may further include an upper crown and a lower ring; the normal horizontal field is used to characterize the crack surface, and the tangential horizontal field is used to characterize the crack leading edge; the crack leading edge is determined by calculating the intersection of the zero isosurfaces of the normal horizontal field and the tangential horizontal field.

[0024] The method for determining whether a through crack has formed can be: calculate the distance between the center point of the crack leading edge and the edge of the blade mesh where the crack is located, that is, the distance from the crack tip to the upper crown or lower ring. When this distance is less than the crack penetration determination threshold, it is determined that a through crack has formed. The crack penetration threshold can be set as the maximum crack propagation length in each crack propagation iteration.

[0025] Furthermore, in step two, the mesh near the crack tip is a tetrahedral mesh; the step of refining the mesh near the crack tip to form the target mesh is as follows: Step A1: Determine the target mesh size based on the crack length, and then determine the number of mesh refinements by combining the original mesh size near the crack leading edge with the target mesh size; Step A2: Use the distance control index field to filter out grid cells that meet the encryption requirements; Step A3: Process the mesh cells that meet the encryption requirements, set the midpoints of each of their edges as newly added mesh points, and divide the original tetrahedral mesh cells into four new mesh cells. Step A4: Determine whether the encryption count has been reached. If yes, complete the encryption; otherwise, return to step A2.

[0026] The area near the crack tip can be defined as the region extending from the crack tip along the crack surface in both the normal and tangential directions to 1.5 to 2 times the crack length.

[0027] Furthermore, in step two, the method for filtering the nodes of the target mesh and identifying the enhanced nodes is as follows: Find the mesh elements in the target mesh that cross the zero isosurface at the crack front, and use the nodes of these mesh elements as candidate nodes; Traverse all candidate nodes and filter out crack-enhancing nodes and crack tip-enhancing nodes; the crack-enhancing nodes are used to characterize the displacement discontinuity of the crack surface; the crack tip-enhancing nodes are used to characterize the stress singularity of the crack tip; the specific filtering method can be: Assign initial values ​​to the numerical range of the tangential level set of candidate nodes; The mesh cells containing candidate nodes are evaluated to determine whether the mesh edges of these mesh cells cross the crack leading edge; For mesh edges that pass through the crack front, the range of tangential level set values ​​of the candidate nodes to which these mesh edges belong is updated by directly taking values ​​or by interpolation. Traverse all candidate nodes, classify nodes whose tangential level set values ​​are still in the initial range as ordinary nodes; classify nodes whose tangential level set values ​​include zero as crack enhancement nodes; classify nodes whose tangential level set values ​​conform to the range near the crack tip as crack tip enhancement nodes. In step two, the method for determining the type of enhancement unit where the enhancement node resides is as follows: A mesh element containing at least one crack reinforcement node is defined as a crack reinforcement element; a mesh element containing at least one crack tip reinforcement node is defined as a crack tip reinforcement element; and a mesh element containing both is defined as a hybrid reinforcement element.

[0028] Furthermore, in step two, the method for solving the time-series displacement field of the turbine runner under transient hydraulic excitation is as follows: Based on the horizontal field, the reinforcement nodes near the crack tip are selected and the reinforcement elements in which they are located are determined; Based on the aforementioned enhancement unit, the unit stiffness matrix and mass matrix are calculated, and then assembled into the overall stiffness matrix and mass matrix of the impeller. Based on fluid dynamics simulation, data on pulsating fluid pressure loads under transient hydraulic excitation are obtained. Based on the overall stiffness matrix and mass matrix of the runner, as well as the pulsating fluid pressure load data, the acoustic-structure interaction dynamic equilibrium equation of the runner is solved to obtain the time-series displacement field of the runner.

[0029] Specifically, this can be done by: based on the horizontal field, selecting crack reinforcement nodes and crack tip reinforcement nodes near the crack front, and determining the crack reinforcement unit, crack tip reinforcement unit, and mixed reinforcement unit in which they are located; For an element containing reinforcing nodes, the element is divided into second-order tetrahedral sub-elements using the zero isosurface of the normal horizontal set as the dividing surface. The Gaussian points and weights of the sub-elements are calculated, and the element stiffness matrix and mass matrix are obtained by integration. Then, the elements are assembled into the overall stiffness matrix and mass matrix of the impeller. Based on fluid dynamics simulation, the pulsating fluid pressure load data under transient hydraulic excitation can be obtained. Specifically, based on fluid dynamics simulation, the pulsating fluid pressure data of the fluid-structure interaction surface inside the volute under transient hydraulic excitation can be obtained under the steady-state operation of the turbine. Within the total duration of transient hydraulic excitation, a fixed duration containing multiple load cycles is extracted and taken as the time history of pulsating fluid pressure load. Extract pulsating fluid pressure load data within the time history of the pulsating fluid pressure load; Fluid dynamics simulation may include the following: Using CFD simulation software, a geometric model of the turbine volute fluid domain and rotating components matching the model from step one was first constructed. The turbulence model uses the inlet velocity of the spiral casing and the outlet pressure of the tailrace pipe, which are monitored under the current operating power of the turbine, as boundary conditions for the computational fluid dynamics model to carry out simulation. The assembled overall stiffness matrix and mass matrix, combined with fluid domain parameters (including fluid domain stiffness matrix and mass matrix), structural and fluid damping matrix, fluid-structure interaction matrix, and pulsating fluid pressure load, are substituted into the acoustic-structure interaction dynamic equilibrium equation of the runner. The displacement data of each node of the runner at each time moment are obtained by solving the equation using the finite element numerical method, thus obtaining the time-series displacement field of the runner.

[0030] Furthermore, in step three, the step of calculating the stress intensity factor and energy release rate of the crack front node based on the time-series displacement field of the turbine is as follows: Step B1: Obtain the location of the crack tip, take several nodes that are evenly distributed near the crack tip as the original calculation nodes, and establish a local coordinate system of the crack tip that matches the geometric direction of the crack tip. Step B2: Obtain the displacement field at each moment based on the time-series displacement field of the rotor, and extract the displacement results on both sides of the crack surface at each moment; Step B3: Transform the displacement results obtained in step B2 into the local coordinate system of the crack front, and then project them onto the original calculation nodes to form projected calculation nodes; Step B4: Based on the projection calculation nodes obtained in Step B3, draw a scatter plot and obtain the stress intensity factor of the crack front node at each time point through the scatter plot. Step B5: Based on the stress intensity factor obtained in step B4, the energy release rate of the crack front node at each time step is calculated using the Evan formula.

[0031] Crack leading edge nodes refer to finite element mesh nodes located on the boundary line of the crack leading edge; The number of original computing nodes can be four to eight; The area near the crack tip can be defined based on the mechanical properties of the crack tip (the core micro-region where stress is most concentrated on the crack front) and the following criteria can be used: it can be a region centered on the crack tip with a distance not exceeding 0.1 times the crack length, or a region centered on the crack tip with a distance not exceeding 4 times the target mesh size.

[0032] The method for establishing a local coordinate system at the crack tip can be as follows: taking the area near the crack tip where the original calculation node is located as the reference, and the crack tip as the origin, an orthogonal local coordinate system is established along the crack surface normal, crack surface tangent and crack tip axis to ensure that subsequent displacement transformation and projection can accurately match the crack force direction. The displacement results can be transformed into the local coordinate system at the crack front by coordinate system transformation.

[0033] The scatter plot can be drawn as follows: the constant coefficient is plotted on the horizontal axis, and the displacement difference between the two sides of the crack surface corresponding to each projection calculation node is plotted on the vertical axis; wherein, the constant coefficient is determined by the crack type (Type I opening type, Type II shear type) and derived from the analytical expression of crack tip stress-displacement; for Type I cracks, the displacement difference perpendicular to the crack surface is selected, and for Type II cracks, the displacement difference along the tangential direction of the crack surface is selected.

[0034] The stress intensity factor can be determined using a scatter plot by the following method: The maximum intercept of the straight line equation formed by two adjacent points is taken as the stress intensity factor. The maximum slope of the straight line formed by each point and the origin is taken as the stress intensity factor. The coefficient of the first-order term is taken as the stress intensity factor by least squares regression calculation. To reduce calculation errors and to take into account a certain degree of conservatism, the stress intensity factor can be taken as the maximum value of the three methods mentioned above as the calculation result.

[0035] In step three, the method for calculating the stress intensity factor and energy release rate of the crack front node based on the time-series displacement field of the rotor can also employ another high-precision calculation method. This method is based on the Lagrange differential of the structural potential energy in the neighborhood of the crack front, utilizes the Sita field describing crack propagation to solve for the energy release rate of the crack front node, and further calculates the stress intensity factor. This is a more accurate calculation method applicable to arbitrary topologies. Specific steps may include: The first step is to take the Lagrange differential principle of structural potential energy near the crack tip as the core and assume a Westerner field describing crack propagation. This Westerner field satisfies the condition that it is adapted to the tangential direction vector of the crack tip. The second step is to define a spatial transformation using the assumed Westerland field, which transforms any spatial point in the computational domain into a transformed spatial point. The third step is to construct a weak form equation involving the energy release rate of the crack front, computational domain, and crack surface, which holds for any Theta field. The fourth step is to select multiple West Tower fields that meet the conditions, and at the same time use the crack front shape function to discretize the energy release rate in the equation, and solve the equation to obtain the energy release rate of the crack front node at each time. Fifth step: For any given Westerland field that meets the conditions, define a bilinear function, which acts to separate the three cracking modes of opening, shearing and tearing. The sixth step is to construct the weak form equation corresponding to the stress intensity factor at each moment of the flow-induced vibration. This equation incorporates the material's elastic modulus, Poisson's ratio, crack front tangential direction vector, crack front geometric information, as well as the temporal displacement field and the displacement fields corresponding to the three cracking types under the assumption of an infinitely large structural plane crack in the neighborhood of the crack front. The seventh step involves discretizing the stress intensity factor using the same shape function as the energy release rate, and using the same Theta field to solve the weak form equation corresponding to the stress intensity factor, thereby obtaining the stress intensity factor at the crack front node at each time step.

[0036] Furthermore, in step three, the step of calculating the crack propagation direction under transient hydraulic excitation is as follows: Step C1: Based on the stress intensity factor of the crack front node at each moment under transient hydraulic excitation, the maximum circumferential stress criterion is used to calculate the deflection angle of each crack front node at each moment; the deflection angle is located in the plane formed by the crack front tangent and the crack surface normal. Step C2: Collect the deflection angle data of all crack leading edge nodes at each time step; Step C3: Calculate the average deflection angle data collected in step C2 over time. The average deflection angle obtained is the crack propagation direction. In step three, the method for calculating the crack propagation rate under transient hydraulic excitation is as follows: Calculate the corresponding equivalent stress intensity factor using the energy release rate of the crack front node at each time point; The equivalent stress intensity factor at each time point is statistically analyzed across the entire range, and the amplitude of the equivalent stress intensity factor corresponding to all load cycles is extracted. Obtain the crack propagation threshold value of the rotating wheel; The amplitude of the equivalent stress intensity factor corresponding to each load cycle is judged. When it is greater than the crack propagation threshold, the crack propagation rate is calculated by the Paris formula; otherwise, the assumed value is forcibly set as the crack propagation rate.

[0037] The deflection angle is located in the plane formed by the tangential direction of the crack leading edge and the normal direction of the crack surface. This plane is the core area where both Type I (opening) and Type II (shearing) cracking modes work together, and it is fully compatible with the crack propagation law of the thin-shell structure of the turbine blade. Considering the characteristics of the thin-shell structure of the turbine blade, the crack torsion problem that may be caused by the Type III (tearing) cracking mode is excluded. The focus is only on the plane formed by the tangential direction of the crack leading edge and the normal direction of the crack surface. All subsequent deflection direction calculations are carried out within this plane.

[0038] The stress intensity factor at the crack front at various times is used as the core basic data. These data cover different time points during the operation of the turbine, ensuring that they can reflect the dynamic changes of the stress intensity factor under transient hydraulic excitation.

[0039] In the process of calculating the deflection angle, the direction attribute of the deflection angle can be clearly defined by the sign function. The direction attribute can include positive and negative directions to accurately reflect the possible deflection trend of the crack.

[0040] Using the obtained average deflection angle as the crack propagation direction ensures that the results can comprehensively reflect the influence of stress state on crack propagation direction at different times, and are more in line with actual operating conditions.

[0041] The equivalent stress intensity factor amplitude corresponding to a load cycle refers to the difference between the maximum and minimum values ​​of the equivalent stress intensity factor within each independent load cycle.

[0042] Forcing a crack propagation rate as the assumed value can induce forced crack propagation. This technique is used for cracks that do not meet the physical propagation conditions under current hydraulic excitation conditions, where the equivalent stress intensity factor amplitude has not reached the fatigue propagation threshold. By setting an assumed crack propagation rate (not a physically meaningful rate), it forces the initiation of iterative crack propagation calculations. Its core purpose is to continuously simulate the evolution of potential subsequent crack morphologies. Without this forced propagation mechanism, cracks that stagnate due to not meeting the propagation conditions cannot have their subsequent morphology obtained through simulation iterations, thus making it impossible to complete the later-stage stability assessment of crack propagation (such as determining whether unstable propagation has occurred) and the prediction of critical crack length.

[0043] The method for determining the assumed value can be: under the premise that the crack propagation threshold is 0, the propagation rate calculated by the Paris formula is the assumed value.

[0044] Further, step four: updating and correcting the horizontal field based on the crack propagation direction and crack propagation rate obtained in step three to obtain the propagated crack step is as follows: Step D1: Based on the tangential level set, divide all meshes into cracked and uncracked regions; For the already cracked area, the propagation rate of the normal horizontal set is fixed to zero, so that the existing crack surface remains unchanged. For the uncracked zone, the crack propagation rate is decomposed into normal propagation rate and tangential propagation rate based on the crack propagation direction obtained in step three. Step D2: Based on the normal propagation rate and the tangential propagation rate, update the horizontal field after crack propagation in the uncracked zone to characterize the updated crack surface and crack front topology, so as to obtain the expanded crack. Step D3: Correct the horizontal field obtained after crack propagation in step D2, ensuring that the absolute value of the new horizontal gradient field is always equal to 1, and that the normal horizontal field and the tangential horizontal field on the crack surface are orthogonal.

[0045] Step D2's core is to transform the physical process of crack propagation in the uncracked zone into numerical changes in the horizontal field, thereby accurately reflecting the crack morphology after propagation. The specific steps of Step D2 can be as follows: First, based on the normal and tangential propagation rates: the tangential rate determines the speed at which the crack propagates forward, and the normal rate determines the speed at which the crack deflects laterally; then, combined with the time span of the current iteration cycle, update the horizontal field values ​​of all mesh nodes in the uncracked zone. The numerical changes reflect the relative position of the nodes and the new crack; for example, if the values ​​of some nodes change from positive to negative, it means that the region where that node is located has changed from uncracked to cracked; finally, update the horizontal field to clearly define the new crack surface range (tangential horizontal field zero isosurface) and the crack leading edge position (normal horizontal field zero isosurface), providing accurate geometric basis for subsequent crack state judgment and the next iteration.

[0046] The specific steps of step D3 can be as follows: First, the gradient field of the normal level set is normalized; The zero isosurface of the tangential level set is then corrected; Finally, the gradient of the tangential level set is normalized so that the absolute value of the gradient field of the tangential level set is always 1.

[0047] Furthermore, step five, the method for predicting the fracture life of the impeller, includes: Extract the crack length, maximum equivalent stress intensity factor, and maximum equivalent stress intensity factor amplitude in each iteration step; Plot the first relationship curve between the maximum equivalent stress intensity factor and the crack length. Based on this first relationship curve and the fracture toughness of the runner material, determine whether the runner blade will undergo brittle fracture. If there is a possibility of brittle fracture, determine the critical crack length for brittle fracture and provide monitoring prompts. A second relationship curve is plotted between the amplitude of the maximum equivalent stress intensity factor and the crack length. Based on this second relationship curve and the fatigue crack propagation threshold value of the impeller material, it is determined whether fatigue fracture will occur due to fatigue crack propagation. When fatigue fracture is possible, the critical crack length at which propagation begins is determined, and the time from the critical crack length at which propagation begins to the formation of a through crack is calculated. The calculation result is the fracture life of the impeller.

[0048] The method for judging the possibility of fatigue fracture is: the amplitude of the maximum equivalent stress strength factor is greater than the fatigue crack propagation threshold value of the wheel material.

[0049] The calculation of the time from the critical crack length at which propagation begins to the formation of a through crack can be based on the number of load cycles required to form a through crack, the number of load cycles corresponding to the critical length, and the time history of the pulsating fluid pressure load; the specific steps can be: First, relevant data are extracted from each crack propagation iteration step, including crack length, maximum equivalent stress intensity factor, maximum equivalent stress intensity factor amplitude, number of load cycles, and cumulative number of load cycles in each iteration step. Specifically, within a single crack propagation iteration step, the current crack length is fixed, and the equivalent stress intensity factor at each moment of transient hydraulic excitation under that crack length is calculated. The maximum value among these values ​​is selected as the maximum equivalent stress intensity factor for that iteration step. Within a single crack propagation iteration step, the amplitude of the equivalent stress intensity factor at the crack front corresponding to each load cycle is obtained after full-amplitude statistics. The maximum value among these values ​​is selected as the amplitude of the maximum equivalent stress intensity factor for that iteration step. Plot a third relationship curve between the cumulative number of load cycles and the crack length. Based on this third relationship curve, determine the number of load cycles corresponding to the formation of a through crack and the number of load cycles corresponding to the critical crack length. Then calculate the difference between the two load cycle counts. This difference is the number of load cycles required for the crack to expand from the critical length to a through crack. By obtaining the duration of each load cycle through the time history of the pulsating fluid pressure load, and multiplying this duration by the required number of load cycles, we obtain the time from the critical crack length at which the crack begins to expand to the formation of a through crack, which is the fracture life of the impeller.

[0050] This embodiment, through multiple iterations, can track the crack's propagation process from its current minor state to full penetration, and can accurately predict the rotor's fracture life based on its propagation pattern.

[0051] Example 3 This embodiment describes a method for predicting crack propagation in a turbine runner under transient hydraulic excitation, including the following steps: Step 1: Draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, including the turbine blades; draw the horizontal concentrated field based on the cracks on the turbine blades, including the normal horizontal concentrated field and the tangential horizontal concentrated field; Specifically, based on the turbine runner design drawings, a geometric model including the runner blades, upper crown, and lower ring can be drawn, and a first-order linear mesh can be divided. The mesh can be appropriately refined at component connections, areas with large geometric curvature, and blades where initial cracks are located, thereby improving the accuracy of finite element displacement calculations at the corresponding locations.

[0052] Based on the known crack configuration information, normal and tangential horizontal fields are created in the finite element mesh to characterize the crack lead and crack surface, respectively. The horizontal field calculation method is as follows: for each mesh node, the values ​​of the normal and tangential horizontal fields are taken as the directed projection distances of the node to the crack lead in the local crack coordinates, i.e., for any mesh node P:

[0053] Point M is the point on the crack front edge that is closest to point P. This represents the vector originating from M and pointing to P. and Let M represent the local coordinate normal and tangential unit vectors of the crack at point M, where the tangential vector is perpendicular to the tangent at the crack tip and points in the opposite direction to the crack surface. Represents the normal level set, This represents the tangential level set.

[0054] Step 2: Determine the crack front based on the horizontal field, refine the mesh near the crack front to form the target mesh; screen the nodes of the target mesh, identify the reinforcing nodes, determine the type of reinforcing element where the reinforcing nodes are located, and solve the time-series displacement field of the runner under transient hydraulic excitation.

[0055] The crack front can be obtained by calculating the intersection of the zero isosurfaces of the normal and tangential horizontal concentrated fields.

[0056] First, based on the known crack configuration information and finite element mesh, the maximum mesh size and mesh refinement level near the crack tip that meet the requirements of fracture mechanics numerical simulation are calculated. Specifically, the calculation method is as follows: Let the original mesh size near the crack tip be... The crack length is Then the target mesh size for mesh refinement should be smaller than The required number of mesh encryption attempts is:

[0057] The final mesh size near the crack leading edge obtained after mesh refinement is: .

[0058] Then, the mesh is refined to ensure that the mesh size meets the requirements of fracture mechanics numerical simulation. Specifically, let the current mesh refinement step number be... Using distance control index field During the control grid densification process, the control index field value for all cell nodes is greater than... The original tetrahedral mesh element is divided into four new mesh elements by taking the midpoints of each edge as the new mesh points, and the nodes of the refined mesh are renumbered. Then return to continue encryption.

[0059] Next, based on the horizontal field, the extended finite element method identifies the reinforcement nodes and their reinforcement types, and then identifies the type of reinforcement element of the element to which they belong. The dynamic process of a cracked turbine blade has significant discontinuous and nonlinear characteristics. Therefore, the extended finite element method defines two types of reinforcement nodes near the crack: crack reinforcement nodes and crack tip reinforcement nodes, based on step functions and crack tip reinforcement functions with singularities. These are used to characterize the displacement discontinuity of the crack surface and the stress singularity of the crack tip in a linear dynamic system, respectively.

[0060] This embodiment proposes two enhanced node identification algorithms with different mesh convergence and time complexity. The main difference between the two algorithms lies in the identification range of enhanced nodes at the crack leading edge.

[0061] The first algorithm is as follows: for any grid node Define condition one:

[0062] in Represents all containing nodes The set of nodes for the grid cells. Initialization. and These are used to record the minimum and maximum values ​​of the tangential level set, respectively. Traverse the edges of the mesh cells containing the node P that satisfies condition one, denoted by endpoints A and B. If there exists an endpoint with a normal level set of zero, update the tangential level set according to the value at that endpoint. and ;like Then, based on the element order, the value of the tangential level set at the intersection point C of AB and the zero isosurface of the normal level set is calculated using linear or quadratic polynomial interpolation, and then updated. and The calculation method for linear interpolation is as follows:

[0063] The calculation method for quadratic interpolation is as follows:

[0064] Where M is an intermediate node. For the equation The solution.

[0065] If the minimum and maximum values ​​of the tangential level set obtained after traversal are still the initial values, then node P is a normal finite element node; if it satisfies Then node P is a crack reinforcement node based on the step function; If satisfied Then node P is a crack tip enhancement node.

[0066] Related research indicates that appropriately increasing the recognition range of crack tip reinforcement nodes helps improve the mesh convergence rate, i.e., reducing the mesh density required to achieve the same computational accuracy. Therefore, the second algorithm defines the geometric recognition range of crack tip reinforcement nodes. The decision condition for the crack tip enhancement node in the first algorithm is replaced with... ,in It is approximately 1 / 10 of the crack length, or 4 times the crack leading edge mesh size after mesh refinement.

[0067] After identifying the enhancement nodes, the enhancement units are defined, specifically including: A crack-enhancing element is defined as an element with at least one crack-enhancing node; a crack tip-enhancing element is defined as an element with at least one crack tip-enhancing node; and a hybrid-enhancing element is the intersection of the two.

[0068] Finally, considering the extended finite element reinforcement nodes near the crack, the structural stiffness and mass matrix are calculated and assembled based on the finite element method. The time-series displacement field is obtained by solving the dynamic equilibrium equation of the turbine runner according to the time history of the pulsating fluid pressure load under the steady-state operation condition of the turbine.

[0069] The time history of pulsating fluid pressure load can be obtained through fluid dynamics simulation, which may include the following: using CFD simulation software, first constructing a geometric model of the turbine volute fluid domain and rotating components that matches step one, and selecting... The turbulence model uses the inlet velocity of the spiral casing and the outlet pressure of the tailrace pipe, which are monitored under the current operating power of the turbine, as boundary conditions for the computational fluid dynamics model to carry out simulation.

[0070] Specifically, in the extended finite element method, the approximate expression for the solution of the displacement field is:

[0071] in Representing three-dimensional coordinates The set of nodes in the grid cell. Represents the set of crack-enhancing nodes. This represents the set of nodes that enhance crack tip reinforcement. Represents the unit shape function, This represents the extended finite element degrees of freedom (i.e., the unknowns to be solved), which are associated with ordinary finite element nodes, crack-enhanced nodes, and crack tip-enhanced nodes, respectively. This represents the step function, used to distinguish the relative positions of nodes and crack surfaces, and characterizes the displacement discontinuity near the crack surface. It is a basis function whose first derivative tends to infinity at the crack tip, characterizing the stress concentration at the crack tip.

[0072] The displacement field can be approximated in matrix form:

[0073] in It is a unit shape function matrix; This is the finite element degree-of-freedom vector. Substituting it into the weak form of the dynamic equilibrium equations, when calculating the structural stiffness and mass matrix, since crack-reinforced elements, crack-tip-reinforced elements, and hybrid-reinforced elements all contain at least one reinforcing node, their shape functions... or Due to its singularity, it cannot be directly integrated numerically using the Gauss-Legend method. Therefore, the algorithm uses the zero isosurface of the normal level set as the dividing surface to segment the reinforcement element into second-order tetrahedral sub-elements, calculates the Gaussian points and their weights of the elements, and uses the Gauss-Legend method to integrate and calculate the element stiffness and mass matrix, i.e.:

[0074]

[0075] in Represents the integration domain of the sub-unit. Represents the gradient matrix of the shape function The transpose of the matrix, Represents the shape function matrix The transpose matrix x, Represents the material elasticity matrix. Indicates the density of the material. Represents the Gaussian point of the sub-unit. Indicates the Gaussian point index. Indicates the Gaussian point weights. This represents the deformation gradient matrix of the sub-unit. This represents the element stiffness matrix and mass matrix.

[0076] The time-series displacement field of the turbine runner is obtained by solving the acoustic-structure interaction dynamic equations of the turbine runner using the finite element method:

[0077] in: The overall stiffness and mass matrix of the impeller are assembled from the element stiffness and mass matrix. The element stiffness matrix and mass matrix can be obtained from the aforementioned formula. During assembly, according to the connection relationship of the element nodes and the global degree of freedom number, the matrix items of each element are accumulated to the corresponding position of the overall matrix. For crack-reinforced elements, crack tip-reinforced elements and hybrid-reinforced elements, their corresponding extended degrees of freedom are uniformly incorporated into the global degree of freedom system and the assembly is completed.

[0078] These are the fluid domain stiffness and mass matrix, which are related to the fluid density and compressibility modulus. The fluid domain mass matrix and stiffness matrix are obtained by discretizing and assembling the fluid acoustic pressure control equations on the fluid domain mesh in the structural finite element solution model.

[0079] These are the structural and fluid damping matrices, which can be constructed using conventional damping models in this field, such as proportional damping models.

[0080] This is the coupling matrix associated with the fluid-structure interaction surface element; the coupling matrix is ​​used to realize the discrete mapping of fluid-structure interface variables (pressure, displacement, etc.), and its structure is determined by the coupling surface mesh and shape function (structure mesh), which belongs to conventional finite element coupled discretization; For the displacement field of the rotating wheel, The two are the acoustic pressure field in the fluid domain and the variables to be solved in the acoustic-structure interaction dynamic equation.

[0081] The blade pressure load is obtained from the pulsating pressure distribution on the fluid-structure interaction surface obtained by CFD simulation. This pressure load is then equated to a structural pressure field through interpolation and projection from the fluid computational grid to the structural grid. The interpolation and projection are conventional load transfer methods in fluid-structure interaction calculations. Specifically, for each structural grid node, its corresponding fluid element is located in the fluid grid, and the node pressure value is obtained by interpolation using the shape function of that fluid element. The pressure values ​​at each node of the structural element are then equated to the structural element pressure load.

[0082] Represents the second and first partial derivatives with respect to time; The above process can be implemented using finite element method (FEM) software.

[0083] Step 3: Based on the time-series displacement field of the turbine runner, calculate the stress intensity factor and energy release rate of the crack front node, and use this as a basis to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation; First, based on the time-series displacement field results of the cracked impeller, the stress intensity factor and energy release rate at the crack front are calculated. This embodiment proposes two calculation methods with different levels of accuracy, efficiency, and applicability: The first calculation method is the displacement extrapolation method: In linear elastic, homogeneous, and isotropic two-dimensional cracking problems, the stress intensity factors of the three cracking modes can be derived from the stress-displacement analytical solution at the crack tip. This method is also applicable to the small neighborhood of cracks in thin-shell structures such as turbine blades, and can be written in limit form:

[0084] in Let M be the type I (opening) or II (shear) stress intensity factor, M be the crack front node, and r be the distance from M along the crack surface. This represents the displacement difference between the two sides of the crack in a direction perpendicular to the crack surface or the crack tip. The difference in tangential displacement between the two sides of the crack. The elastic modulus of the material. It is Poisson's ratio.

[0085] Based on the displacement field obtained in step two, the maximum computation distance is taken as 4 to 6 times the mesh size near the crack tip. This distance is then divided into 4 to 8 nodes as computation nodes. The displacement results on both sides of the crack surface are extracted, transformed to the local coordinate system of the crack front, and projected onto the computation nodes. The stress intensity factor is calculated using three different methods: Method 1: Drawing about The scatter plot, where These are constant coefficients derived from analytical expressions, for example, for type I cracks. The maximum intercept of the straight line equation formed by two adjacent points is taken as the stress intensity factor. Here, U represents the displacement difference across the crack surface.

[0086] Method 2: Drawing about The scatter plot is used to determine the stress intensity factor by taking the maximum slope of the straight line formed by each point and the origin.

[0087] Method 3: Drawing about The scatter plot was used, and the coefficient of the first term was taken as the stress intensity factor by least squares regression calculation.

[0088] To reduce calculation errors and to consider a certain degree of conservatism, the maximum stress intensity factor obtained from the three methods is taken as the final calculation result. Finally, the energy release rate at each crack front node is calculated according to the Evan formula.

[0089] in, It is a type I (opening) or II (shear) stress intensity factor; It is a Type III (tear) stress intensity factor.

[0090] The second method is the G-theta method: based on the Lagrange differential of the structural potential energy in the neighborhood of the crack front, it uses the Theta field describing crack propagation to solve for the energy release rate at each computational node of the crack front. This is a high-precision calculation method with higher accuracy and applicability to arbitrary topologies, but its computational efficiency is lower than that of the displacement extrapolation method. Specifically, using the assumed Theta field, a spatial transformation is defined:

[0091] in Representation and parameters Related spatial transformations, For any point in space, For the computational domain, for Transformed spatial points, To indicate the western field, it must satisfy the condition at the crack leading edge. , Let be the tangential direction vector at the crack tip. Solve the weak-form equation satisfied by the energy release rate:

[0092] in , and This represents the crack tip, computational domain, and crack surface. Indicates the energy release rate. This represents the tangential direction vector at the crack tip. Indicates strain energy. Representing volume forces and surface forces, This represents the displacement field after spatial transformation. The equation holds for any Sita field, and the right-hand side is related to the solution of the displacement field. This is achieved by selecting P Sita fields that satisfy the conditions, and simultaneously selecting N crack front shape functions to transform the left-hand side... Discrete, in satisfying Under these conditions, the energy release rate at each calculated node at the crack front can be obtained.

[0093] To calculate the stress intensity factor at the crack tip, a bilinear function is defined for any given Westerland field that satisfies the following conditions:

[0094] in For time-series displacement fields, This represents the auxiliary displacement field used to extract the stress intensity factor of the corresponding cracking mode. Specifically, the auxiliary displacement field is the analytical displacement field of the corresponding cracking mode within the crack front neighborhood under the assumption of an infinitely large structural plane crack. The cracking modes are represented by values ​​I, II, and III. A bilinear function is chosen to separate the three cracking modes: opening, shearing, and tearing. At each moment of flow-induced vibration, the weak-form equation satisfied by the stress intensity factor is solved:

[0095] in The elastic modulus of the material. Poisson's ratio, This represents the tangential direction vector at the crack tip. Indicates the leading edge of the crack. It is a type I (opening) or type II (shear) stress intensity factor. For time-series displacement fields, This represents the displacement field corresponding to the three cracking types under the assumption of a planar crack within an infinitely large structure in the neighborhood of the crack front. s represents the crack surface. For Xitachang, It is a bilinear function constructed from the generalized function of energy release rate.

[0096] The stress intensity factor was discretized using the same shape function as the energy release rate, and the same Westerland field was selected to obtain the stress intensity factor of each calculation node at the crack front.

[0097] Next, the crack deflection direction is calculated based on the stress intensity factor at each moment at the crack tip. Considering the thin-shell structure of the turbine blade, crack torsion caused by Type III (tear-off) cracking mode is not considered; only crack deflection in the plane formed by the tangential direction at the crack tip and the normal direction of the crack surface is considered, and the deflection angle is calculated according to the maximum circumferential stress criterion. The specific calculation method is as follows:

[0098] in Indicates the deflection angle. Indicates the sign function, The time value is used to represent the crack propagation direction. The time average of the deflection angles at each calculation node is taken as the final crack propagation direction.

[0099] Then, the crack propagation rate under transient hydraulic excitation is calculated through load cycle statistics. Specifically, the equivalent stress intensity factor is first calculated using the energy release rate:

[0100] in It is the equivalent stress intensity factor. Energy release rate, The elastic modulus of the material. Poisson's ratio, Indicates time.

[0101] Subsequently, the rainflow counting method was used to perform full-range statistical analysis of the equivalent stress intensity factors at each calculation node of the crack front, and the crack propagation rate was calculated according to the Paris formula:

[0102] in Indicates the crack length. This represents the number of hydraulic excitation load cycles used to calculate the time-series displacement field. The constant parameter of the Paris formula, This represents the first result obtained by the rainflow counting method. The amplitude of the equivalent stress intensity factor This is the crack propagation threshold value. To assess the crack propagation path and crack stability at each stage of propagation under the current operating condition, when... hour, This strategy forces crack propagation even when the stress intensity factor does not meet the fatigue crack propagation condition. The impact of this calculation strategy is that the propagation rate calculated using the Paris formula during forced crack propagation is an assumed rate and has no physical meaning. Furthermore, considering that the hydraulic excitation load on turbine blades generally has approximately proportional loading characteristics, meaning the change in the principal stress direction in the crack front neighborhood with time is negligible, it is intuitively reflected in the crack deflection angle. The approximate time independence of the crack makes it possible to obtain the correct crack propagation path even for cracks that would not normally propagate under the current conditions, by applying a hypothetical forced crack propagation. This helps to evaluate whether unstable cracks will occur in the later stages of crack propagation or to predict the critical crack length at which cracks actually begin to propagate through simulation.

[0103] Step 4: Based on the crack propagation direction and crack propagation rate obtained in Step 3, update and correct the horizontal field to obtain the propagated crack; based on the propagated crack, determine whether a through crack has formed; if a through crack has formed, proceed to Step 5; otherwise, return to Step 2 for iteration. First, in the mesh of the turbine runner blades with cracks, the horizontal set propagation rate of all mesh nodes is calculated based on the crack propagation rate and deflection direction angle obtained in step three. Specifically, to ensure that existing cracks remain unchanged in the current crack propagation iteration step, the tangential horizontal set symbol is used to distinguish the cracked regions ( ) and uncracked areas ( ), calculate the normal and tangential crack propagation rates respectively. For the already cracked region, fix the normal horizontal propagation rate to zero, thus keeping the existing crack surface unchanged, and calculate the tangential horizontal propagation rate as follows:

[0104] in This represents the tangential crack propagation rate at the crack front projection of mesh node M. This represents the gradient vector of the tangential horizontal field at node M. For the uncracked region, the horizontal set propagation rate is calculated as follows:

[0105] in , For the normal and tangential horizontal set propagation rates, This represents the normal and tangential crack propagation rates at point P, where mesh node M is projected at the crack front. and Let represent the normal and tangential local coordinate unit vectors at point P at the crack tip. Indicates the maximum load cycle count. The maximum crack propagation length is fixed for each iteration.

[0106] Secondly, update the crack propagation. The horizontal set of lengthened parameters is used to characterize the updated crack surface and crack front topology. Specifically, the updated normal or tangential horizontal set is calculated as follows:

[0107] in Indicates a normal or tangential horizontal field; The propagation rate field representing the normal or tangential horizontal field; This represents the gradient vector of the level set.

[0108] Then, the horizontal field after crack propagation is modified to ensure that the absolute value of the new horizontal gradient field is always equal to 1. Furthermore, the normal and tangential horizontal field on the crack surface are orthogonal. This serves two purposes: first, to ensure that the normal and tangential direction vectors and displacements of the crack front can be correctly obtained from the horizontal field during the calculation of stress intensity factor and energy release rate in the next crack propagation calculation; second, to more accurately deduce the horizontal propagation rate of the grid nodes of the turbine blade with crack from the crack propagation rate in the next crack propagation calculation.

[0109] Considering that it is impossible to simultaneously satisfy the normalization and orthogonality of the horizontal set gradient field for cracks with large surface curvature, the normalization condition should be satisfied first, while orthogonality should be satisfied in the neighborhood of the crack tip. Therefore, the following technical solution is adopted: The gradient field of the normal level set is normalized using a fast-progression method. Starting from the crack surface nodes, partial differential equations are solved using a spatial finite-difference scheme towards both sides of the crack surface. :

[0110] in This represents the spacing between adjacent grid nodes in the x, y, and z directions. This represents the corrected normal level set values ​​for the mesh nodes with indices i, j, and k in the x, y, and z directions, respectively. The constants related to i, j, k are calculated as follows:

[0111] Based on the characteristics of level sets Thus ensuring the above regarding The quadratic polynomial equation has a unique solution.

[0112] The zero isosurface of the tangential level set is corrected. Specifically, based on the updated tangential level set field, each grid node near its zero isosurface is selected, and the tangential level set value is recalculated according to the corrected normal level set and the updated crack front, so that the zero isosurfaces of the tangential and normal level sets are orthogonal at the crack front.

[0113] Starting from the node near the zero isovalue of the original tangential level set after correction, the tangential level set values ​​of the remaining grid nodes are corrected using the same fast-advance method as the normal level set, so that the absolute value of the gradient field of the tangential level set is always 1.

[0114] Finally, based on the expanded crack, a determination is made as to whether a through crack has formed. The distance between the updated crack leading edge center point and the edge of the blade mesh containing the crack is calculated, i.e., the distance from the crack tip to the upper cap or lower ring. If this distance is less than... If a through crack has been formed, proceed to step five; otherwise, return to step two for the next crack propagation iteration.

[0115] in, This refers to the maximum crack propagation length at each iteration.

[0116] Step 5: Predict the fracture life of the turbine; this includes the following steps: First, the crack length in each crack propagation iteration step is extracted and calculated. Maximum equivalent stress intensity factor Maximum equivalent stress intensity factor variation Number of loops and its cumulative value . It is the equivalent stress intensity factor; Indicates the first The amplitude of the equivalent stress intensity factor.

[0117] Secondly, the stability of crack propagation is determined and the fracture life of the runner blade is calculated.

[0118] Specifically, first draw The curve shows that the maximum equivalent stress strength factor is greater than the fracture toughness of the turbine material. ( Using this as a criterion, it is determined whether the turbine blade will undergo brittle fracture during the initial crack and each stage of crack propagation. If the curve indicates that the turbine blade has the potential for brittle fracture under the current calculated operating conditions, further action is required based on... The curve determines the critical crack length at which brittle fracture occurs. In actual operation, the crack distribution is identified and monitored in a timely manner to ensure that the crack length is much smaller than the critical crack length at which brittle fracture occurs.

[0119] Then draw and The curve is such that the amplitude of the maximum equivalent stress intensity factor is greater than the fatigue crack propagation threshold of the wheel material. ( Using fatigue crack propagation as a criterion, it is determined whether fatigue fracture will occur due to fatigue crack propagation. When the curve indicates the possibility of fatigue fracture, the critical crack length at which propagation begins is determined, and the runner fracture life is calculated, which is the time required for the crack to propagate from the critical crack length to the formation of a through crack. ,in This indicates the number of cycles required for a through-crack to form. express The number of cycles corresponding to the critical crack length in the curve. This indicates the physical duration of the pulsating fluid pressure load. During the actual operation of a water turbine, existing cracks that have reached the critical crack length for fatigue fracture should be repaired within their fracture lifespan. A reasonable predictive maintenance plan should be developed to avoid the formation of through-cracks.

[0120] This embodiment addresses the long-term fracture stability and life prediction problem of turbine runner blades under transient hydraulic excitation. It proposes a method for crack propagation simulation and fracture life prediction without requiring the creation of an extended finite element method (FE) framework with cracked meshes. The algorithm uses an acoustic-structure interaction dynamics method to solve the displacement response of the cracked runner blade under transient hydraulic excitation. Time-varying stress intensity factors and energy release rates are calculated using displacement extrapolation or the G-theta method. Crack deflection angle and propagation rate are calculated using the maximum circumferential stress criterion and rainflow counting method. Then, through level set propagation and updating, and with the normalization and orthogonality correction of the level set gradient field using a fast propagation algorithm, iterative solutions and predictions of the entire crack propagation process are achieved. Finally, the critical crack length and fracture life predictions for brittle and fatigue fractures are given, providing a valid reference for the monitoring, management, and predictive maintenance planning of turbine runner cracks.

[0121] Using the crack propagation prediction method for turbine runners under transient hydraulic excitation described in this embodiment, crack propagation simulation and fracture life analysis were conducted on the runner blades of a mixed-flow turbine at a hydropower station in the Fenshui River Basin of Zhejiang Province. Specifically, the crack propagation simulation and fracture life analysis were performed on the initial crack on the crown side of the blade at rated power. Due to its compact structure, high efficiency, and wide range of applicable heads, the mixed-flow turbine is one of the most widely used turbine types in the world. This embodiment selects a 30MW mixed-flow turbine runner at a hydropower station in the Fenshui River Basin of Zhejiang Province for flow-induced vibration simulation analysis at rated power. The turbine has a diameter of 1.7m, and the runner is a cast-welded structure, formed by pressing S135 (ZG0Cr13Ni5Mo) steel plate. The upper crown and lower ring are made of ZG0Cr16Ni5Mo. The turbine has 24 movable guide vanes and 15 runner blades, a rated speed of 500 rpm, a design head of 186.08m, and a design flow rate of... Due to its relatively thin blade thickness and the large shear moment it bears, an initial crack approximately 2 cm long was found on the outlet edge of blade No. 2, about 10 cm from the upper crown. The fracture toughness of the S135 high-strength steel used in the turbine runner blades is approximately 118. The crack propagation threshold is approximately 3.65. Specifically, the original mesh size near the crack front is The crack length is The target mesh size for mesh encryption is Calculations show that 5 mesh encryption cycles are required.

[0122] With each crack propagation, the mesh needs to be refined accordingly. In this embodiment, the refined mesh result of the last crack propagation calculation is as follows: Figure 2 As shown.

[0123] In this embodiment, the fluid pressure field is obtained through computational fluid dynamics simulation. The turbulence model uses the inlet velocity of the spiral casing monitored at the current operating power of the turbine. The tailrace outlet pressure was set to 0 MPa (relative to atmospheric pressure) as the boundary condition for the computational fluid dynamics model. A time step of 0.001 s was selected, and the steady flow pressure and velocity were solved within a total calculation time of 3.6 s (approximately 30 rotation cycles). The pressure results at the fluid-structure interaction surface between 1.2 s and 3.6 s were used as the hydraulic excitation load for the cracked turbine blades.

[0124] Paris formula constant parameters The maximum crack propagation length is fixed at each iteration. .

[0125] The method described in this embodiment is used to analyze and statistically analyze the mixed-flow turbine of a hydropower station, and the results are... curves and curve; Curves, such as Figure 3 As shown, the maximum equivalent stress intensity factor is approximately 14.4 during the process from the initial crack to the formation of a through crack. The fracture toughness is much smaller than that of S135 steel, indicating that under the rated operating conditions of the turbine, the crown side crack on the blade outlet edge analyzed in this embodiment will not undergo brittle fracture, so there is no need to calculate the critical crack length for brittle fracture.

[0126] Curves, such as Figure 4 As shown, the maximum equivalent stress intensity factor varies by approximately 1.22 from the initial crack to the formation of a through crack. The value is still lower than the fatigue crack propagation threshold of S135 steel, indicating that under the rated operating conditions of the turbine, the crown side crack on the blade outlet edge analyzed in this embodiment will not exhibit fatigue crack propagation within any observed crack length range, and there is no possibility of fatigue fracture. Therefore, it is not necessary to calculate the critical crack length and fracture life.

[0127] In summary, the method described in this embodiment, through crack propagation simulation and fracture life analysis of the initial crown crack on the blade outlet edge of the mixed-flow turbine under rated power, assesses that the risk of fracture failure at the dynamic stress concentration point of the runner blade during the entire life cycle of the turbine is low under rated operating conditions. When formulating operation and maintenance plans, more attention can be paid to the fracture failure mode of the runner blade under non-rated operating conditions, such as transitional operating conditions, 50% power operation conditions, and load shedding conditions, thus providing a scientific and effective reference for predictive operation and maintenance.

[0128] Example 4 A crack propagation prediction system for a turbine runner under transient hydraulic excitation includes the following modules: The mesh generation module is used to draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, the rotating components including the runner blades; and to draw a horizontal concentrated field based on the cracks on the runner blades, the horizontal concentrated field including a normal horizontal concentrated field and a tangential horizontal concentrated field. The temporal displacement field generation module is used to determine the crack front based on the horizontal concentrated field, refine the mesh near the crack front to form the target mesh, filter the nodes of the target mesh, identify the reinforcing nodes, determine the type of reinforcing element where the reinforcing nodes are located, and solve the temporal displacement field of the runner under transient hydraulic excitation. The crack propagation direction and rate generation module is used to calculate the stress intensity factor and energy release rate of the crack front node based on the time-series displacement field of the runner, and to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation based on these. Crack propagation module: Used to update and correct the horizontal field based on the crack propagation direction and crack propagation rate to obtain the propagated crack.

[0129] The module is an object that uses physical or virtual representation to form an objective description of its morphological structure. The object is not the same as a physical object, and is not limited to physical or virtual. It can be a software program, electronic hardware, circuit module, or simulation object.

[0130] Finally, it should be noted that the above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit them. The scope of protection of the present invention is not limited thereto. Those skilled in the art should understand that any person skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features within the scope of the technology disclosed in the present invention, and these should all be covered within the scope of protection of the present invention.

Claims

1. A method for predicting crack propagation in a turbine runner under transient hydraulic excitation, characterized in that, Includes the following steps: Step 1: Draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, including the turbine blades; draw the horizontal concentrated field based on the cracks on the turbine blades, including the normal horizontal concentrated field and the tangential horizontal concentrated field; Step 2: Determine the crack front edge based on the horizontal field concentration, and refine the mesh near the crack front edge to form the target mesh; The nodes of the target mesh are screened to identify the reinforcement nodes, the type of reinforcement element to which the reinforcement nodes are located is determined, and the time-series displacement field of the runner under transient hydraulic excitation is solved. Step 3: Based on the time-series displacement field of the turbine runner, calculate the stress intensity factor and energy release rate of the crack front node, and use this as a basis to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation; Step 4: Based on the crack propagation direction and crack propagation rate obtained in Step 3, update and correct the horizontal field to obtain the propagated crack.

2. The method according to claim 1, characterized in that, It also includes step five: predicting the fracture life of the impeller; Step four further includes: determining whether a through crack has formed based on the expanded crack; if a through crack is formed, proceed to step five; otherwise, return to step two for iteration.

3. The method according to claim 1 or 2, characterized in that, In step two, the mesh near the crack tip is a tetrahedral mesh; the steps for refining the mesh near the crack tip to form the target mesh are as follows: Step A1: Determine the target mesh size based on the crack length, and then determine the number of mesh refinements by combining the original mesh size near the crack leading edge with the target mesh size; Step A2: Use the distance control index field to filter out grid cells that meet the encryption requirements; Step A3: Process the mesh cells that meet the encryption requirements, set the midpoints of each of their edges as newly added mesh points, and divide the original tetrahedral mesh cells into four new mesh cells. Step A4: Determine whether the encryption count has been reached. If yes, complete the encryption; otherwise, return to step A2.

4. The method according to claim 2, characterized in that, In step two, the method for filtering the nodes of the target mesh and identifying the enhanced nodes is as follows: Find the mesh elements in the target mesh that cross the zero isosurface at the crack front, and use the nodes of these mesh elements as candidate nodes; Iterate through all candidate nodes and filter out crack-enhanced nodes and crack tip-enhanced nodes; The crack reinforcement nodes are used to characterize the displacement discontinuity of the crack surface; The crack tip reinforcement node is used to characterize the stress singularity of the crack tip; Or / and in step two, the method for determining the type of enhancement unit where the enhancement node resides is as follows: A mesh element containing at least one crack reinforcement node is defined as a crack reinforcement element; a mesh element containing at least one crack tip reinforcement node is defined as a crack tip reinforcement element; and a mesh element containing both is defined as a hybrid reinforcement element.

5. The method according to claim 4, characterized in that, In step two, the method for solving the time-series displacement field of the turbine runner under transient hydraulic excitation is as follows: Based on the horizontal field, the reinforcement nodes near the crack tip are selected and the reinforcement elements in which they are located are determined; Based on the aforementioned enhancement unit, the unit stiffness matrix and mass matrix are calculated, and then assembled into the overall stiffness matrix and mass matrix of the impeller. Based on fluid dynamics simulation, data on pulsating fluid pressure loads under transient hydraulic excitation are obtained. Based on the overall stiffness matrix and mass matrix of the runner, as well as the pulsating fluid pressure load data, the acoustic-structure interaction dynamic equilibrium equation of the runner is solved to obtain the time-series displacement field of the runner.

6. The method according to claim 5, characterized in that, In step three, the steps for calculating the stress intensity factor and energy release rate at the crack front node based on the time-series displacement field of the turbine are as follows: Step B1: Obtain the location of the crack tip, take several nodes that are evenly distributed near the crack tip as the original calculation nodes, and establish a local coordinate system of the crack tip that matches the geometric direction of the crack tip. Step B2: Obtain the displacement field at each moment based on the time-series displacement field of the rotor, and extract the displacement results on both sides of the crack surface at each moment; Step B3: Transform the displacement results obtained in step B2 into the local coordinate system of the crack front, and then project them onto the original calculation nodes to form projected calculation nodes; Step B4: Based on the projection calculation nodes obtained in Step B3, draw a scatter plot and obtain the stress intensity factor of the crack front node at each time point through the scatter plot. Step B5: Based on the stress intensity factor obtained in step B4, the energy release rate of the crack front node at each time step is calculated using the Evan formula.

7. The method according to claim 6, characterized in that, In step three, the step of calculating the crack propagation direction under transient hydraulic excitation is as follows: Step C1: Based on the stress intensity factor of the crack leading edge node at each moment under transient hydraulic excitation, the maximum circumferential stress criterion is used to calculate the deflection angle of each crack leading edge node at each moment. The deflection angle lies within the plane formed by the tangential direction of the crack leading edge and the normal direction of the crack surface; Step C2: Collect the deflection angle data of all crack leading edge nodes at each time step; Step C3: Calculate the average deflection angle data collected in step C2 over time. The average deflection angle obtained is the crack propagation direction. Or / and in step three, the method for calculating the crack propagation rate under transient hydraulic excitation is as follows: Calculate the corresponding equivalent stress intensity factor using the energy release rate of the crack front node at each time point; The equivalent stress intensity factor at each time point is statistically analyzed across the entire range, and the amplitude of the equivalent stress intensity factor corresponding to all load cycles is extracted. Obtain the crack propagation threshold value of the rotating wheel; The amplitude of the equivalent stress intensity factor corresponding to each load cycle is judged. When it is greater than the crack propagation threshold, the crack propagation rate is calculated by the Paris formula. Otherwise, the assumed value is forced to be the crack propagation rate.

8. The method according to claim 7, characterized in that, Step four: Based on the crack propagation direction and crack propagation rate obtained in step three, update and correct the horizontal field to obtain the propagated crack. Step D1: Based on the tangential level set, divide all meshes into cracked and uncracked regions; For the already cracked area, the propagation rate of the normal horizontal set is fixed to zero, so that the existing crack surface remains unchanged. For the uncracked zone, the crack propagation rate is decomposed into normal propagation rate and tangential propagation rate based on the crack propagation direction obtained in step three. Step D2: Based on the normal propagation rate and the tangential propagation rate, update the horizontal field after crack propagation in the uncracked zone to characterize the updated crack surface and crack front topology, so as to obtain the expanded crack. Step D3: Correct the horizontal field obtained after crack propagation in step D2, ensuring that the absolute value of the new horizontal gradient field is always equal to 1, and that the normal horizontal field and the tangential horizontal field on the crack surface are orthogonal.

9. The method according to claim 7, characterized in that, Step five: The method for predicting the fracture life of the runner includes: Extract the crack length, maximum equivalent stress intensity factor, and maximum equivalent stress intensity factor amplitude in each iteration step; Plot the first relationship curve between the maximum equivalent stress intensity factor and the crack length. Based on this first relationship curve and the fracture toughness of the runner material, determine whether the runner blade will undergo brittle fracture. If there is a possibility of brittle fracture, determine the critical crack length for brittle fracture and provide monitoring prompts. A second relationship curve is plotted between the amplitude of the maximum equivalent stress intensity factor and the crack length. Based on this second relationship curve and the fatigue crack propagation threshold value of the impeller material, it is determined whether fatigue fracture will occur due to fatigue crack propagation. When fatigue fracture is possible, the critical crack length at which propagation begins is determined, and the time from the critical crack length at which propagation begins to the formation of a through crack is calculated. The calculation result is the fracture life of the impeller.

10. A crack propagation prediction system for a turbine runner under transient hydraulic excitation, characterized in that... Includes the following modules: The mesh generation module is used to draw the finite element mesh of the fluid domain and rotating components inside the turbine casing, the rotating components including the runner blades; and to draw a horizontal concentrated field based on the cracks on the runner blades, the horizontal concentrated field including a normal horizontal concentrated field and a tangential horizontal concentrated field. The temporal displacement field generation module is used to determine the crack front based on the horizontal concentrated field and to refine the mesh near the crack front to form the target mesh. The nodes of the target mesh are screened to identify the reinforcement nodes, the type of reinforcement element to which the reinforcement nodes are located is determined, and the time-series displacement field of the runner under transient hydraulic excitation is solved. The crack propagation direction and rate generation module is used to calculate the stress intensity factor and energy release rate of the crack front node based on the time-series displacement field of the runner, and to calculate the crack propagation direction and crack propagation rate under transient hydraulic excitation based on these. Crack propagation module: Used to update and correct the horizontal field based on the crack propagation direction and crack propagation rate to obtain the propagated crack.

Images