Dynamic response rapid prediction method under complex excitation of double-variable cycle engine

Abstract

The invention discloses a dynamic response rapid prediction method under complex excitation of a double-variable cycle engine. The dynamic response rapid prediction method comprises the following steps: calculating a steady force on the surface of a turbine blade disc in a mode switching process; according to the dynamic response rapid prediction method under the complex excitation of the double-variable cycle engine, a modal compression method is used for achieving reduction of model dimensions, and the calculation model obtained through the method can greatly reduce resources needed by calculation and improve the calculation efficiency. Compared with an existing full-size model calculation method, the calculation efficiency is improved, the model scale is remarkably reduced, errors are extremely small, and the model scale directly affects the solving efficiency of the finite element model. Besides, in the previous research, an analysis method for the dynamic response of the bladed disc under mode switching of the double-variable cycle engine is not provided yet, and the blank of the response prediction technology of the bladed disc under sudden change excitation is filled.

Description

Technical Field

[0001] The present invention relates to the technical field of dual variable cycle engines, and specifically to a method for rapidly predicting the dynamic response of a dual variable cycle engine under complex excitations. Background Art

[0002] Due to the sudden change in aerodynamic load caused by the switching of the working mode during wartime and the wider load coverage frequency of the turbine disk of the dual variable cycle engine, higher requirements are put forward for the structural design of the turbine disk and its vibration suppression. Therefore, developing the coupled dynamic analysis technology of the turbine disk under the condition of sudden load change to ensure that the whole machine rotor system has sufficient safety margin is one of the key objectives in the design stage of the variable cycle engine;

[0003] To realize the simulation analysis of the dynamic characteristics of the turbine disk of the dual variable cycle engine, it is first necessary to solve the problem of characterizing the unsteady aerodynamic load of the turbine disk under the service conditions of a complex flow field environment, so as to accurately simulate the characteristics such as the sudden change of blade surface excitation induced by the switching of the working mode, and provide input for the subsequent response calculation.

[0004] At present, no calculation method for predicting the response under complex loads of a dual variable cycle engine has been developed, and general response prediction methods cannot meet the needs of such sudden loads. In addition, since the calculation process of the dynamic response of the turbine disk will be very time-consuming; therefore, an efficient method for predicting the disk response under sudden loads needs to be proposed to achieve an accurate analysis of the dynamic response of the turbine disk of the dual variable cycle engine. Summary of the Invention

[0005] In view of the problems existing in the prior art, the present invention discloses a method for rapidly predicting the dynamic response of a dual variable cycle engine under complex excitations. The technical solution adopted is as follows: Step 1: Calculation of the steady force on the surface of the turbine disk during the mode switching process; establish an air flow load model as the input for response calculation, establish an aerodynamic load model through CFD numerical simulation, follow the mass conservation equation, momentum conservation equation, and energy conservation equation, integrate the control volume, obtain a discrete equation set, and find the dependent variables of the nodes. The discrete equation set is as follows: Among them, ρ is density, p is momentum, t is time, U is the velocity of the fluid, τ is the stress tensor, λ is the thermal conductivity of the fluid, h tol is the total enthalpy, S M and S E are source terms; In the Cartesian coordinate system, the integral form of the control volume equation is as follows: where \(V\) is the volume of the control volume, \(S\) is the surface area of the control volume, and \(\mathrm{d}\vec{n}\) j is the unit normal vector perpendicular to the control volume surface; when performing unsteady calculations, a high-order solution mode is adopted for the solution, the second-order backward Euler format is selected for the time term, and the standard \(k-\omega\) turbulence model is adopted for the turbulence model, which includes the turbulent kinetic energy \(k\) equation and the turbulent frequency \(\omega\) equation: \(\rho\) is the density, \(t\) is the time, \(U\) is the velocity component, \(x\) is the spatial coordinate component, \(\sigma\) is the turbulent Prandtl number, \(\beta\) is the turbulent kinetic energy dissipation coefficient, \(\mu\) is the turbulent dynamic viscosity, and \(\mu\) t is the turbulent viscosity, and \(P\) is the potential energy; Step 2: Steady force TSM model on the blisk surface; the external excitation distribution in the cylindrical coordinate system is described by the Tyler - Sofrin excitation model (TSMs), the unsteady aerodynamic force is decomposed in time, and the harmonic components at a single engine order (EO) are extracted through the following formula: The complex Fourier coefficients \(P(h)\) extracted at the blade surface grid points are further spatially decomposed circumferentially to determine the rotational perturbations that will excite the corresponding traveling wave modes; the time - space Fourier coefficients are extracted at the corresponding points along the grid surface as follows: where \(m\) is the circumferential order, and \(P\) h m is the tangential distribution of the discrete - time Fourier coefficients at the corresponding surface points of the blade, \(N\) is the number of blades. In the current spatial Fourier decomposition, the circumferential DFT is performed on a sampling set that only takes \(N\) samples on each blade at the grid points corresponding to the surface positions. Therefore, according to the Nyquist theorem, when the number of blades is odd, the extracted rotational pressure components have multiple traveling waves from \(-N / 2\) to \(+N / 2\) or from \(-(N - 1) / 2\) to \(+(N - 1) / 2\). This method also takes into account the aliasing phenomenon experienced by the blade row when the traveling wave excitation of the excitation order is higher than \(N / 2\); By equating the unsteady aerodynamic force data of all nodes to the excitation of a few key nodes, these key nodes, according to the equivalence principle, are equivalent to the unsteady aerodynamic forces of each node obtained by computational fluid dynamics based on the dynamic response. Reducing them to a very small number of excitation nodes can greatly reduce the dimension of the computational model. During the equivalence process, first calculate the airflow excitation response at this rotational speed, then calculate the initial response on the simplified nodes, and judge whether the responses are equivalent according to the root mean square error between the response of the simplified nodes and the airflow excitation response; if not, update and iterate the numerical value of the nodal force according to the sensitivity coefficient of the nodal force to the response function until the response function meets the requirements. During the equivalence process, calculate the sensitivity matrix by the finite difference method, perturb the excitation force and calculate the response to obtain the sensitivity: Among them, R(F) is the response solution function, F is the excitation force. After forming the sensitivity matrix S through sensitivity, the update formula for the nodal force amplitude is based on the least squares method, using the pseudo-inverse of the sensitivity matrix S to minimize the response error: Among them, k is the error iteration number. For the dual variable cycle model, select three nodes with different chord lengths to equivalent the excitation force. Step 3: The main node dimensionality reduction model of the turbine disk and blade. Calculate the dynamic response of the disk and blade system by the model dimensionality reduction method of the main nodes. During the calculation process, select the equivalent excitation nodes, response nodes, and contact interface nodes on the surface as the main nodes. The modal dimensionality reduction method based on the main nodes is more suitable for the application scenarios when there are many contact nodes and non-linear calculations are carried out; for the multi-harmonic motion equation, only retain the "non-linear" degrees of freedom, that is, the degrees of freedom applying non-linear forces, and reduce all other "linear" degrees of freedom. The calculation can be completed by simplifying the model. During this process, the conversion from the physical domain to the modal domain is realized through the FRF matrix, and the dimensionality reduction of the model can be achieved. The FRF matrix of the multi-harmonic motion equation is as follows: Q j =A(m j ω)(P j -F j (Q)) Q j is the overall harmonic expansion matrix, A is the overall frequency response function matrix, P j is the dynamic stiffness matrix, F j is the amplitude of the excitation force in the frequency domain, Q is the harmonic component, and then divide the non-linear harmonic component Q n and the linear harmonic component Q l , and for each harmonic, it can be written in the following form: A ll is the linear frequency response function matrix, Anl With A ln is the frequency response function matrix of linear and nonlinear coupling, A nn Frequency response function FRF matrix of the nonlinear term, Q n is the harmonic component of the nonlinear term, and the equation can be simplified to the following form: where is the modal displacement vector of the nonlinear degrees of freedom calculated without nonlinear contact interaction forces; using the modal characteristics of the linear structure, the FRF matrix is usually approximated by its modal expansion where, A nn is the processed frequency response function FRF matrix of the nonlinear term, Φ j is the modal shape, ω is the frequency of the calculated order, N is the number of retained modal orders, ω j is the frequency at this order, and η is the modal damping ratio; Step 4: The reduced-order model of the main nodes of the turbine disk. When the dual variable cycle engine switches from the turbojet to the turbofan mode, the disk system generally goes through three stages in chronological order, namely the turbojet windmill stage, the instantaneous mode switching stage, and the turbofan mode stable stage. During the mode switching process, the change in the airflow excitation amplitude is closely related to the rotational speed change; under the high-pressure airflow excitation of the main combustion chamber, the amplitude of the airflow will gradually increase and be proportional to the rotational speed as the disk speed increases; with the continuous action of the high-pressure airflow, the rotational speed of the disk continues to rise, and the excitation amplitude of the airflow will also increase accordingly; at this time, the air excitation amplitude increases with the increase in rotational speed; until the high-pressure turbine disk reaches the turbofan mode speed and the rotational speed change stops, the airflow excitation amplitude tends to be stable, and the dynamic equation of the blade is as follows: And the disk will be hardened by the centrifugal force during rotation, making the overall stiffness of the disk show an increasing trend with the rotational speed. At this time, the dynamic equation can be expressed as: where, Ω is the disk rotational speed, M is the mass matrix, C is the damping matrix, K is the stiffness matrix, K Ω is the centrifugal stress stiffness, f is the excitation force amplitude, θ is the phase, t is the time parameter, NB is the number of blades. After obtaining the above equation, the direct integration method or the mode superposition method can be considered to solve the dynamic response in the modal coordinates; here, the Newmark-β method is used for solution, and the above equation is transformed into the following form: In the formula Both α and δ are the integration constants of the Newmark-β method. c0 to c 5 are dimensionless coefficients in the calculation process, K is the equivalent stiffness, different subscripts represent the equivalent stiffness at different times, Ω is the frequency, u is the displacement, M and C are the mass matrix and damping matrix of the system respectively, and F is the equivalent force amplitude at this moment.

[0006] Furthermore, in the first step, multiple coordinate systems are adopted in the calculation process. Among them, the flow field of the high-pressure turbine disk is calculated in the rotating coordinate system, and the flow fields of the other stator guide vane flow components are calculated in the stationary coordinate system; the total temperature, total pressure, and velocity direction are given at the inlet of the high-pressure turbine; the back pressure is given at the outlet; the adiabatic and no-slip boundary conditions are selected for the wall surface; the mixing plane method is selected for the interface when performing steady-state calculations; the moving / static blade slip boundary method is selected for the interface when performing unsteady calculations.

[0007] Furthermore, in the second step, the time-space Fourier decomposition can obtain the variation model of the unsteady aerodynamic force in the circumferential direction and time of the disk, but the variation of the aerodynamic force on the blade surface is ignored. Therefore, the TSMs model expresses the unsteady aerodynamic force in the following cylindrical coordinate form: where r is the radial coordinate from the hub to the tip, and θ is the circumferential coordinate of the disk.

[0008] Furthermore, in the second step, in order to appropriately simplify the equivalent excitation force, the blade excitation forces before and after the mode switching moment are calculated respectively to simulate the sudden change characteristics of the airflow under mode switching. Then, it is assumed that the excitation force on the blade changes with the rotational speed during the variable rotational speed process after mode switching. After calculating the excitation amplitudes at N rotational speeds respectively, the excitation force amplitude at any rotational speed during the variable rotational speed process can be calculated by cubic spline approximation:

[0009] Furthermore, in the fourth step, when the turbojet mode switches to the turbofan mode, the engine begins to enter the mode switching stage; in this stage, the combustion process of the main combustion chamber starts rapidly, generating a large amount of high-pressure airflow; the sudden load change in the mode switching stage will cause the dynamic response of the disk system to have shock fluctuations. According to the equivalent airflow model established in the second step, this excitation force model is characterized by the following formula. The amplitude of this equivalent model is different at different times, especially at the mode switching moment, the excitation amplitude increases rapidly: where C k , t k are the coefficients of cubic spline interpolation, θ(t) is a function of phase-time, EO is the excitation order, am is the amplitude of the exciting force at this order.

[0010] Furthermore, in the fourth step, due to the high-pressure nature of the air flow, the blade will experience a large acceleration process, the rotational speed will increase rapidly, and a relatively high rotational speed will be reached in a short time. The function θ(t) describing the phase change of the exciting force caused by the rotation of the disk can be calculated in the following form: where α is the rotational speed acceleration, t0 is the starting time, t is the time variable, ω0 is the starting rotational speed, and θ(t) is the phase-time function.

[0011] Furthermore, in the fourth step, the response solution formula is: where u is the displacement at this moment, F is the equivalent force at this moment, and K is the equivalent stiffness at this moment; thus, the dynamic response u at this time can be obtained.

[0012] Advantages of the present invention: This method uses the modal compression method to reduce the model dimension. The calculation model obtained by this method can greatly reduce the resources required for calculation, improve the calculation efficiency compared with the existing full-scale model calculation method. The comparison of the model scale and error is shown in Table 1 and Table 2, and the model scale directly affects the solution efficiency of the finite element model. In addition, in previous studies, no analysis method for the dynamic response of the disk under the mode switch of the dual variable cycle engine has been proposed. The method proposed in this patent fills the gap in the response prediction technology of the disk under sudden excitation. Brief Description of the Drawings

[0013] Figure 1 is the flow chart of the present invention;

[0014] Figure 2 Schematic diagram of establishing the CFD model of the variable cycle high-pressure turbine disk;

[0015] Figure 3 is the schematic diagram of the unsteady force in different regions of the blade;

[0016] Figure 4 is the schematic diagram of the force characterization model of the equivalent node during the mode switch;

[0017] Figure 5 is the schematic diagram of selecting the main node of the contact interface;

[0018] Figure 6 is the schematic diagram of the dynamic response of the disk considering different rotational speed accelerations during the mode switch;

[0019] Figure 7 is the schematic diagram of selecting the equivalent node of the turbine blade response;

[0020] Figure 8 Schematic diagram of the correction process curve of the equivalent excitation force at different excitation orders;

[0021] Figure 9 Schematic diagram of the response comparison before and after the equivalence of the excitation force at each excitation order. Specific implementation manner

[0022] Example 1

[0023] The present invention discloses a method for rapidly predicting the dynamic response under complex excitations of a dual variable cycle engine. The technical solution adopted includes the following steps: Step 1: Calculation of the steady force on the surface of the turbine disk during the mode switching process; establishing an air flow load model as the input for response calculation, and establishing an aerodynamic load model through CFD numerical simulation. Following the mass conservation equation, momentum conservation equation, and energy conservation equation, integrating the control volume to obtain a discrete equation set, and finding the dependent variables of the nodes. The discrete equation set is as follows: Among them, ρ is the density, p is the momentum, t is the time, U is the fluid velocity, τ is the stress tensor, λ is the thermal conductivity of the fluid, h tol is the total enthalpy, S M and S E are source terms; In the Cartesian coordinate system, the integral form of the control volume equation is as follows: Among them, V is the control volume, S is the control surface area, dn j is the outer normal vector perpendicular to the control surface; when performing unsteady calculations, a high-order solution mode is adopted for solution, the time term selects the second-order backward Euler format, and the turbulence model adopts the standard k-ω turbulence model, which includes the turbulent kinetic energy k (turbulent kinetic energy) equation and the turbulent frequency ω (turbulent frequency) equation: ρ is the density, t is the time, U is the velocity component, x is the spatial coordinate component, σ is the turbulent Prandtl number, β is the turbulent kinetic energy dissipation coefficient, μ is the turbulent dynamic viscosity, μ t is the turbulent viscosity, P is the potential energy; The established calculation model is as Figure 2As shown, multiple coordinate systems are adopted in the calculation process. Among them, the flow field of the high-pressure turbine disk is calculated in the rotating coordinate system, and the flow fields of the other stator guide vane flow components are calculated in the stationary coordinate system. The total temperature, total pressure, and velocity direction are given at the inlet of the high-pressure turbine; the back pressure is given at the outlet; the adiabatic and no-slip boundary conditions are selected for the wall surface; the mixing plane method is selected for the interface when performing steady-state calculations; the moving / stationary blade slip boundary method is selected for the interface when performing unsteady calculations; Step 2: Steady force TSM model on the disk surface; The distribution of external excitation in the cylindrical coordinate system is described by the Tyler-Sofrin excitation model (TSMs). The unsteady aerodynamic force is decomposed in time, and the harmonic components at a single engine order (EO) are extracted through the following formula: The complex Fourier coefficient P(h) extracted at the grid points on the blade surface is further spatially decomposed circumferentially to determine the rotational perturbation that will excite the corresponding traveling wave mode; The time-space Fourier coefficients are extracted at the corresponding points along the grid surface as follows: where m is the circumferential order, and P h m is the tangential distribution of the discrete-time Fourier coefficients at the corresponding surface points of the blade. N is the number of blades. In the current spatial Fourier decomposition, the circumferential DFT is performed on a sampling set that only takes N samples on each blade at the grid points at the corresponding positions on the surface; Therefore, according to the Nyquist theorem, when the number of blades is odd, the extracted rotational pressure components have multiple traveling waves from -N / 2 to +N / 2 or from -(N - 1) / 2 to +(N - 1) / 2; This method also takes into account the aliasing phenomenon experienced by the blade row when the excitation order is higher than the traveling wave excitation of N / 2; The time-space Fourier decomposition can obtain the variation model of the unsteady aerodynamic force in the circumferential direction and time of the disk, but ignores the variation of the aerodynamic force on the blade surface; Therefore, the TSMs model expresses the unsteady aerodynamic force in the following cylindrical coordinate form: Among them, r is the radial coordinate from the hub to the tip, and θ is the circumferential coordinate of the disk. Here, the blade is regarded as a simple cantilever beam, considering the variation form of the aerodynamic excitation from the hub to the tip of the disk with the radial coordinate r. However, the geometric characteristics of different blade heights and chord lengths are closely related to the local flow field change factors. Analyzing the change of the excitation force solely from the radial coordinate will ignore the significant differences in the aerodynamic load distribution in different regions of the blade. Therefore, by subdividing the blade into multiple segments along the chord length direction, the distribution change law of the aerodynamic load in each region of the blade can be simulated more precisely, the aerodynamic excitation of the disk under actual working conditions can be reflected more comprehensively, and a more accurate load model can be provided for the vibration analysis and optimization of the turbine disk of the dual variable cycle engine. By equivalenting the unsteady aerodynamic force data of all nodes to the excitation of a few key nodes, according to the equivalent principle, the unsteady aerodynamic forces of each node obtained by using computational fluid dynamics are equivalent based on the dynamic response, and reducing them to a very small number of excitation nodes can greatly reduce the dimension of the calculation model. During the equivalent process, first calculate the airflow excitation response at this rotational speed, then calculate the initial response on the simplified nodes, and judge whether the response is equivalent according to the root mean square error between the simplified node response and the airflow excitation response. If not equivalent, update and iterate the numerical value of the node force according to the sensitivity coefficient of the node force to the response function until the response function meets the requirements. During the equivalent process, calculate the sensitivity matrix by the finite difference method, perturb the excitation force and calculate the response to obtain the sensitivity: Among them, R(F) is the response solution function, F is the excitation force. After forming the sensitivity matrix S through the sensitivity, the update formula of the node force amplitude is based on the least square method, and the pseudo-inverse of the sensitivity matrix S is used to minimize the response error: Among them, k is the error iteration number. For the dual variable cycle model, three nodes with different chord lengths are selected to equivalent the excitation force; the node selection is as Figure 7 shown, and the error iteration descent diagram under each order of excitation is as Figure 8 shown. Finally, the equivalent results of the excitation force under each excitation order are as Figure 9 shown; In order to appropriately simplify the equivalent excitation force, calculate the blade excitation forces before and after the mode switching moment respectively to simulate the sudden change characteristics of the airflow under mode switching. Then assume that the excitation force on the blade changes with the rotational speed during the variable rotational speed process after mode switching. After calculating the excitation amplitudes at N rotational speeds respectively, the excitation force amplitude at any rotational speed during the variable rotational speed process can be calculated through cubic spline approximation: Finally, the interpolated excitation is as Figure 3 shown; Step 3: The main node reduced-order model of the turbine disk calculates the dynamic response of the disk system through the model reduction method of the main nodes. During the calculation, the equivalent excitation nodes, response nodes, and contact interface nodes on the surface are selected as the main nodes. The selection of the main nodes is as shown in Figure 5 and the model scales before and after reduction are as follows: Table 1 Model scales of the high-pressure turbine disk finite element model before and after reduction Table 2 Comparison of the natural frequencies of the model before and after reduction The credibility of the model obtained after reduction is analyzed by using the frequency difference to analyze the correlation between the original model and the reduced model. Among them, the frequency differences of the natural frequencies of the node diameter modes concerned before and after reduction are shown in Table 2. It can be seen that the maximum frequency difference is only 3.89%, and the average frequency difference is less than 1%, which proves the accuracy of the constructed main node model; The modal reduction method based on the main nodes is more suitable for application scenarios with more contact nodes and when performing nonlinear calculations; for the multi-harmonic motion equation, only the "nonlinear" degrees of freedom are retained, that is, the degrees of freedom applying the nonlinear force, while reducing all other "linear" degrees of freedom; the calculation can be completed by simplifying the model. During this process, the conversion from the physical domain to the modal domain is realized through the FRF matrix, and the reduction of the model can be achieved. The FRF matrix of the multi-harmonic motion equation is as follows: Q j = A(m j ω)(P j - F j (Q)) Q j is the overall harmonic expansion matrix, A is the overall frequency response function matrix, P j is the dynamic stiffness matrix, F j is the amplitude of the excitation force in the frequency domain, Q is the harmonic component, and then the nonlinear harmonic component Q n and the linear harmonic component Q l are divided. For each harmonic, it can be written in the following form: A ll is the linear frequency response function matrix, A nl and A ln are the frequency response function matrices of the linear and nonlinear couplings, A nn is the frequency response function FRF matrix of the nonlinear term, Q n is the harmonic component of the nonlinear term. The equation can be simplified into the following form: Among them, It is the modal displacement vector of the nonlinear degrees of freedom calculated without nonlinear contact interaction forces. Utilizing the modal characteristics of the linear structure, the FRF matrix is usually approximated by its modal expansion where, A nn is the processed FRF matrix of the nonlinear term frequency response function, Φ j is the modal shape, ω is the frequency of the calculated order, N is the number of retained modal orders, ω j is the frequency at this order, and η is the modal damping ratio;

[0024] Step 4: The reduced-order model of the main nodes of the turbine disk. When the dual variable cycle engine switches from the turbojet mode to the turbofan mode, the disk system roughly experiences three stages in chronological order, namely the turbojet windmill stage, the instantaneous mode switching stage, and the turbofan mode stable stage. When the turbojet mode switches to the turbofan mode, the engine begins to enter the mode switching stage; in this stage, the combustion process in the main combustion chamber starts rapidly, generating a large amount of high-pressure air flow; this switch will trigger a significant transient impact, and the disk load instantaneously switches from the original windmill air flow to the high-pressure air flow generated by the main combustion chamber; this sudden change in load will cause the dynamic response of the disk system to have shock fluctuations. According to the equivalent air flow model established in Step 2 as Figure 4 shown, the excitation in the figure increases rapidly during the mode switching process. Characterize this excitation force model as the following formula. The amplitude of this equivalent model is different at different times, especially at the moment of mode switching, the excitation amplitude increases rapidly: where, C k , t k is the coefficient of cubic spline interpolation, θ(t) is a function of phase-time, EO is the excitation order, a m is the excitation force amplitude at this order; Due to the high-pressure nature of the air flow, the blades will experience a large acceleration process, the rotational speed will increase rapidly, and reach a relatively high rotational speed in a short time. The function θ(t) describing the phase change of the excitation force caused by the rotation of the disk can be calculated in the following form: where, α is the rotational speed acceleration, t0 is the starting time, t is the time variable, ω0 is the starting rotational speed, and θ(t) is a function of phase-time; During the process of mode switching, the change in the amplitude of the airflow excitation is closely related to the change in rotational speed; under the excitation of the high-pressure airflow in the main combustion chamber, the amplitude of the airflow will gradually increase and be proportional to the rotational speed as the rotational speed of the disk increases; with the continuous action of the high-pressure airflow, the rotational speed of the disk keeps rising, and the excitation amplitude of the airflow will also increase accordingly; at this time, the air excitation amplitude increases with the increase in rotational speed; until the high-pressure turbine disk reaches the turbofan mode rotational speed and the rotational speed change stops, the airflow excitation amplitude tends to be stable. The dynamic equation of the blade is as follows: During the rotation of the disk, the hardening effect of the centrifugal force will be exerted, making the overall stiffness of the disk show an increasing trend with the rotational speed. At this time, the dynamic equation can be expressed as: where, Ω is the rotational speed of the disk, M is the mass matrix, C is the damping matrix, K is the stiffness matrix, K Ω is the centrifugal stress stiffness, f is the excitation force amplitude, θ is the phase, t is the time parameter, NB is the number of blades. After obtaining the above equation, the direct integration method or the mode superposition method can be considered to solve the dynamic response in the modal coordinates; here, the Newmark-β method is used for solution, and the above formula is transformed into the following form: In the formula where both α and δ are the integration constants of the Newmark-β method. To ensure the unconditional stability of the algorithm, c0 to c 5 are dimensionless coefficients in the calculation process, K is the equivalent stiffness, different subscripts represent the equivalent stiffness at different times, Ω is the frequency, u is the displacement, M and C are the mass matrix and damping matrix of the system, and F is the equivalent force amplitude at this moment. The response solution formula is: where, u is the displacement at this moment, F is the equivalent force at this moment, and K is the equivalent stiffness at this moment. α and δ can be taken as 0.25 and 0.5 respectively, and then the dynamic response u at this time can be obtained. The response solution results under different rotational speed accelerations are as Figure 6 shown. In summary, this patent first establishes a spatio-temporal distribution model of the airflow load, considering the time-varying and mutation effects of the airflow load during the mode switching process of the variable cycle engine. By means of response node equivalence, the airflow load is applied to the key positions of the bladed disk to form a high-precision airflow load excitation model, ensuring that the model can accurately describe the excitation effect of the airflow load on the bladed disk. The modal compression method is used to simplify the bladed disk model by compressing high-order modes and retaining key nodes as master nodes. This step can minimize the dimension of the model and retain the key features of the dynamic response of the bladed disk. The time-domain integration method is used to calculate the simplified bladed disk dynamic model. According to the constructed airflow load excitation model, the transient non-linear dynamic response of the bladed disk under mode switching is calculated by the time-domain integration method. During this process, considering the mutation effect of the excitation load caused by mode switching, the transient response of the bladed disk under different working modes is obtained, and its dynamic characteristics are analyzed.

[0025] Components not described in detail in this article are prior art.

[0026] Although the specific embodiments of the present invention have been described in detail above, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those of ordinary skill in the art, various changes can be made without departing from the gist of the present invention, and modifications or deformations that do not involve creative labor are still within the protection scope of the present invention.

Claims

1. A fast prediction method for the dynamic response of a dual-cycle engine under complex excitation, characterized in that: The following steps are involved: Step 1: Calculate the steady force on the turbine blade surface during mode switching; establish the airflow load model as the input of the response calculation, and establish the aerodynamic load model through CFD numerical simulation. Follow the mass conservation equation, momentum conservation equation and energy conservation equation, integrate the control volume, obtain the discrete equation group, and find the dependent variables of the node. The discrete equation group is as follows: Where ρ is density, p is momentum, t is time, U is the velocity of the fluid, τ is the stress tensor, λ is the thermal conductivity of the fluid, and h tol is the total enthalpy, S M and S E is the source term; In the Cartesian coordinate system, the integral form of the control volume equation is as follows: Where V is the volume of the control body, S is the surface area of ​​the control body, and dn j is the external normal vector perpendicular to the surface of the control body; when performing unsteady calculations, a high-order solution mode is used for solving, the time term selects the second-order backward Euler format, and the turbulence model adopts the standard k-ω turbulence model, which includes the turbulent kinetic energy k (turbulent kinetic energy) equation and the turbulent frequency ω (turbulent frequency) equation: k equation: Omega equation: ρ is density, t is time, U is velocity component, x is spatial coordinate component, σ is turbulent Prandtl number, β is turbulent kinetic energy dissipation coefficient, μ is turbulent dynamic viscosity, μ t is the turbulent viscosity, P is the potential energy; Step 2: TSM model of steady force on the blade surface; the distribution of external excitation in the cylindrical coordinate system is described by the Tyler-Sofrin excitation model (TSMs), the unsteady aerodynamic force is decomposed in time, and the harmonic components at a single engine order (EO) are extracted by the following formula: The complex Fourier coefficients P(h) extracted at the blade surface grid points are further spatially decomposed along the circumferential direction to determine the rotational disturbance that will excite the corresponding traveling wave mode; the time-space Fourier coefficients are extracted along the corresponding points of the grid surface as follows: Where m is the circumferential order, P h m is the tangential distribution of the discrete-time Fourier coefficients at the corresponding surface point of the blade, N is the number of blades, and in the current spatial Fourier decomposition, the circumferential DFT is performed on a sampling set of only N samples on each blade at the grid point at the corresponding position on the surface; By equating the unsteady aerodynamic data of all nodes to the excitation of a few key nodes, these key nodes are equivalent to the unsteady aerodynamics of each node obtained by computational fluid dynamics according to the equivalent principle. Reducing them to a very small number of excitation nodes can greatly reduce the dimension of the computational model. In the equivalent process, the airflow excitation response at the speed is first calculated, and then the initial response on the simplified node is calculated. The root mean square error between the simplified node response and the airflow excitation response is used to determine whether the response is equivalent; if not equivalent, the value of the node force is updated and iterated according to the sensitivity coefficient of the node force to the response function until the response function meets the requirements. In the equivalent process, the sensitivity matrix is ​​calculated by the finite difference method, the excitation force is disturbed and the response is calculated to obtain the sensitivity: Among them, R(F) is the response solution function, F is the exciting force, and after the sensitivity matrix S is constructed through sensitivity, the update formula of the node force amplitude is based on the least squares method, and the pseudo-inverse of the sensitivity matrix S is used to minimize the response error: Where k is the number of error iterations. For the double variable cycle model, three nodes with different chord lengths are selected to equivalently excite the force; Step 3: The main node dimensionality reduction model of the turbine blade disk. The dynamic response of the blade disk system is calculated by the main node model dimensionality reduction method. In the calculation process, the equivalent excitation nodes, response nodes and contact interface nodes of the surface are selected as the main nodes. The modal dimensionality reduction method based on the main node is more suitable for application scenarios with many contact nodes and nonlinear calculations. For the multi-harmonic motion equation, only the "nonlinear" degree of freedom is retained, that is, the degree of freedom for applying nonlinear forces, and all other "linear" degrees of freedom are reduced. The calculation can be completed by simplifying the model. In this process, the conversion from the physical domain to the modal domain is realized through the FRF matrix, which can achieve the dimensionality reduction of the model. The FRF matrix of the multi-harmonic motion equation is as follows: Q j =A(m j ω)(P j -F j (Q)) Q j is the overall harmonic expansion matrix, A is the overall frequency response function matrix, P j is the dynamic stiffness matrix, F j is the frequency domain exciting force amplitude, Q is the harmonic component, and then the nonlinear harmonic component Q is divided n and linear harmonic components Q l , for each harmonic, it can be written as follows: A ll is the linear frequency response function matrix, A nl With A ln is the frequency response function matrix of linear and nonlinear coupling, A nn Nonlinear frequency response function FRF matrix, Q n is the harmonic component of the nonlinear term, the equation can be simplified to the following form: in, is the modal displacement vector of the nonlinear degrees of freedom calculated in the absence of nonlinear contact interaction forces; using the modal properties of the linear structure, the FRF matrix is ​​usually approximated by its modal expansion Among them, A nn is the processed nonlinear frequency response function FRF matrix, Φ j is the modal vibration shape, ω is the calculated order frequency, N is the number of retained modal orders, ω j is the frequency at this order, η is the modal damping ratio; Step 4: Main node dimensionality reduction model of turbine blade disk. When the dual-cycle engine switches from turbojet to turbofan mode, the blade disk system roughly goes through three stages in chronological order, namely, the turbojet windmill stage, the mode switching instantaneous stage, and the turbofan mode stable stage. During the mode switching process, the change in the airflow excitation amplitude is closely related to the speed change; under the excitation of the high-pressure airflow in the main combustion chamber, the airflow amplitude will gradually increase, and will be proportional to the speed as the blade disk speed increases; with the continuous action of the high-pressure airflow, the blade disk speed continues to rise, and the airflow excitation amplitude will also increase accordingly; at this time, the air excitation amplitude increases with the increase in speed; until the high-pressure turbine blade disk reaches the turbofan mode speed, the speed change stops, the airflow excitation amplitude tends to be stable, and the blade dynamics equation is as follows: The blade disk will be hardened by the centrifugal force during the rotation process, so that the overall stiffness of the blade disk tends to increase with the rotation speed. At this time, the dynamic equation can be expressed as: Among them, Ω is the blade speed, M is the mass matrix, C is the damping matrix, K is the stiffness matrix, K Ω is the centrifugal stress stiffness, f is the exciting force amplitude, θ is the phase, t is the time parameter, NB is the number of blades, and the above equation is obtained; the Newmark-β method is used here to solve, and the above equation is converted into the following form: In the formula Where α and δ are the integral constants of the Newmark-β method, c0 to c5 is the dimensionless coefficient in the calculation process, K is the equivalent stiffness, different subscripts represent the equivalent stiffness at different times, Ω is the frequency, u is the displacement, M, C are the mass matrix and damping matrix of the system, and F is the equivalent amplitude at this moment.

2. The method for rapid prediction of dynamic response of a dual-cycle engine under complex excitation according to claim 1 is characterized in that: In the step 1, multiple coordinate systems are used in the calculation process, wherein the flow field of the high-pressure turbine blade disk is calculated in a rotating coordinate system, and the flow field of the remaining stator guide vane flow components is calculated in a stationary coordinate system; the total temperature, total pressure and velocity direction are given at the inlet of the high-pressure turbine; the back pressure is given at the outlet; the wall surface is selected to be adiabatic and has no slip boundary conditions; when performing steady calculations, the interface selects the mixed plane method; when performing unsteady calculations, the interface selects the moving / stationary blade sliding boundary method.

3. The method for rapid prediction of dynamic response of a dual-cycle engine under complex excitation according to claim 1 is characterized in that: In step 2, the time-space Fourier decomposition can obtain the variation model of the unsteady aerodynamic force in the circumferential direction of the blade and in time, but ignores the variation of the aerodynamic force on the blade surface; therefore, the TSMs model expresses the unsteady aerodynamic force in the following cylindrical coordinate form: Among them, r is the radial coordinate from the hub to the blade tip, and θ is the circumferential coordinate of the blade disk.

4. The method for rapid prediction of dynamic response of a dual-cycle engine under complex excitation according to claim 1 is characterized in that: In order to appropriately simplify the equivalent exciting force in step 2, the blade exciting force before and after the mode switching moment is calculated to simulate the sudden change characteristics of the airflow under mode switching. Then, it is assumed that the exciting force on the blade changes with the speed during the variable speed process after the mode switching. After calculating the excitation amplitudes at N speeds, the exciting force amplitude at any speed during the variable speed process can be calculated by cubic spline approximation:

5. According to the method for rapid prediction of dynamic response of a dual-cycle engine under complex excitation according to claim 1, in the step 4, when the turbojet mode is changed to the turbofan mode, the engine begins to enter the mode switching stage; the sudden load change in the mode switching stage will cause the dynamic response of the blade disk system to have impact fluctuations, and according to the equivalent airflow model established in step 2, the exciting force model is characterized by the following formula, and the amplitude of this equivalent model is different at different times, especially at the mode switching moment, the excitation amplitude increases rapidly: in, C k , t k is the coefficient of cubic spline interpolation, θ(t) is the phase-time function, EO is the excitation order, a m is the exciting force amplitude at this order.

6. According to the method for rapid prediction of dynamic response of a dual-cycle engine under complex excitation according to claim 1, in the step 4, due to the high-pressure nature of the airflow, the blades will experience a large acceleration process, the speed will increase rapidly, and will reach a higher speed in a short time. The function θ(t) describing the phase change of the exciting force caused by the rotation of the blade disk can be calculated in the following form: in, α is the speed acceleration, t0 is the starting time, t is the time variable, ω0 is the starting speed, and θ(t) is the phase-time function.

7. The method for rapid prediction of dynamic response of a dual-cycle engine under complex excitation according to claim 1 is characterized in that: In step 4, the response solution formula is: Among them, u is the displacement at this moment, F is the equivalent force at this moment, and K is the equivalent stiffness at this moment; the dynamic response u at this moment can be obtained.

Images