Method and system for determining folding steering response of high-order small-amplitude-error aircraft

Abstract

The invention provides a high-order small-amplitude-error aircraft folding steering response determination method and system, and belongs to the technical field of high-dimensional nonlinear dynamic response analysis. The problems that according to an existing single-step single-solution implicit dynamic response determination method, consistent second-order precision, controllable numerical value high-frequency dissipation and zero-order overshoot cannot be achieved at the same time, the dynamic response of the aircraft folded rudder cannot be effectively solved, and then follow-up flutter and structural optimization analysis of the folded rudder are hindered are solved. According to the method, a high dissipation mechanism is introduced, numerical instability caused by nonlinearity, high dimensionality and false high-frequency components is effectively reduced, and meanwhile, the calculation precision can be remarkably improved on the premise that unconditional stability is guaranteed in combination with the low dispersion characteristic; the method is particularly suitable for determining the dynamic response of the folding rudder of the aircraft in aerospace, and has good adaptability and advantages.

Description

Technical Field

[0001] The present invention relates to the technical field of high-dimensional non-linear dynamic response analysis in aerospace, and in particular, to a method and system for determining the dynamic response of a folding rudder of an aircraft with high-order small-amplitude errors. Background Art

[0002] In the high-dimensional non-linear dynamic response analysis in aerospace, how to accurately and quickly calculate the dynamic response of the system has always been a difficult problem in the engineering field. Traditional numerical methods, such as explicit time-domain integrators, although having certain advantages in the single-step calculation speed, when dealing with the non-linearity (rigidity) of the folding rudder of an aircraft, due to their low numerical stability and strong dispersion effect, it is impossible to efficiently obtain accurate dynamic responses. Implicit time-domain integrators, with their good stability and high accuracy, can better handle high-dimensional non-linear problems. Although the implicit time-domain integrator has a higher computational cost within a single step, it can reduce the total computational cost of solving the folding rudder of an aircraft by selecting a larger integration step size with unconditional stability.

[0003] Traditional single-step single-solution implicit integration algorithms (such as the Newmark-β method [1] , the Crank-Nicolson method [2] etc.) are widely used in dealing with dynamic response problems, but there are also certain limitations and deficiencies. Typically, among all the algorithm members of the Newmark family, only the trapezoidal formula is optimal with respect to the accuracy order and overshoot, but it does not converge when solving strong non-linear dynamic problems [3] . The Newmark method cannot introduce numerical high-frequency dissipation when achieving second-order accuracy, and there is a first-order overshoot of displacement and velocity when introducing numerical high-frequency dissipation by reducing the accuracy. To solve the drawback that the second-order Newmark method cannot introduce numerical high-frequency dissipation, scholars have made a series of improvements to the Newmark-β method from different algorithm design perspectives. Domestic and foreign scholars have successively proposed three classic α-family methods, namely HHT-α [4] , WBZ-α [5] and TPO / G-α [6,7] . All three α methods have a first-order overshoot of displacement and velocity and only achieve first-order accuracy for acceleration. In addition, the HHT-α method cannot achieve the full-process change of numerical high-frequency dissipation, that is, it cannot achieve the continuous change from zero dissipation to asymptotic elimination characteristics. To improve the first-order overshoot of the three α methods, Yu Kaiping [8]A class of single-solution implicit integration methods with seven parameters is designed, and a zero-order overshoot algorithm corresponding to three types of α methods is proposed through accuracy, stability, and overshoot analysis. However, due to incomplete accuracy analysis (lacking accuracy analysis of velocity and acceleration), the developed zero-order overshoot algorithm has a degradation in accuracy order, that is, it cannot achieve second-order accuracy in solving damped forced vibrations. In addition, some scholars have proposed some single-solution implicit dissipative algorithms from different algorithm design perspectives, such as the Wilson-θ method [9] , ρ-method

[10] , Θ format

[11] , SSpj method

[12] , single-step Houbolt algorithm

[13] (SSH) and optimal collocation format

[14] etc. None of these algorithms have achieved the six-point algorithm properties mentioned above. For example, the Wilson-θ method and the optimal collocation format achieve second-order accuracy in displacement, velocity, and acceleration, but there are second-order displacement and first-order velocity overshoots; at the same time, these two types of algorithms do not achieve the full-course variation of numerical high-frequency dissipation, and the high-frequency dissipation is not easy to adjust. The ρ-method achieves the full-course variation of numerical high-frequency dissipation and the high-frequency dissipation amount can be precisely adjusted through its parameter ρ, but it only achieves first-order accuracy and there is first-order velocity overshoot. Although the author who proposed the Θ format claims that this algorithm has achieved the six-point algorithm properties mentioned above, this algorithm actually only achieves first-order accuracy and there are first-order displacement and first-order velocity overshoots. At the same time, the numerical high-frequency dissipation of the Θ format is not adjustable. The SSH method also has the disadvantages of first-order overshoot and first-order accuracy.

[0004] Literature

[15] has proved that traditional single-step implicit integrators cannot simultaneously achieve consistent second-order accuracy, controllable numerical high-frequency dissipation, and zero-order overshoot, thus also restricting the further development of traditional single-step implicit integrators and their effective application in aircraft folding rudders. These limitations and problems have promoted the research on new methods for determining dynamic responses. Especially in the analysis of high-dimensional nonlinear dynamic responses such as aircraft folding rudders, how to improve the computational efficiency and accuracy of algorithms by introducing mechanisms with high dissipation and low dispersion has become one of the current research focuses. Summary of the Invention

[0005] The technical problem to be solved by the present invention is:

[0006] To solve the problem that the existing single-step single-solution implicit dynamic response determination method cannot simultaneously achieve consistent second-order accuracy, controllable numerical high-frequency dissipation, and zero-order overshoot, and thus cannot effectively solve the dynamic response of aircraft folding rudders, which hinders the subsequent flutter and structural optimization analysis of folding rudders.

[0007] The technical solution adopted by the present invention to solve the above technical problem:

[0008] The present invention provides a method for determining the dynamic response of a folding rudder of an aircraft with high-order small-amplitude errors, comprising the following steps:

[0009] S100. The partial differential equation governing the dynamic response of the folding rudder of the aircraft is degenerated into a system of ordinary differential equations in the following form after spatial discretization:

[0010]

[0011] In the formula, M represents the global matrix; f represents the contributions of external and internal forces; u(t) represents the nodal state vector of the folding rudder of the aircraft; represents the time derivative of the nodal state vector; u0 represents the value of the nodal state vector at the initial time t0;

[0012] S200. Assume that the integration duration for analyzing the folding rudder of the aircraft is t ∈ [0, T]. First, divide this integration interval into discrete time intervals, i.e., where N represents the total number of integration steps. Then, construct an integration scheme for the small integration interval t ∈ [t n , t n+1 to update the dynamic response at time t n to time t n+1 . Finally, recursively iterate this process to solve the dynamic response of the folding rudder of the aircraft over the entire integration interval;

[0013] S300. Introduce constraint conditions to determine the design parameters to be determined during the process of solving the dynamic response of the folding rudder of the aircraft over the entire integration interval in step S200, so as to make the system have optimal spectral characteristics, consistent second-order accuracy, controllable numerical high-frequency dissipation ability, and unconditional stability.

[0014] Further, in step S200, it specifically includes:

[0015] S210. Divide the integration interval t ∈ [t n , t n+1 into three small intervals:

[0016] t ∈ [t n , t n + γDt] ∪ [t n + γDt, t n + ζDt] ∪ [t n + ζDt, t n+1 (2)

[0017] In the formula, γ and ζ represent the two sub-step size division rates of the folding rudder of the aircraft within the current integration interval, and satisfy γ ≤ ζ; Dt represents the integration step size for solving the folding rudder of the aircraft;

[0018] S220. Utilize t nKnown node state vector u at a moment n Solve for t n+γ := t n + γ Node state quantity at Dt moment:

[0019]

[0020] S230. After completing step S220, use the known node state vector u n and node derivative vector Solve for t n+ζ := t n + ζ Node state quantity at Dt moment:

[0021]

[0022] S240. After completing steps S220 and S230, use the known node state vector u n and node derivative vector Solve for t n+1 Node state quantity at moment:

[0023]

[0024] S250. After completing the above calculations, use the known node state vector u n and node derivative vector Calculate t n+1 Node state quantity at moment:

[0025]

[0026] Perform iterative calculations step by step to obtain the structural dynamic response of the folding rudder of the aircraft within the entire integration interval.

[0027] Furthermore, in step S300, it specifically includes: determining the design parameters γ, ζ, α, μ, σ, β1, and β2 in step S200 by introducing constraint conditions

[0028] When the design parameters satisfy formula (7), the dynamic response determination method of the folding rudder of the aircraft realizes a consistent effective stiffness matrix, which reduces the computational cost of solving the folding rudder of the aircraft in the linear setting and realizes the optimal spectral characteristics at this time;

[0029] γ = α = σ (7)

[0030] The design parameters are solved to achieve the consistent second-order accuracy of the folding rudder of the aircraft by satisfying the following conditions:

[0031]

[0032] To effectively reduce numerical instability, a controllable numerical high-frequency dissipation ability is required, that is:

[0033]

[0034] In the formula, ρ ∞ ∈[0, 1] represents a user-specified parameter; the strongest and weakest numerical high-frequency dissipations are introduced when ρ ∞ = 0 and ρ ∞ = 1 respectively. When ρ ∞ = 1, it is called zero dissipation. When solving the folding rudder of the aircraft subsequently, it is usually required that the user selects a ρ ∞ value less than 1;

[0035] The parameter μ satisfies the following conditions so that the numerical solution predicted by formula (6) and the solution predicted by formula (5) are second-order accurate:

[0036]

[0037] To determine the redundant integration parameters, it is required that the sizes of the first two sub-steps of the dynamic response determination method of the folding rudder of the aircraft are equal, that is:

[0038] ζ = 2γ. (11)

[0039] Solving formulas (7) to (11) gives:

[0040]

[0041] When the parameters γ and ρ ∞ satisfy the inequality constraint (13), the dynamic response determination method has unconditional stability:

[0042] 24(ρ ∞ + 1)γ 4 - 8(ρ ∞ + 9)γ 3 + 48γ 2 - 12γ + 1 ≤ 0 (13)

[0043] After the user gives the ρ ∞ value, the value range of the parameter γ, γ ∈ [γ min , γ max , is obtained by solving formula (13) to achieve unconditional stability.

[0044] Furthermore, first set ρ ∞ ∈[0, 1] to determine the numerical high-frequency dissipation amount, then solve the inequality (13) to obtain the value range of the parameter γ, γ ∈ [γ min , γ max , and select the γ value in this range through the following steps;

[0045] The main terms of the amplitude (δ) and phase (ε) errors calculated by the amplitude and phase error calculation method are as follows:

[0046]

[0047] In the formula, c4 is the value of the irrelevant parameter; ω is the natural frequency of the linear single-degree-of-freedom system; ξ represents the damping ratio of the linear single-degree-of-freedom system; O represents the high-order quantity;

[0048] Formula (14) shows that when analyzing the undamped problem and the parameter values are γ = γ min or γ = γ max at this time, the amplitude error δ is increased in order, reaching O(Dt 5 ), and thus the reduction of the amplitude error is:

[0049]

[0050] The method for determining the dynamic response of the folding rudder of an aircraft obtains a higher-order amplitude error order and the minimum phase error when the parameter values are as in formula (12) and the parameter γ value is γ min at this time.

[0051] Furthermore, it is used to solve the nine-point acceleration response value of the folding rudder system of the aircraft for clearance nonlinearity.

[0052] A system for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error, the system has a program module corresponding to the above steps, and executes the steps in the method for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error when running.

[0053] A computer-readable storage medium, the computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the method for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error when called by a processor.

[0054] Compared with the prior art, the beneficial effects of the present invention are:

[0055] The method and system for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error of the present invention effectively reduce the numerical instability caused by nonlinearity, high dimension and false high-frequency components by introducing a high-dissipation mechanism, and at the same time combine the low-dispersion characteristics, so that on the premise of ensuring unconditional stability, the calculation accuracy can be significantly improved. This method is particularly suitable for determining the dynamic response of the folding rudder of an aircraft in aerospace, and has good adaptability and advantages. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 It is a schematic diagram of the sub-step division of the method of the present invention in the embodiment of the present invention;

[0057] Figure 2 is the curve graph of the parameter γ varying with ρ in the embodiment of the present invention ∞ ;

[0058] Figure 3 is the curve graph of the main item of the phase error of the method of the present invention in the embodiment of the present invention

[0059] Figure 4 is the comparison graph of the spectral characteristics of the method of the present invention in the embodiment of the present invention. Among them, (a) is the comparison graph of the spectral radius change under different numerical high-frequency dissipation amounts, (b) is the comparison graph of the numerical damping ratio change under different numerical high-frequency dissipation amounts, and (c) is the comparison graph of the relative period error change under different numerical high-frequency dissipation amounts

[0060] Figure 5 is the comparison graph of the convergence rate change of the method of the present invention in the embodiment of the present invention. Among them, (a) is the comparison graph of the global displacement error change under different numerical high-frequency dissipation amounts, (b) is the comparison graph of the global velocity error change under different numerical high-frequency dissipation amounts, and (c) is the comparison graph of the global acceleration error change under different numerical high-frequency dissipation amounts

[0061] Figure 6 is the overshoot response change graph of the method of the present invention in the embodiment of the present invention. Among them, (a) is the comparison graph of the displacement response change under different numerical high-frequency dissipation amounts, and (b) is the comparison graph of the velocity response change under different numerical high-frequency dissipation amounts

[0062] Figure 7 is the comparison graph of the amplitude and phase error of the SS-VA method, SUCI2 method, the method of the present invention, and the OESS method under the same sub-step size. Among them, (a) is the comparison graph of the amplitude error of each implicit integration method under the same sub-step size, and (b) is the comparison graph of the phase error of each implicit integration method under the same sub-step size

[0063] Figure 8 is the flow chart of the method of the present invention for solving the first-order dynamic system in the embodiment of the present invention

[0064] Figure 9 is the flow chart of the method of the present invention for solving the second-order dynamic system in the embodiment of the present invention

[0065] Figure 10 is the schematic diagram of the folding rudder system of the aircraft with clearance nonlinearity in the embodiment of the present invention. Among them, (a) is the schematic diagram of the model, and (b) is the schematic diagram of the finite element mesh of the model

[0066] Figure 11Schematic diagram of acceleration response of nine test points when the method of the present invention is used to solve the nonlinear clearance of an aircraft folding rudder in an embodiment of the present invention, wherein (a)-(i) are schematic diagrams of acceleration response of the nine test points respectively. DETAILED DESCRIPTION

[0067] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0068] The present invention aims to solve the problem that the existing single-step single-solution implicit step-by-step integration method cannot simultaneously achieve consistent second-order accuracy, controllable numerical high-frequency dissipation and zero-order overshoot, and therefore cannot effectively solve the dynamic response of the aircraft folding rudder, thereby hindering the subsequent flutter and structural optimization analysis of the folding rudder. Based on this, the present invention proposes a high-order small amplitude error dynamic response determination method for predicting the dynamic response results of the aircraft folding rudder. The present invention can provide an accurate numerical solution for the prediction of high-dimensional nonlinear dynamic systems in aerospace, thereby accelerating the transformation of structures from theoretical design to practical application.

[0069] Specific implementation plan 1: Combine Figures 1 to 3 As shown, the present invention provides a method for determining the folding rudder dynamic response of an aircraft with a high-order small amplitude error, comprising the following steps:

[0070] S100. The partial differential equation governing the dynamic response of the folding rudder structure of the aircraft is degenerated into a set of ordinary differential equations in the following form after spatial discretization:

[0071]

[0072] Where M is the global matrix; f represents the contribution of external and internal forces; u(t) represents the node state vector of the folding rudder of the aircraft; represents the time derivative of the node state vector; u0 represents the node state vector value at the initial time t0; assuming that the integral time length of the folding rudder of the analyzed aircraft is t∈[0,T], the determination method of the present invention first divides the integral interval into discrete time intervals, namely Where N represents the total number of integration steps, and then for the small integration interval t∈[t n ,t n+1 ] Construct a suitable integral format to convert t n The dynamic response at time t is updated to n+1 moment, and finally, by recursively performing this process step by step, the dynamic response of the aircraft's folding rudder in the entire integral interval can be solved;

[0073] S200, the method for determining the aircraft folding rudder dynamic response of high-order small amplitude error is implemented according to the following steps:

[0074] S210, fold the aircraft rudder in the current integration interval t∈[tn , t n+1 is refined into three small intervals as shown in Figure 1 :

[0075] t ∈ [t n , t n + γDt] ∪ [t n + γDt, t n + ζDt] ∪ [t n + ζDt, t n+1 (2)

[0076] In the formula, γ and ζ are the two sub-step size division rates of the folding rudder of the aircraft in the current integration interval and satisfy γ ≤ ζ; Dt represents the integration step size for solving the folding rudder of the aircraft;

[0077] S220. Use the known node state vector u n at time t n to solve the node state quantity at t n+γ : = t n + γDt:

[0078]

[0079] S230. After completing step S220, use the known node state vector u n and the node derivative vector to solve the node state quantity at t n+ζ : = t n + ζDt:

[0080]

[0081] S240. Similarly, after completing steps S220 and S230, use the known node state vector u n and the node derivative vector to solve the node state quantity at t n+1 :

[0082]

[0083] S250. After completing the above calculations, use the known node state vector u n and the node derivative vector to calculate the node state quantity at t n+1 :

[0084]

[0085] Steps S210 to S220 describe how the present invention calculates the dynamic response (u n at time t n and ) Update to the dynamic response (u n+1 ) at time t n+1 and ), and it belongs to the implicit time domain integration method; by such step - by - step iterative calculation, the present invention can obtain the structural dynamic response of the folding rudder of the aircraft within the entire integration interval; note that the dynamic response determination method of the present invention does not require solving the derivative vector at the initial time to start the integration process, so it is directly self - starting;

[0086] S300. The dynamic response determination method of the present invention introduces undetermined design parameters: γ, ζ, α, μ, σ, β1, and β2. Therefore, constraint conditions need to be introduced to determine the unknown undetermined parameters;

[0087] When the design parameters to be determined satisfy formula (7), the dynamic response determination method of the present invention realizes a consistent effective stiffness matrix; at this time, it not only reduces the computational cost of solving the folding rudder of the aircraft in the linear setting, but also realizes the optimal spectral characteristics;

[0088] γ = α = σ (7)

[0089] The dynamic response determination method of the present invention needs to achieve a consistent second - order accuracy in solving the folding rudder of the aircraft, that is, the parameters need to satisfy the conditions:

[0090]

[0091] Although the method of the present invention is directly self - starting, it can still output the derivative vector n+1 at time t and the accuracy is still second - order accuracy;

[0092] In order to effectively reduce the numerical instability caused by non - linearity, high - dimensionality, and spurious high - frequency components, the method of the present invention needs to have controllable numerical high - frequency dissipation ability, that is:

[0093]

[0094] where ρ ∞ ∈[0, 1] is a user - specified parameter; the strongest and weakest numerical high - frequency dissipations are introduced when ρ ∞ = 0 and ρ ∞ = 1 respectively; usually, the algorithm is called zero - dissipation when ρ ∞ = 1. When solving the folding rudder of the aircraft subsequently, it is usually required that the user selects a ρ ∞ value less than 1;

[0095] The method of the present invention not only requires that the numerical solution predicted by formula (6) is second - order accurate, but also requires that the solution predicted by formula (5) is also second - order accurate; therefore, the parameter μ should satisfy:

[0096]

[0097] To determine the redundant integration parameters, the present invention requires that the sizes of the first two sub - steps of the aircraft folding rudder dynamic response determination method be equal, that is:

[0098] ζ = 2γ. (11)

[0099] Solving equations (7) to (11) gives:

[0100]

[0101] At this time, only the parameter γ is unknown in the method of the present invention; the implicit method should generally achieve unconditional stability to ensure that a larger integration step size can be selected when solving the rigid problem of the aircraft folding rudder; the unconditional stability analysis shows that when the parameters γ and ρ ∞ satisfy the inequality constraint (13), the new implicit integration method achieves unconditional stability:

[0102] 24(ρ ∞ + 1)γ 4 - 8(ρ ∞ + 9)γ 3 + 48γ 2 - 12γ + 1 ≤ 0 (13)

[0103] After the user gives a value of ρ ∞ , the value range of the parameter γ, γ ∈ [γ min , γ max , can be obtained by solving equation (13) to achieve unconditional stability, as Figure 2 shown; note that when the parameter γ takes the value of γ min or γ max , the equality in inequality (13) holds;

[0104] S400. In the analysis of the aircraft folding rudder, the user should first set ρ ∞ ∈ [0, 1] to determine the numerical high - frequency dissipation, and then solve inequality (13) to obtain the value range of the parameter γ, γ ∈ [γ min , γ max and select an appropriate value of γ within this range; by default, the present invention recommends using γ = γ min ;

[0105] Using the amplitude and phase error calculation method, the main terms of the amplitude (δ) and phase (ε) errors of the new implicit integration method can be calculated as:

[0106]

[0107] where c4 is a certain complex expression and is not important, so it is not explicitly given here; ω is the natural frequency of the linear single-degree-of-freedom system; ξ represents the damping ratio of the linear single-degree-of-freedom system; O represents a higher-order quantity (mathematical usage); Equation (14) shows that when analyzing the undamped problem (ξ = 0) and the parameter takes the value γ = γ min or γ = γ max the amplitude error δ can greatly increase its order, reaching O(Dt 5 ), and thus reduces the amplitude error;

[0108]

[0109] When the parameter takes the value γ = γ min and γ = γ max the main term of the phase error of the new implicit integration algorithm is as Figure 3 shown; from Figure 3 it can be seen that the algorithm provides an obviously minimum phase error when taking γ = γ min and the phase error is very close to zero; therefore, the present invention recommends that users use γ = γ min to reduce the amplitude and phase errors;

[0110] In summary, the method for determining the dynamic response of the folding rudder of an aircraft can obtain a higher-order amplitude error order and the minimum phase error when the parameter takes the value of Equation (12) and the parameter γ takes the value of γ min .

[0111] Specific implementation plan two: A system for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error according to the present invention, the system has program modules corresponding to the above steps, and executes the steps in the method for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error when running.

[0112] Other combinations and connection relationships in this implementation plan are the same as those in the first specific implementation plan.

[0113] Specific implementation plan three: A computer-readable storage medium according to the present invention, the computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the method for determining the dynamic response of the folding rudder of an aircraft with high-order small amplitude error when called by a processor.

[0114] Other combinations and connection relationships in this implementation plan are the same as those in the first specific implementation plan.

[0115] Simulation experiment 1

[0116] Figure 4 shows the spectral characteristics of the method of the present invention under undamped conditions. As Figure 4 (a) shows, the user specifies the parameter ρ ∞ to precisely control the spectral radius in the high-frequency region (Ω ∈ [102 , 10 3 )'s spectral radius value, thus achieving a controllable numerical high-frequency dissipation ability. For example, the new algorithm of the present invention has a spectral radius that ultimately reaches zero when ρ ∞ = 0.0. Since Figure 4 the spectral radius shown in (a) is always less than or equal to 1, the new implicit algorithm of the present invention achieves unconditional stability. Figure 4 (b) shows that the numerical high-frequency dissipation increases gradually with the decrease of the parameter ρ ∞ , which also demonstrates the adjustable ability of the numerical high-frequency dissipation of the method of the present invention. Figure 4 (c) shows that the relative period error of the present method increases with the decrease of the parameter ρ ∞ , and the difference in relative period error between different parameters ρ ∞ is very small.

[0117] Figure 5 gives the convergence rate of the method of the present invention for solving the standard single-degree-of-freedom system. As can be seen from Figure 5 , the method of the present invention achieves consistent second-order accuracy and the difference between algorithms for different ρ ∞ values is small. Figure 6 The overshoot response of the method of the present invention is analyzed, where the two short light black dashed lines are the exact solution boundaries. If the numerical response exceeds the boundary of the exact solution, there is overshoot behavior. Obviously, the numerical solutions predicted by the new implicit integration algorithm do not exceed the exact solution boundary, so there is no overshoot behavior.

[0118] Table 1 compares the amplitude and phase errors of each implicit integration algorithm without damping. It can be found that the amplitude error of the conventional second-order implicit integration method only reaches O(Dt 3 ), while the amplitude error of the method of the present invention can reach O(Dt 5 ), which is two orders of magnitude higher than that of the conventional implicit algorithm. At the same time, although the phase errors reach the same order of magnitude, the new implicit integration method also gives a smaller leading error term, as shown in Figure 7 (b). Figure 7 shows that the method of the present invention is superior to the existing integration methods in both amplitude and phase errors; due to achieving higher-order small quantities, the leading amplitude error term of the method of the present invention is always zero on the scale of Δδ3 / ω 4 , demonstrating the advantages of the present invention over the existing implicit integration algorithms.

[0119] Table 1: Comparison of amplitude and phase errors of each implicit integration algorithm without damping (ξ = 0)

[0120]

[0121] a The parameter γ of the SUCI2 algorithm takes the value of b The parameter γ of the method of the present invention takes the value of γ min .

[0122] Simulation Experiment 2

[0123] (1) The present invention can be used to solve a first-order dynamic system in the following form:

[0124]

[0125] where t0 represents the initial time. Figure 8 The flowchart of the present invention for solving Equation (17) is given.

[0126] (2) The present invention can be used to solve a second-order dynamic system in the following form:

[0127]

[0128] where t0 represents the initial time. Figure 9 The flowchart of the present invention for solving Formula (18) is given.

[0129] (3) The present invention can be used to solve dynamic systems of any order. Its implementation scheme either reduces the high-order dynamic system to a first-order dynamic system through dimension elevation and order reduction for solution, or directly generalizes the implicit integration method, such as solving the second-order dynamic system (18).

[0130] (4) The present invention can be used to solve the (high-dimensional nonlinear) folding rudder system of an aircraft.

[0131] The present invention can be used to solve the folding rudder system of an aircraft with clearance nonlinearity as shown in Figure 10 to overcome the problem that the conventional time-domain integration method cannot stably solve this system. Figure 11 The nine-point acceleration response values when the method of the present invention solves this rudder system are given. As can be seen from the figure, the present invention can accurately and stably solve this clearance nonlinear folding rudder system, demonstrating its advantage in dynamic response solution.

[0132] Although the present invention is disclosed as above, the protection scope of the present invention is not limited thereto. Those skilled in the art of the present invention can make various changes and modifications without departing from the spirit and scope of the present invention, and these changes and modifications will all fall within the protection scope of the present invention.

[0133] References

[0134] [1] NEWMARK N M. A method of computation for structural dynamics[J]. Journal of Engineering Mechanic Division, 1959, 85(3): 67 - 94.

[0135] [2] HUGHES T J R. The finite element method: Linear static and dynamic finite element analysis[M]. Dover Publications, 2000.

[0136] [3] BATHE K J, BAIG M M I. On a composite implicit time integration procedure for nonlinear dynamics[J]. Computers&Structures, 2005, 83(31–32): 2513 - 2524.

[0137] [4] HILBER H M, HUGHES T J R, TAYLOR R L. Improved numerical dissipation for time integration algorithms in structural dynamics[J]. Earthquake Engineering&Structural Dynamics, 1977, 5(3): 283--292.

[0138] [5] WOOD W L, BOSSAK M, ZIENKIEWICZ O C. An alpha modification of Newmark’s method[J]. International Journal for Numerical Methods in Engineering, 1980, 15(10): 1562--1566.

[0139] [6]SHAO H,CAI C.A three parameters algorithm for numerical integration of structural dynamic equations[J].Chinese Journal of Applied Mechanics,1988,5(4):76-81.

[0140] [7]CHUNG J,HULBERT G M.A time integration algorithm for structural dynamics with improved numerical dissipation:The generalized-a method[J].Journal of Applied Mechanics,1993,60(2):371-375.

[0141] [8]YU K.A new family of generalized-α time integration algorithms without overshoot for structural dynamics[J].Earthquake Engineering & Structural Dynamics,2008,37(12):1389--1409.

[0142] [9]WILSON E L,FARHOOMAND I,BATHE K J.Nonlinear dynamic analysis of complex structures[J].Earthquake Engineering & Structural Dynamics,1972,1(3):241-252.

[0143]

[10] BAZZI G,ANDERHEGGEN E.The ρ-family of algorithms for time-step integration with improved numerical dissipation[J].Earthquake Engineering & Structural Dynamics,1982,10(4):537-550.

[0144]

[11] HOFF C, PAHL P J. Development of an implicit method with numerical dissipation from a generalized single-step algorithm for structural dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 67(3): 367-385.

[0145]

[12] ZIENKIEWICZ O C, WOOD W L, HINE N W, et al. A unified set of single-step algorithms. Part 1: General formulation and applications[J]. International Journal for Numerical Methods in Engineering, 1984, 20(8): 1529-1552.

[0146]

[13] CHUNG J, HULBERT G M. A family of single-step houbolt time integration algorithms for structural dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1–2): 1-11.

[0147]

[14] HILBER H M, HUGHES T J R. Collocation, dissipation and ‘overshoot’ for time integration schemes in structural dynamics[J]. Earthquake Engineering & Structural Dynamics, 1978, 6(1): 99-117.

[0148]

[15] Li Jinze. Design and analysis of self-starting step-by-step time integration algorithms for structural dynamics problems[D]. Harbin: Harbin Institute of Technology, 2023.

Claims

1. A method for determining the dynamic response of an aircraft folding rudder with a high-order small amplitude error, characterized in that: The following steps are involved: S100. The partial differential equation governing the dynamic response of the folding rudder structure of the aircraft degenerates into a set of ordinary differential equations of the following form after spatial discretization: Where M represents the global matrix; f represents the contribution of external and internal forces; u(t) represents the node state vector of the folding rudder of the aircraft; represents the time derivative of the node state vector; u0 represents the value of the node state vector at the initial time t0; S200, assuming that the integral time length of the folding rudder of the analyzed aircraft is t∈[0,T], first divide the integral interval into discrete time intervals, namely Where N represents the total number of integration steps, and then for the small integration interval t∈[t n ,t n+1 ] Construct the integral format and convert t n The dynamic response at time t is updated to n+1 moment, and finally, by recursively carrying out this process step by step, the dynamic response of the aircraft's folding rudder in the entire integral interval can be solved; S300, introducing constraints to determine the design parameters that appear in the process of solving the dynamic response of the aircraft folding rudder in the entire integration interval in step S200, so as to make the system have optimal spectral characteristics, consistent second-order accuracy, controllable numerical high-frequency dissipation capability and unconditional stability.

2. The method for determining the dynamic response of an aircraft folding rudder with a high-order small amplitude error according to claim 1, characterized in that: In step S200, it specifically includes: S210, the integral interval t∈[t n ,t n+1 ] is divided into three small intervals: t∈[t n ,t n +γDt]∪[t n +γDt,t n +ζDt]∪[t n +ζDt,t n+1 ](2)In the formula, γ and ζ represent the two sub-step size division rates of the aircraft folding rudder in the current integration interval, and satisfy γ≤ζ; Dt represents the integral step size for solving the aircraft folding rudder; S220, using t n The known node state vector u at time n Solving for t n+γ :=t n +γDt time node state quantity: S230, after completing step S220, using the known node state vector u n and the node derivative vector Solving for t n+ζ :=t n +ζDt node state quantity: S240, after completing steps S220 and S230, using the known node state vector u n and the node derivative vector Solving for t n+1 Node state quantity at the moment: S250, after completing the above calculation, using the known node state vector u n and the node derivative vector Calculate t n+1 Node state quantity at the moment: By performing iterative calculation step by step, the structural dynamic response of the aircraft folding rudder in the entire integration interval is obtained.

3. The method for determining the dynamic response of an aircraft folding rudder with a high-order small amplitude error according to claim 2, characterized in that: In step S300, it specifically includes: introducing constraint conditions for the design parameters γ, ζ, α, μ, σ, β1 and β2 in step S200 to determine, When the design parameters satisfy formula (7), the method for determining the dynamic response of the aircraft folding rudder realizes a consistent effective stiffness matrix, which not only reduces the computational cost of solving the linear setting of the aircraft folding rudder but also realizes the optimal spectral characteristics. γ=α=σ (7) The parameters to be designed can achieve consistent second-order accuracy in solving the folding rudder of the aircraft by satisfying the following conditions: In order to effectively reduce numerical instability, it is necessary to have controllable numerical high-frequency dissipation capability, namely: In the formula, ρ ∞ ∈[0,1] represents the user-specified parameter; ∞ =0 and ρ ∞ = 1, the strongest and weakest numerical high-frequency dissipation is introduced. ∞ = 1 is called zero dissipation. When solving the folding rudder of the aircraft in the subsequent process, the user is usually required to select a ρ less than 1. ∞ value; The parameter μ satisfies the following condition, so that the numerical solution predicted by formula (6) and the solution predicted by formula (5) are is second-order accurate: In order to determine the redundant integral parameters, the first two sub-steps of the method for determining the folding rudder dynamic response of the aircraft are required to be equal in size, that is: ζ=2γ. (11) Solving formula (7) to formula (11) yields: When the parameters γ and ρ ∞ When the inequality constraint (13) is satisfied, the dynamic response determination method has unconditional stability: 24(p ∞ +1)c 4 -8(r ∞ +9)c 3 +48c 2 -12γ+1≤0 (13) The user has a given ∞ After the value is calculated, the value interval of parameter γ is obtained by solving formula (13) γ∈[γ min ,γ max ] to achieve unconditional stability.

4. The method for determining the dynamic response of an aircraft folding rudder with a high-order small amplitude error according to claim 3, characterized in that: First, set ρ ∞ ∈[0,1] to determine the numerical high-frequency dissipation, and then solve inequality (13) to obtain the value interval of parameter γ ∈ [γ min ,γ max ], select the value of γ within this interval by the following steps; The amplitude and phase error calculation method is used to calculate the main terms of amplitude (δ) and phase (ε) errors: Where c4 is the irrelevant parameter value; ω is the natural frequency of the linear single degree of freedom system; ξ represents the damping rate of the linear single degree of freedom system; O represents the high-order quantity; Formula (14) shows that when analyzing the undamped problem and the parameter value is γ = γ min Or γ=γ max When the amplitude error δ is increased in order, it reaches O(Dt 5 ), and thus reduce the amplitude error to: The method for determining the dynamic response of the aircraft folding rudder is as follows: the parameter value is formula (12) and the parameter γ is γ min A higher order amplitude error and a minimum phase error are obtained.

5. The method for determining the dynamic response of an aircraft folding rudder with a high-order small amplitude error according to claim 4, characterized in that: Nine-point acceleration response values ​​of the aircraft folding rudder system used to solve the nonlinearity of clearance.

6. A high-order small amplitude error aircraft folding rudder dynamic response determination system, characterized by: The system has a program module corresponding to the steps of any one of claims 1 to 5, and executes the steps in the method for determining the folding rudder dynamic response of an aircraft with a high-order small amplitude error during operation.

7. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and the computer program is configured to implement the steps of the method for determining the folding rudder dynamic response of an aircraft with a high-order small amplitude error according to any one of claims 1 to 5 when called by a processor.

Images