Anisotropic single crystal blade multi-scale dynamics calculation method considering fretting wear influence

Abstract

The invention relates to the technical field of blade vibration suppression and high-cycle fatigue prevention, and discloses an anisotropic single crystal blade multi-scale dynamics calculation method considering fretting wear influence, and the method specifically comprises the following steps: 1, building a finite element model of a single crystal blade under specific crystal orientation; through the method, the dynamic behavior of the single crystal blade under fretting wear can be accurately described, the fatigue life prediction precision of the blade is improved, through comprehensive consideration of the anisotropy, fretting wear and multi-scale effect of the single crystal blade, accurate simulation and analysis of the vibration characteristics of the blade can be realized, and the fatigue life prediction precision of the blade is improved. The method has the advantages that the method is simple and convenient to operate, an important basis is provided for blade design and optimization, scientific guidance can be provided for blade maintenance and replacement, economic loss and safety risks caused by blade failure are reduced, and the method can be popularized and applied to other materials and structures with similar characteristics and has important engineering application value.

Description

Technical Field

[0001] The present invention relates to the technical field of suppressing blade vibration and preventing high cycle fatigue, and in particular to a multi-scale dynamic calculation method for anisotropic single crystal blades taking into account the influence of micro-motion wear. Background Art

[0002] After years of development, dry friction damping technology has become one of the important technical means to suppress blade vibration and prevent high-cycle fatigue. Its principle is to convert the vibration energy into heat energy through the friction between the damper and the blade, and reduce the response amplitude of the blade through friction work. In this process, the contact surface between the blade and the damper will produce alternating normal stress and shear stress, and the material will be sheared under stress. The contact surface will then wear. Wear changes the profile of the contact interface, resulting in a weakened vibration reduction effect of the damper on the blade, and may even cause catastrophic failure of the blade.

[0003] However, due to technical defects, there is a deviation between the stacking axis of the actual single crystal blades and the growth axis of the crystal crystallization. In addition, single crystal materials have orthogonal anisotropy, and their basic physical properties such as elastic modulus, shear modulus and Poisson's ratio are directional. Therefore, the deflection of the crystal axis will significantly change the mechanical properties of the material, and thus affect the dynamic characteristics of the blade, which directly leads to changes in the nonlinear response considering wear.

[0004] To this end, we propose a multi-scale dynamic calculation method for anisotropic single crystal blades considering the influence of fretting wear. Summary of the invention

[0005] The present invention mainly solves the technical problems existing in the above-mentioned prior art and provides a multi-scale dynamic calculation method for anisotropic single crystal blades taking into account the influence of micro-motion wear.

[0006] In order to achieve the above object, the present invention adopts the following technical scheme, which is a multi-scale dynamic calculation method of anisotropic single crystal blade considering the influence of fretting wear, and specifically includes the following steps:

[0007] Step 1: Establish a finite element model of a single crystal blade under a specific crystal orientation;

[0008] Step 2: Use the multidimensional harmonic balance method (MHBM) to calculate the nonlinear dynamic response of the single crystal blade under different crystal orientation angles;

[0009] Step 3: Use the nonlinear response obtained in step 2 to solve the wear depth of the contact interface;

[0010] Step 4: First determine whether the wear depth has reached the maximum wear depth limit. If not, use the wear energy and wear depth calculated in step 3 to update the contact parameters and contact pressure of the contact interface.

[0011] Step 5: Using the contact parameters and contact preload calculated in step 4, repeat the process of steps 3 and 4 until the set maximum wear depth requirement is reached, and the wear nonlinear response and wear distribution of the contact interface under this crystal orientation angle can be obtained.

[0012] Preferably, the specific method of step 1 is as follows: the nickel-based single crystal high temperature alloy has a face-centered cubic structure (FCC), and its overall structure presents a unidirectional ordered lattice arrangement. This single crystal structure helps to improve the high temperature performance, corrosion resistance and mechanical properties of the alloy. Due to the special structure, the nickel-based single crystal alloy has orthogonal anisotropy, so its Young's modulus, shear modulus and Poisson's ratio in the XYZ direction are the same. The general strain-stress relationship of the orthogonal anisotropic material in the rectangular coordinate system can be expressed by Hooke's law:

[0013] ε=Dσ

[0014] Where ε and σ are the elastic strain tensor and elastic stress tensor, C is the flexibility matrix of the single crystal blade. For single crystal orthotropic materials, the elastic modulus in each direction is E, Poisson's ratio is μ, and the shear modulus is G. The elastic matrix D is:

[0015]

[0016] D 44 =D 55 =D 66 =G

[0017] Through the above-mentioned elastic-plastic constitutive relationship of nickel-based single crystal alloy, the dynamic equation of nickel-based single crystal blade can be written as follows:

[0018] ∑ρ∫∫∫ e [V] T [N]dv{δ(t)}+∑∫∫∫ e ([B] e ) T [D][B] e dv{δ(t)}=0

[0019] During the casting process of nickel-based turbine blades, the crystal growth direction [0,0,1] is along the stacking axis (Z axis) of the blade. However, during the actual blade growth process, the angle of crystallization inevitably produces a certain angle deviation from the initial stacking axis. Usually, the angle difference between the main axis of the crystal angle and the stacking axis is 0 to 15°, and is randomly distributed in a conical area with the stacking axis as the rotation axis. This angle is called the primary crystal orientation angle θ. In the plane perpendicular to the crystal growth direction [0,0,1], the growth direction of the crystal is defined as [0,1,0] or [1,0,0] orientation. The crystal growth direction in this plane rotates around the stacking axis of the blade. The obtained angle is defined as the secondary crystal angle orientation angle β. In order to analyze the influence of crystal orientation on the nonlinear response of the blade considering wear, the crystal orientation position of the single crystal blade must be accurately described. Therefore, the Euler angle rotation is used to establish the conversion relationship between the material coordinates and the sample coordinates to describe the crystal angle position. The Euler angle rotation is demonstrated through the common single crystal morphology structure of nickel-based single crystal alloy: first, the crystal rotates β around the Z axis of the sample coordinate system; the second time, it rotates θ around the axis of the coordinate system obtained by the first rotation; the third time, it rotates α around the z axis of the coordinate system obtained by the second rotation. In the Cartesian coordinate system, the rotation matrices about the coordinate axes X and Z are respectively:

[0020]

[0021] Where α represents any rotation angle. The rotation matrix R can be obtained by multiplying the rotation matrices obtained three times in sequence. The specific form is:

[0022]

[0023] Where l, m, n are the direction cosines between different coordinates, and the stiffness matrix transformation between the two coordinate systems is given by:

[0024] D * =TDT T

[0025] Where D* is the stiffness matrix after coordinate system transformation, and T is the stiffness transformation matrix, which is defined as:

[0026]

[0027] The stiffness matrix after crystal angle transformation can be calculated as:

[0028] [K] = ∑∫∫∫ e ([B] e ) T [D * ][B] e dv{δ(t)}

[0029] The crystal angle orientation of the single crystal blade will affect the wear behavior of the contact interface, and ultimately affect the nonlinear dynamic response after wear. The above method is used to express the dynamic equations of the single crystal blade under different crystal orientations.

[0030] Preferably, the specific method of step 2 is: the vibration differential equation of the dry friction system can be expressed as follows:

[0031]

[0032] Among them, M, C, and K are the system mass matrix, damping matrix, and stiffness matrix, respectively, with dimensions of Nd×Nd, Nd is the number of degrees of freedom, u is the system displacement, which is the unknown quantity to be solved, fnl represents the nonlinear force, and fsin(wt) is the external excitation force at the excitation frequency. The steady-state response displacement solution of the nonlinear system can be expressed as the following Fourier series:

[0033]

[0034] Where Nh is the truncated harmonic order, w is the external force excitation frequency, U0 is the Fourier series corresponding to the "0" order harmonic, Uck and Usk are the Fourier series corresponding to the k-order harmonic, and assuming that the system degree of freedom is Nd, then during the calculation process, the Nd(2Nh+1)-dimensional Fourier series vector is arranged as follows:

[0035]

[0036] At the same time, a displacement-coefficient conversion matrix is ​​defined as follows:

[0037]

[0038] Where I is the identity matrix of Nd×Nd dimensions, represents the Kronecker product, the displacement formula can be rewritten as follows:

[0039] u(t)=H(ωt i )

[0040] Introduce a frequency derivative operator:

[0041]

[0042] Then the velocity and acceleration can be expressed using the frequency derivative operator:

[0043]

[0044] Similarly, nonlinear forces and external excitation forces can also be expressed as the product of the transformation matrix and the corresponding Fourier coefficients:

[0045] f(t i )=H(ωt i )×[F 0 ,F c1 ,F s1 ,...F ck ,F sk ,...] T

[0046] =H(ωt i )F

[0047]

[0048] From this we can get:

[0049] G(ω,U)=P(ω)U+F nl (U)-F=0

[0050] P(ω)=ω 2 N M ▽ 2 +ωN C ▽+N K

[0051] Among them, P is the dynamic stiffness matrix, NM, NC, NK are block diagonal matrices composed of mass matrix M, mass matrix C, and stiffness matrix K, which are solved by the Newton-Raphson iterative method.

[0052] Preferably, the specific method of step three is: the essence of wear is energy conversion, and the energy consumed by a single hysteresis cycle during the vibration of the dry friction damper is:

[0053] E=Tdx

[0054] Among them, E represents the energy consumed in each hysteresis cycle, T represents the dry friction force of the contact node, and dx represents the relative displacement between the nodes at the selected analysis frequency. This value can be obtained by solving the frequency response function in step 2. The wear depth caused by each hysteresis cycle can be obtained by using the friction energy consumption, that is:

[0055]

[0056] Among them, α is defined as the energy correlation coefficient, which needs to be obtained in combination with experiments. Ex and Ey represent the energy consumed in two directions under dry friction, A represents the area of ​​the node contact surface, △h represents the wear depth of each vibration cycle, and the standard wear depth parameter νw is introduced. Usually νw is obtained through experimental tests, and different materials have different values. Combined with the maximum wear depth obtained in each vibration cycle, the wear acceleration coefficient can be obtained using the standard wear depth parameter:

[0057]

[0058] The calculated Zw is the wear acceleration coefficient. Since the wear depth of a single micro-motion wear is very small, and the wear depth of several micro-motion wears within a certain limit is almost the same, the product of the wear acceleration coefficient and the micro-motion wear depth is used to express the depth of several micro-motion wears to accelerate the wear process.

[0059] Preferably, the specific method of step 4 is as follows: the contact parameters include but are not limited to contact stiffness and friction coefficient. During the wear process, factors that usually cause surface cracks or surface damage will cause the morphology of the contact surface to undergo irreversible changes as the wear process proceeds, thereby causing changes in the tangential stiffness and friction coefficient of the contact surface. According to simulation research and test results, the following definitions are made for the changes in the tangential stiffness and friction coefficient with the friction dissipation energy:

[0060]

[0061] μ=1+b μ log(a μ E+1)

[0062] Among them, akt and bkt represent the stiffness variation coefficients, aμ and bμ represent the wear factor variation coefficients. As the wear of the dry friction contact surface evolves, the tangential stiffness, friction coefficient and contact pressure between the contact nodes will change with the advancement of the wear process and the loss of wear energy, thereby affecting the relative friction motion state of the contact surface. As the wear process proceeds, the most intuitive change is the continuous increase in the wear depth of the contact surface. The generation of wear depth will directly affect the pre-pressure between the contact surfaces, and then directly affect the motion state between the contact nodes in the nonlinear response. The change of the contact surface pre-pressure with the wear depth is quantitatively described as:

[0063]

[0064] Preferably, in addition to limiting the maximum number of iterations by the maximum wear depth, the step five may also control the number of iterations by limiting the maximum number of wear cycles.

[0065] The present invention provides a multi-scale dynamic calculation method for anisotropic single crystal blades taking into account the influence of fretting wear. It has the following beneficial effects:

[0066] 1. The multi-scale dynamic calculation method of anisotropic single crystal blades considering the influence of micro-motion wear can accurately describe the dynamic behavior of single crystal blades under micro-motion wear, and improve the fatigue life prediction accuracy of blades. By comprehensively considering the anisotropy, micro-motion wear and multi-scale effects of single crystal blades, the method of the present invention can realize accurate simulation and analysis of blade vibration characteristics, and provide an important basis for blade design and optimization. At the same time, the method can also provide scientific guidance for the maintenance and replacement of blades, and reduce the economic losses and safety risks caused by blade failure. In addition, the method of the present invention can also be extended to other materials and structures with similar characteristics, and has important engineering application value.

[0067] 2. The multi-scale dynamic calculation method of anisotropic single crystal blades considering the influence of micro-wear can further analyze the wear mechanism and the specific impact of wear on blade performance through the wear nonlinear response and wear distribution of the contact interface obtained in step 4, combined with the microstructure characteristics of the single crystal blade. At the same time, this method can also consider factors such as load conditions, temperature environment and lubrication status under different working conditions, so as to more comprehensively evaluate the durability and reliability of single crystal blades in actual operation. In addition, this method can also be combined with other multidisciplinary optimization methods, such as structural optimization design, material performance optimization, etc., to achieve the optimization of the overall performance of single crystal blades.

[0068] 3. The multi-scale dynamic calculation method of anisotropic single crystal blades considering the influence of micro-motion wear can reveal the deep influence of micro-motion wear on the dynamic characteristics of blades by fine-tuning the wear process of single crystal blades. This method can capture the slight changes in the contact interface during the wear process, such as the evolution of contact stiffness and friction coefficient, and how these changes further affect the vibration mode and stress distribution of the blades. This fine-tuned simulation not only helps to deeply understand the wear mechanism, but also provides more precise guidance for the maintenance and care of blades. By predicting the dynamic response of blades under different degrees of wear, a more scientific maintenance plan can be formulated, thereby avoiding unnecessary downtime and economic losses. In addition, this method can also provide strong support for the research and development of new single crystal blade materials. By simulating the performance of different materials under micro-motion wear, materials with better durability and reliability can be screened out, providing a scientific basis for the replacement of blades. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 is a flow chart of the method of the present invention;

[0070] Figure 2 A schematic diagram of a single crystal blade model is used for the calculation examples of the present invention;

[0071] Figure 3 It is a schematic diagram showing the change of the response of the present invention with the number of iterations;

[0072] Figure 4 It is a schematic diagram of the wear distribution prediction of the contact interface of the present invention;

[0073] Figure 5 It is a schematic diagram of the response before and after wear under different crystal orientations of the present invention. DETAILED DESCRIPTION

[0074] Embodiment 1: Multi-scale dynamic calculation method of anisotropic single crystal blade considering the influence of fretting wear, such as Figure 1-Figure 5 As shown, the specific steps include: Step 1: Establish a finite element model of a single crystal blade under a specific crystal orientation; Step 2: Use the multidimensional harmonic balance method (MHBM) to calculate the nonlinear dynamic response of the single crystal blade under different crystal orientation angles; Step 3: Use the nonlinear response obtained in Step 2 to solve the wear depth of the contact interface; Step 4: First determine whether the wear depth has reached the maximum wear depth limit. If it has not reached, use the wear energy and wear depth calculated in Step 3 to update the contact parameters and contact pressure of the contact interface; Step 5: Use the contact parameters and contact prepressure calculated in Step 4, repeat Step 3 and Step 4 until the set maximum wear depth requirement is reached, and the wear nonlinear response and wear distribution of the contact interface under this crystal orientation angle can be obtained. In addition to limiting the maximum number of iterations by the maximum wear depth, the number of iterations can also be controlled by limiting the maximum number of wear cycles. This method can accurately describe the dynamic behavior of single-crystal blades under micro-motion wear, and improve the accuracy of fatigue life prediction of blades. By comprehensively considering the anisotropy, micro-motion wear and multi-scale effects of single-crystal blades, the method of the present invention can achieve accurate simulation and analysis of blade vibration characteristics, providing an important basis for blade design and optimization. At the same time, the method can also provide scientific guidance for blade maintenance and replacement, reducing economic losses and safety risks caused by blade failure. In addition, the method of the present invention can also be extended to other materials and structures with similar properties, and has important engineering application value.

[0075] Embodiment 2: Based on Embodiment 1, Figure 2 As shown, the specific method of step one is: the nickel-based single crystal high-temperature alloy has a face-centered cubic structure (FCC), and its overall structure presents a unidirectional ordered lattice arrangement. This single crystal structure helps to improve the high temperature performance, corrosion resistance and mechanical properties of the alloy. Due to the special structure, the nickel-based single crystal alloy has orthogonal anisotropy, so its Young's modulus, shear modulus and Poisson's ratio in the XYZ direction are the same. The general strain-stress relationship of the orthogonal anisotropic material in the rectangular coordinate system can be expressed by Hooke's law:

[0076] ε=Dσ

[0077] Where ε and σ are the elastic strain tensor and elastic stress tensor, C is the flexibility matrix of the single crystal blade. For single crystal orthotropic materials, the elastic modulus in each direction is E, Poisson's ratio is μ, and the shear modulus is G. The elastic matrix D is:

[0078]

[0079] D 44 =D 55 =D 66 =G

[0080] Through the above-mentioned elastic-plastic constitutive relationship of nickel-based single crystal alloy, the dynamic equation of nickel-based single crystal blade can be written as follows:

[0081] ∑ρ∫∫∫ e [V] T [N]dv{δ(t)}+∑∫∫∫ e ([B] e ) T [D][B] e dv{δ(t)}=0

[0082] During the casting process of nickel-based turbine blades, the crystal growth direction [0,0,1] is along the stacking axis (Z axis) of the blade. However, during the actual blade growth process, the angle of crystallization inevitably produces a certain angle deviation from the initial stacking axis. Usually, the angle difference between the main axis of the crystal angle and the stacking axis is 0 to 15°, and is randomly distributed in a conical area with the stacking axis as the rotation axis. This angle is called the primary crystal orientation angle θ. In the plane perpendicular to the crystal growth direction [0,0,1], the growth direction of the crystal is defined as [0,1,0] or [1,0,0] orientation. The crystal growth direction in this plane rotates around the stacking axis of the blade. The obtained angle is defined as the secondary crystal angle orientation angle β. In order to analyze the influence of crystal orientation on the nonlinear response of the blade considering wear, the crystal orientation position of the single crystal blade must be accurately described. Therefore, the Euler angle rotation is used to establish the conversion relationship between the material coordinates and the sample coordinates to describe the crystal angle position. The Euler angle rotation is demonstrated through the common single crystal morphology structure of nickel-based single crystal alloy: first, the crystal rotates β around the Z axis of the sample coordinate system; the second time, it rotates θ around the axis of the coordinate system obtained by the first rotation; the third time, it rotates α around the z axis of the coordinate system obtained by the second rotation. In the Cartesian coordinate system, the rotation matrices about the coordinate axes X and Z are respectively:

[0083]

[0084] Where α represents any rotation angle. The rotation matrix R can be obtained by multiplying the rotation matrices obtained three times in sequence. The specific form is:

[0085]

[0086] Where l, m, n are the direction cosines between different coordinates, and the stiffness matrix transformation between the two coordinate systems is given by:

[0087] D * =TDT T

[0088] Where D* is the stiffness matrix after coordinate system transformation, and T is the stiffness transformation matrix, which is defined as:

[0089]

[0090] The stiffness matrix after crystal angle transformation can be calculated as:

[0091] [K] = ∑∫∫∫ e ([B] e ) T [D * ][B] e dv{δ(t)}

[0092] The crystal angle orientation of the single crystal blade will affect the wear behavior of the contact interface, and ultimately affect the nonlinear dynamic response after wear. The above method is used to express the dynamic equation of the single crystal blade under different crystal orientations. Through the wear nonlinear response obtained in step 4 and the wear distribution of the contact interface, combined with the microstructure characteristics of the single crystal blade, the wear mechanism and the specific impact of wear on the blade performance can be further analyzed. At the same time, this method can also consider factors such as load conditions, temperature environment and lubrication status under different working conditions, so as to more comprehensively evaluate the durability and reliability of the single crystal blade in actual operation. In addition, this method can also be combined with other multidisciplinary optimization methods, such as structural optimization design, material performance optimization, etc., to achieve the optimization of the overall performance of the single crystal blade.

[0093] Embodiment 3: Based on Embodiment 1 and Embodiment 2, Figure 3 As shown, the specific method of step 2 is: the vibration differential equation of the dry friction system can be expressed as follows:

[0094]

[0095] Among them, M, C, and K are the system mass matrix, damping matrix, and stiffness matrix, respectively, with dimensions of Nd×Nd, Nd is the number of degrees of freedom, u is the system displacement, which is the unknown quantity to be solved, fnl represents the nonlinear force, and fsin(wt) is the external excitation force at the excitation frequency. The steady-state response displacement solution of the nonlinear system can be expressed as the following Fourier series:

[0096]

[0097] Where Nh is the truncated harmonic order, w is the external force excitation frequency, U0 is the Fourier series corresponding to the "0" order harmonic, Uck and Usk are the Fourier series corresponding to the k-order harmonic, and assuming that the system degree of freedom is Nd, then during the calculation process, the Nd(2Nh+1)-dimensional Fourier series vector is arranged as follows:

[0098]

[0099] At the same time, a displacement-coefficient conversion matrix is ​​defined as follows:

[0100]

[0101] Where I is the identity matrix of Nd×Nd dimensions, represents the Kronecker product, the displacement formula can be rewritten as follows:

[0102] u(t)=H(ωt i )

[0103] Introduce a frequency derivative operator:

[0104]

[0105] Then the velocity and acceleration can be expressed using the frequency derivative operator:

[0106]

[0107] Similarly, nonlinear forces and external excitation forces can also be expressed as the product of the transformation matrix and the corresponding Fourier coefficients:

[0108] f(t i )=H(ωt i )×[F 0 ,F c1 ,F s1 ,...F ck ,F sk ,...] T

[0109] =H(ωt i )F

[0110]

[0111] From this we can get:

[0112] G(ω,U)=P(ω)U+F nl (U)-F=0

[0113] P(ω)=ω 2 NM ▽ 2 +ωN C ▽+N K

[0114] Among them, P is the dynamic stiffness matrix, NM, NC, NK are block diagonal matrices composed of mass matrix M, mass matrix C, and stiffness matrix K, which are solved by the Newton-Raphson method iteration method. Through the refined simulation of the wear process of single-crystal blades, the deep-seated influence of micro-motion wear on the dynamic characteristics of blades can be revealed. This method can capture the slight changes in the contact interface during the wear process, such as the evolution of contact stiffness and friction coefficient, and how these changes further affect the vibration mode and stress distribution of the blades. This refined simulation not only helps to deeply understand the wear mechanism, but also provides more precise guidance for the maintenance and maintenance of blades. By predicting the dynamic response of blades under different degrees of wear, a more scientific maintenance plan can be formulated to avoid unnecessary downtime and economic losses. In addition, this method can also provide strong support for the research and development of new single-crystal blade materials. By simulating the performance of different materials under micro-motion wear, materials with better durability and reliability can be screened out, providing a scientific basis for the replacement of blades.

[0115] Embodiment 4: Based on Embodiment 1, Embodiment 2 and Embodiment 3, Figure 4 As shown, the specific method of step three is: the essence of wear is energy conversion, and the energy consumed by a single hysteresis cycle during the vibration of the dry friction damper is:

[0116] E=Tdx

[0117] Among them, E represents the energy consumed in each hysteresis cycle, T represents the dry friction force of the contact node, and dx represents the relative displacement between the nodes at the selected analysis frequency. This value can be obtained by solving the frequency response function in step 2. The wear depth caused by each hysteresis cycle can be obtained by using the friction energy consumption, that is:

[0118]

[0119] Among them, α is defined as the energy correlation coefficient, which needs to be obtained in combination with experiments. Ex and Ey represent the energy consumed in two directions under dry friction, A represents the area of ​​the node contact surface, △h represents the wear depth of each vibration cycle, and the standard wear depth parameter νw is introduced. Usually νw is obtained through experimental tests, and different materials have different values. Combined with the maximum wear depth obtained in each vibration cycle, the wear acceleration coefficient can be obtained using the standard wear depth parameter:

[0120]

[0121] The calculated Zw is the wear acceleration coefficient. Since the wear depth of a single micro-motion wear is very small, and the wear depth of several micro-motion wears within a certain limit is almost the same, the product of the wear acceleration coefficient and the micro-motion wear depth is used to express the depth of several micro-motion wears to accelerate the wear process.

[0122] Embodiment 5: Based on Embodiment 1, Embodiment 2, Embodiment 3 and Embodiment 4, Figure 5 As shown, the specific method of step 4 is: contact parameters include but are not limited to contact stiffness and friction coefficient. During the wear process, factors that usually cause surface cracks or surface damage will cause the morphology of the contact surface to undergo irreversible changes as the wear process progresses, which in turn causes changes in the tangential stiffness and friction coefficient of the contact surface. Based on simulation research and test results, the following definitions are made for the changes in tangential stiffness and friction coefficient with friction dissipation energy:

[0123]

[0124] μ=1+b μ log(a μ E+1)

[0125] Among them, akt and bkt represent the stiffness variation coefficients, aμ and bμ represent the wear factor variation coefficients. As the wear of the dry friction contact surface evolves, the tangential stiffness, friction coefficient and contact pressure between the contact nodes will change with the advancement of the wear process and the loss of wear energy, thereby affecting the relative friction motion state of the contact surface. As the wear process proceeds, the most intuitive change is the continuous increase in the wear depth of the contact surface. The generation of wear depth will directly affect the pre-pressure between the contact surfaces, and then directly affect the motion state between the contact nodes in the nonlinear response. The change of the contact surface pre-pressure with the wear depth is quantitatively described as:

[0126]

[0127] Working principle of the present invention: The present invention uses Figure 2 The nickel-based single crystal shrouded blade shown in the figure is calculated. In this model, the shroud parts are in contact with each other, causing wear. Therefore, a friction contact pair is established between adjacent shrouds. The root part of this blade is fixed. In actual working conditions, the blade is excited by the airflow, resulting in a larger amplitude near the tip. Therefore, the excitation point is arranged at the position of the blade close to the shroud. At the same time, in order to simulate the preload of adjacent blades during assembly, two static loads are applied to both sides of the shroud of the blade. First, a shrouded blade model with an ideal crystal orientation is calculated separately. Using the method proposed in the present invention, its maximum wear depth is controlled to be 0.1 mm. The response result after calculation is as follows Figure 3 The wear distribution of the contact interface is shown in Figure 4As shown, for different crystal angle orientations, the maximum wear cycle is controlled to be 2 million times, and the method proposed in the present invention is used for calculation iteration. The response before and after wear is as follows Figure 5 shown.

[0128] The above shows and describes the basic principles and main features of the present invention and the advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments. The above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention to be protected. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A multi-scale dynamic calculation method for anisotropic single crystal blades considering the influence of fretting wear, characterized in that: The specific steps include: Step 1: Establish a finite element model of a single crystal blade under a specific crystal orientation; Step 2: Use the multidimensional harmonic balance method (MHBM) to calculate the nonlinear dynamic response of the single crystal blade under different crystal orientation angles; Step 3: Use the nonlinear response obtained in step 2 to solve the wear depth of the contact interface; Step 4: First determine whether the wear depth has reached the maximum wear depth limit. If not, use the wear energy and wear depth calculated in step 3 to update the contact parameters and contact pressure of the contact interface. Step 5: Using the contact parameters and contact preload calculated in step 4, repeat the process of steps 3 and 4 until the set maximum wear depth requirement is reached, and the wear nonlinear response and wear distribution of the contact interface under this crystal orientation angle can be obtained.

2. The multi-scale dynamic calculation method for anisotropic single crystal blades considering the influence of fretting wear according to claim 1 is characterized in that: The specific method of step 1 is as follows: the nickel-based single crystal high temperature alloy has a face-centered cubic structure (FCC), and its overall structure presents a unidirectional ordered lattice arrangement. This single crystal structure helps to improve the high temperature performance, corrosion resistance and mechanical properties of the alloy. Due to the special structure, the nickel-based single crystal alloy has orthogonal anisotropy, so its Young's modulus, shear modulus and Poisson's ratio in the XYZ direction are the same. The general strain-stress relationship of the orthogonal anisotropic material in the rectangular coordinate system can be expressed by Hooke's law: ε=Dσ Where ε and σ are the elastic strain tensor and elastic stress tensor, C is the flexibility matrix of the single crystal blade. For single crystal orthotropic materials, the elastic modulus in each direction is E, Poisson's ratio is μ, and the shear modulus is G. The elastic matrix D is: D 44 =D 55 =D 66 =G Through the above-mentioned elastic-plastic constitutive relationship of nickel-based single crystal alloy, the dynamic equation of nickel-based single crystal blade can be written as follows: ∑ρ∫∫∫ e [V] T [N]dv{δ(t)}+∑∫∫∫ e ([B] e ) T [D][B] e dv{δ(t)}=0 During the casting process of nickel-based turbine blades, the crystal growth direction [0,0,1] is along the stacking axis (Z axis) of the blade. However, during the actual blade growth process, the angle of crystal crystallization inevitably produces a certain angle deviation from the initial stacking axis. Usually, the angle difference between the main axis of the crystal angle and the stacking axis is 0-15°, and is randomly distributed in a conical area with the stacking axis as the rotation axis. This angle is called the primary crystal orientation angle θ. In the plane perpendicular to the crystal growth direction [0,0,1], the crystal growth direction is defined as [0,1,0] or [1,0,0] orientation. The angle obtained by rotating the crystal growth direction of this plane around the stacking axis of the blade is defined as the secondary crystal angle orientation angle β. In order to analyze the influence of crystal orientation on the nonlinear response of the blade considering wear, the crystal orientation position of the single crystal blade must be accurately described. Therefore, the Euler angle rotation is used to establish the conversion relationship between the material coordinates and the sample coordinates to describe the crystal angle position. The Euler angle rotation is demonstrated through the common single crystal morphology structure of nickel-based single crystal alloy: first, the crystal is rotated around the Z axis of the sample coordinate system by β; The second time, the axis of the coordinate system obtained by the first rotation is rotated by θ; the third time, the axis of the coordinate system obtained by the second rotation is rotated by α. In the Cartesian coordinate system, the rotation matrices about the coordinate axes X and Z are respectively: Where α represents any rotation angle. The rotation matrix R can be obtained by multiplying the rotation matrices obtained three times in sequence. The specific form is: Where l, m, n are the direction cosines between different coordinates, and the stiffness matrix transformation between the two coordinate systems is given by: D * =TDT T Where D* is the stiffness matrix after coordinate system transformation, and T is the stiffness transformation matrix, which is defined as: The stiffness matrix after crystal angle transformation can be calculated as: [K]=∑∫∫∫ e ([B] e ) T [D * ][B] e dv{δ(t)} The crystal angle orientation of the single crystal blade will affect the wear behavior of the contact interface, and ultimately affect the nonlinear dynamic response after wear. The above method is used to express the dynamic equations of the single crystal blade under different crystal orientations.

3. The multi-scale dynamic calculation method of anisotropic single crystal blade considering the influence of fretting wear according to claim 1 is characterized in that: The specific method of step 2 is: the vibration differential equation of the dry friction system can be expressed as follows: Among them, M, C, and K are the system mass matrix, damping matrix, and stiffness matrix, respectively, with dimensions of Nd×Nd, Nd is the number of degrees of freedom, u is the system displacement, which is the unknown quantity to be solved, fnl represents the nonlinear force, and fsin(wt) is the external excitation force at the excitation frequency. The steady-state response displacement solution of the nonlinear system can be expressed as the following Fourier series: Where Nh is the truncated harmonic order, w is the external excitation frequency, U0 is the Fourier series corresponding to the "0" order harmonic, Uck and Usk are the Fourier series corresponding to the k-order harmonic, and assuming that the system degree of freedom is Nd, then during the calculation process, the Nd(2Nh+1)-dimensional Fourier series vector is arranged as follows: At the same time, a displacement-coefficient conversion matrix is ​​defined as follows: Where I is the identity matrix of Nd×Nd dimensions, represents the Kronecker product, the displacement formula can be rewritten as follows: u(t)=H(ωt i )U Introduce a frequency derivative operator: Then the velocity and acceleration can be expressed using the frequency derivative operator: Similarly, nonlinear forces and external excitation forces can also be expressed as the product of the transformation matrix and the corresponding Fourier coefficients: f(t i )=H(ωt i )×[F 0 ,F c1 ,F s1 ,...F ck ,F sk ,...] T =H(ωt i )F From this we can get: G(ω,U)=P(ω)U+F nl (U)-F=0 P(ω)=ω 2 N M ▽ 2 +ωN C ▽+N K Among them, P is the dynamic stiffness matrix, NM, NC, NK are block diagonal matrices composed of mass matrix M, mass matrix C, and stiffness matrix K, which are solved by the Newton-Raphson iterative method.

4. The multi-scale dynamic calculation method of anisotropic single crystal blade considering the influence of fretting wear according to claim 1 is characterized in that: The specific method of step three is: the essence of wear is energy conversion, and the energy consumed by a single hysteresis cycle during the vibration of the dry friction damper is: E=Tdx Among them, E represents the energy consumed in each hysteresis cycle, T represents the dry friction force of the contact node, and dx represents the relative displacement between the nodes at the selected analysis frequency. This value can be obtained by solving the frequency response function in step 2. The wear depth caused by each hysteresis cycle can be obtained by using the friction energy consumption, that is: Among them, α is defined as the energy correlation coefficient, which needs to be obtained in combination with experiments. Ex and Ey represent the energy consumed in two directions under dry friction, A represents the area of ​​the node contact surface, △h represents the wear depth of each vibration cycle, and the standard wear depth parameter νw is introduced. Usually νw is obtained through experimental tests, and different materials have different values. Combined with the maximum wear depth obtained in each vibration cycle, the wear acceleration coefficient can be obtained using the standard wear depth parameter: The calculated Zw is the wear acceleration coefficient. Since the wear depth of a single micro-motion wear is very small, and the wear depth of several micro-motion wears within a certain limit is almost the same, the product of the wear acceleration coefficient and the micro-motion wear depth is used to express the depth of several micro-motion wears to accelerate the wear process.

5. The multi-scale dynamic calculation method of anisotropic single crystal blade considering the influence of fretting wear according to claim 1 is characterized in that: The specific method of step 4 is as follows: contact parameters include but are not limited to contact stiffness and friction coefficient. During the wear process, factors that usually cause surface cracks or surface damage will cause the morphology of the contact surface to undergo irreversible changes as the wear process progresses, thereby causing changes in the tangential stiffness and friction coefficient of the contact surface. Based on simulation research and test results, the following definitions are made for the changes in tangential stiffness and friction coefficient with friction dissipation energy: µ1+b μ log ( a ) μ E+1) Among them, akt and bkt represent the stiffness variation coefficients, aμ and bμ represent the wear factor variation coefficients. As the wear of the dry friction contact surface evolves, the tangential stiffness, friction coefficient and contact pressure between the contact nodes will change with the advancement of the wear process and the loss of wear energy, thereby affecting the relative friction motion state of the contact surface. As the wear process proceeds, the most intuitive change is the continuous increase in the wear depth of the contact surface. The generation of wear depth will directly affect the pre-pressure between the contact surfaces, and then directly affect the motion state between the contact nodes in the nonlinear response. The change of the contact surface pre-pressure with the wear depth is quantitatively described as:

6. The multi-scale dynamic calculation method of anisotropic single crystal blade considering the influence of fretting wear according to claim 1 is characterized in that: In addition to limiting the maximum number of iterations by the maximum wear depth, the step five may also control the number of iterations by limiting the maximum number of wear cycles.

Images