CN106844914A - A kind of rapid simulation method of re-entry space vehicle wing flutter response - Google Patents

Abstract

The invention discloses a kind of Fast simulation algorithm of re-entry space vehicle wing flutter response, belong to field of aerospace technology, be related to a kind of Vibration Simulation technology of aeroelasticity structure.By to re-entry space vehicle, wing fluid-solid coupled motion problem in flow field is modeled during high-speed flight, and proposes that new differential equation method is solved.Quantitatively study the dynamic characteristics such as flutter, the limit cycle of re-entry space vehicle binary wing in High Speed Flow Field.For accurate reliable Mathematical Modeling is set up in the design of re-entry space vehicle wing, and a kind of new fast and effectively time discrete method is applied in the research of the kinetic model of non-linear binary wing, and then studies the various dynamic characteristics of wing.The method can be emulated suitable for strong nonlinearity binary wing flutter, with efficient, high-precision feature.

Description

A kind of rapid simulation method of re-entry space vehicle wing flutter response

Technical field

The invention belongs to field of aerospace technology, it is related to a kind of air-foil aeroelastic vibration emulation technology, specifically It is related to a kind of rapid simulation method of re-entry space vehicle wing flutter response.

Background technology

Re-entry space vehicle is a kind of reusable world shuttle system of new generation, is had for reducing space transoportation cost There is important meaning.One of its feature is freely to rise and fall and fly between the outer space and atmosphere.Develop empty day flight Device, will solve following key technology:Have super light weight, tough, heat-resisting and shock resistance empty day device structure skill high concurrently Art, hypersonic technology, Highly maneuverable aircraft technology, empty day airmanship, stealth technique high, precisely hit and reliability over long distances Technology.In these key technologies, Aeroelastic Problems are the underlying issues for designing, manufacturing re-entry space vehicle, and with the problem The design of closely related aircraft wing is then current study hotspot and technological difficulties.

Vibration mode of the wing of re-entry space vehicle during high-speed flight directly affects the intensity of wing structure body.One As for, wing vibration mode that may be present includes：Stable state, limit-cycle oscillation, periodic vibration, chaotic motion and diverging. The design of wing is tried one's best makes it that stable state is in whole flight course.Strength damage so to wing structure body is most Small.Limit-cycle oscillation refers to that wing does class simple harmonic oscillation.And periodic vibration refers to the somewhat complicated still vibration of vibration of wing Amplitude has repeatability.Under different airfoil profiles parameter and flying speed, flutter of aerofoil border and vibration mode are different.Therefore have Necessity is analyzed to airfoil flutter system.Original wing model is rich in integral term.Because the presence of integral term, either right For direct numerical solution algorithm or perturbation method or weighted residual method, solve get up to be required for carrying out in advance it is substantial amounts of Symbol processing, consumes substantial amounts of computing resource.Existing method such as harmonic wave equilibrium method, perturbation method etc. tend not to quick, comprehensive Analysis system response.This largely have impact on the parameter designing of wing structure.

The content of the invention

It is an object of the invention to provide a kind of rapid simulation method of re-entry space vehicle wing flutter response, on overcoming The defect of prior art presence is stated, the present invention is based on time discrete method, by the discrete acquisition aircraft wing gas in time-domain The algebraic equation of bullet system, it is to avoid substantial amounts of symbolic operation in general nonlinear system analysis method, while simplify being System equation so that solution efficiency obtains larger raising.

To reach above-mentioned purpose, the present invention is adopted the following technical scheme that：

A kind of rapid simulation method of re-entry space vehicle wing flutter response, comprises the following steps：

1) binary flutter of aerofoil kinetic model of the re-entry space vehicle under fluid structurecoupling is set up；

2) by carrying out solution setup time discrete method Algebraic Equation set to flutter of aerofoil kinetic model；

3) the time discrete method Algebraic Equation set of compact is set up；

4) above-mentioned compact time discrete method Algebraic Equation set is solved, the vibration response curve of aerofoil system is obtained.

Further, step 1) it is specially：

Contain aeroelastic model during Hookean spring with reference to binary wing, and consider aerofoil system in pitching and sink-float two The structural nonlinear of the individual free degree, sets up aerofoil system equation：

Wherein, x_αIt is dimensionless distance of the wing gas bullet coordinate origin to barycenter, ξ=h/b is dimensionless sink-float amount, () represents the derivative to nondimensional time τ, and τ=Ut/b, t are time, U^*It is nondimensional velocity, is defined as U^*=U/ (b ω_α), U is air speed of incoming flow,Wherein ω_ξAnd ω_αIt is respectively that coupled wave equation sink-float is not intrinsic with pitch freedom Frequency, ζ_ξAnd ζ_αIt is damping ratio, r_αIt is the torque around elastic shaft, α is the angle of pitch, and h is deflection angle, μ=m/ π ρ b², m is wing matter Amount, ρ is atmospheric density, and b is half chord length of wing；

M (α) and G (ξ) are respectively the nonlinear terms of pitching and the sink-float free degree, and expression formula is：

M (α)=alpha+beta α³, G (ξ)=ξ+γ ξ³, (3)

Wherein β and γ is nonlinear terms coefficient；

C_L(τ) and C_M(τ) is linear aerodynamic force and aerodynamic moment, and expression formula is：

Wherein Wagner function phis (τ) areψ₁, ψ₂, ε₁, ε₂It is Wagner constants, a_hIt is dimensionless distance of the wing axis to the gas bullet origin of coordinates；

Then, by introducing one group of integral transform, by C in above formula_LAnd C_MComprising integral term eliminate so that will integration The differential equation is converted into differential equation group, is designated as following form：

Wherein, c₀=1+1/ μ, c₁=x_α-a_h/ μ,

c₃=[1+ (1-2a_h)(1-ψ₁-ψ₂)]/μ, c₄=2 (ε₁ψ₁+ε₂ψ₂)/μ,

c₅=2 [1- ψ₁-ψ₂+(1/2-a_h)(ε₁ψ₁+ε₂ψ₂)]/μ,

c₆=2 ε₁ψ₁[1-ε₁(1/2-a_h)]/μ, c₇=2 ε₂ψ₂[1-ε₂(1/2-a_h)]/μ,

Wherein, r_αIt is the radius of gyration of elastic coordinate system.

Further, step 2) it is specially：

By step 1) 6 α (τ) to be found a function in the differential equation group (4) that obtains, ξ (τ), ω_i(τ), i=1,2,3, 4, assume initially that into Fourier progression forms：

Wherein, α_j, ξ_j, ω_j, j=0 ..., 2N is corresponding harmonic constant；

The approximate solution (5) that it will be assumed is substituted into differential equation group (4), obtains residual error function：

Wherein, R_jRepresent

Finally, residual error function R is forced_j2N+1 equidistant time point τ over one period_iOn be zero obtain final product to the time from Arching pushing Algebraic Equation set, the aerofoil system contains 6 × (2N+1) individual equations ,+1 unknown number of 6 × (2N+1).

Further, step 3) it is specially：

Time discrete treatment is carried out to time discrete method Algebraic Equation set each, to α (τ) in 2N+1 time point τ_iFrom Dissipate, can obtain：

Wherein θ_i=ω τ_i, above formula is written as matrix form：

It is represented simply as：

Wherein,

Q_α=[α₀,α₁,...,α_2N]^T

SimilarWhereinQ_ξIt is respectively relevant with ξ With a harmonic coefficient,WithRespectively and ω_iIt is relevant with a harmonic coefficient；

Then, it is rightDiscrete obtaining is carried out at 2N+1 time point：

Above formula is written as matrix form：

Order matrix

Therefore have：

Wherein

Similarly have：WhereinWithPoint It is not on ξ and ω_jDerivative with point；

Then, it is rightIt is discrete at 2N+1 time point, have：

Above formula (9) is written as matrix form：

Square formation FA in above formula²It is expressed as

Therefore,

Similarly haveWhereinWithIt is respectively on ξ and ω_jSecond dervative with point；

Thus following form is obtained after time discrete method Algebraic Equation set is transformed：

Wherein D=FAF^-1, and

WillWithWithRepresent, then they are updated in nonlinear equation and are obtained：

Wherein A_1α, B_1ξ, A_2α, B_2ξRespectively：

A_1α=c₁ω²D²+c₃ωD+c₅I+c₆(ωD+∈₁I)^-1+c₇(ωD+∈₂I)^-1

B_1ξ=c₀ω²D²+c₂ωD+(c₄+c₁₀)I+c₈(ωD+∈₁I)^-1+c₉(ωD+∈₂I)^-1

A_2α=d₁ω²D²+d₃ωD+d₅I+d₆(ωD+∈₁I)^-1+d₇(ωD+∈₂I)^-1

B_2ξ=d₀ω²D²+d₂ωD+d₄I+d₈(ωD+∈₁I)^-1+d₉(ωD+∈₂I)^-1.

Formula (11) is the compact time discrete method Algebraic Equation set with time domain variable as variable, when using-frequency conversion pass It is formulaWithAbove formula is converted to the compact time discrete with frequency variable as variable Method Algebraic Equation set (12)：

Formula (12) containing 2N+2 unknown number, wherein, A₂=A_2α, A₁=A_1α, B₂=B_2ξ, B₁=B_1ξ。

Further, step 4) the middle use scalar Homotopy above-mentioned compact time discrete method Algebraic Equation set of solution.

Compared with prior art, the present invention has following beneficial technique effect：

The inventive method is based on time discrete method, by the generation of the discrete acquisition aircraft wing gas bullet system in time-domain Number equations, it is to avoid substantial amounts of symbolic operation in general nonlinear system analysis method, while simplifying system equation so that Solution efficiency obtains larger raising, coordinates by by time discrete scheme and parameter scanning method, successfully avoid general analysis side The non-physical solution being likely to occur in method, improves the reliability of wing flutter emulation, and the steady-state period of aerofoil system can be rung Direct estimation should be carried out, and without being simulated to transient oscillation, specific aim is stronger, computing resource is saved, while also allowing for The situation of change of analysis system dynamic response under variable element.

Brief description of the drawings

Fig. 1 is the flow chart of the inventive method；

Fig. 2 is binary wing model schematic diagram；

Fig. 3 is the wing flutter frequency and flying speed response curve calculated using low order time discrete method, wherein a) is The result obtained using 9 rank time domain point collocations, b) is the result that is obtained using 11 rank time domain point collocations；

Fig. 4 is the wing flutter frequency and flying speed response curve calculated using high order time discrete method, wherein a) is The result obtained using 15 rank time domain point collocations, b) is the result that is obtained using 21 rank time domain point collocations.

Specific embodiment

The present invention is described in further detail below：

A kind of rapid simulation method of the re-entry space vehicle wing flutter response based on time discrete method, the spy of the method Point is, using Integral Transformation Method, re-entry space vehicle air-foil flutter Integral-differential Equations to be converted into pure differential side Journey.Further, choose the Fourier space containing unknowm coefficient and respond approximate solution as flutter, bring the approximate solution into flutter micro- Point equation, then use time discrete method by the differential equation it is discrete be using Fourier coefficient and vibration frequency as unknown variable Algebraic Equation set.Finally, coefficient to be asked is solved using scalar Homotopy Method, so that it is determined that the cycle of re-entry space vehicle flutter of aerofoil Solution.The semi analytic periodic solution of the explanation flutter system that the method is given, the method has computational efficiency and precision very high.

Specifically include following steps：

1) re-entry space vehicle flutter of aerofoil kinetic model is set up：According to Lagrange's equation, it is considered to structural nonlinear and Piston theory, sets up binary wing flutter original integral differential equation of the re-entry space vehicle under fluid structurecoupling, and carry out immeasurable Guiding principle treatment.One group of integral transform is introduced, the integral term in integral differential equation is eliminated, so as to obtain pure differential equation form Kinetic model.

2) use time discrete method carries out sliding-model control to flutter of aerofoil kinetics equation：Treating in wing model is asked Function, i.e. the wing angle of pitch and sink-float displacement, and the form of Fourier space is assumed to be by the variable that integral transform is introduced, And the coefficient to be asked of Fourier space is assembled into vector.The approximate solution that it will be assumed is substituted into above-mentioned differential equation group, obtains residual Difference function.Then, it is over one period zero Algebraic Equation set for obtaining time discrete method on limited point to force residual error function.

3) the time discrete method Algebraic Equation set of compact is set up：Because the equation array dimension set up in step 2 is too big, no Beneficial to calculating, therefore enter row equivalent conversion to it, to obtain the algebraic equation of compact form.Using assume function for Fourier level Several this feature, the matrix relationship set up between discrete point and harmonic constant, then by matrixing, setup time single order is led Number and the relation between the discrete point and harmonic constant of second dervative.Based on matrixing relation, what residual error function was constituted is non- System of linear equations is processed item by item, and the linear segment in gained equation group is brought into non-linear partial, it becomes possible to Simple, the relatively low compact time discrete method Algebraic Equation set of dimension to form.

4) the time discrete method Algebraic Equation set of above-mentioned compact is solved using scalar homotopy Method：Time discrete method is to first Value is very sensitive, and arbitrarily given initial value typically cannot be solved truly, therefore, it is offer with the insensitive scalar Homotopy Method of initial value Reasonable initial value obtains semi analytic periodic solution.

In order to objects and advantages of the present invention are better described, present invention is done into one with example below in conjunction with the accompanying drawings Step explanation：

1. the Mathematical Modeling of binary wing flutter model is set up：

Fig. 2 is the binary wing flutter model containing pitching and sink-float two-freedom.Sink-float is represented with h, is downwards pros To.Pitching on elastic shaft is represented with α, past to face upward as just.Distance of the elastic shaft away from airfoil center is a_hB, barycenter is from elasticity The distance of axle is x_ab；Two positive directions of distance point to trailing edge.

Consider that pitching is with the Dimensionless Form system equation of the structural nonlinear of sink-float:

Wherein x_αIt is dimensionless distance of the wing gas bullet coordinate origin to barycenter, ξ=h/b is dimensionless sink-float amount, () represents the derivative to nondimensional time τ, and wherein τ=Ut/b, t are time, U^*It is nondimensional velocity, is defined as U^*=U/ (b ω_α), U is air speed of incoming flow.Wherein ω_ξAnd ω_αIt is respectively that coupled wave equation does not rise and fall and pitch freedom Intrinsic frequency.ζ_ξAnd ζ_αIt is damping ratio, r_αIt is the torque around elastic shaft, α is the angle of pitch, and h is deflection angle, μ=m/ π ρ b², m is machine Wing quality, ρ is atmospheric density, and b is half chord length of wing.M (α) and G (ξ) are respectively pitching and the non-linear of the free degree that rise and fall , their expression is：

M (α)=alpha+beta α³, G (ξ)=ξ+γ ξ³, (3)

Wherein β and γ is nonlinear terms coefficient.C_L(τ) and C_M(τ) is linear aerodynamic force and aerodynamic moment：

Wherein Wagner function phis (τ) areψ₁, ψ₂, ε₁, ε₂It is Wagner constants, a_hIt is dimensionless distance of the wing axis to the gas bullet origin of coordinates.For the ease of model solution, we introduce Following four integral transformation relational expression, full scale equation group can equivalence transformation be following form：

Wherein：

c₀=1+1/ μ, c₁=x_α-a_h/ μ,

c₃=[1+ (1-2a_h)(1-ψ₁-ψ₂)]/μ, c₄=2 (ε₁ψ₁+ε₂ψ₂)/μ,

c₅=2 [1- ψ₁-ψ₂+(1/2-a_h)(ε₁ψ₁+ε₂ψ₂)]/μ,

c₆=2 ε₁ψ₁[1-ε₁(1/2-a_h)]/μ, c₇=2 ε₂ψ₂[1-ε₂(1/2-a_h)]/μ,

Wherein, r_αIt is the radius of gyration of elastic coordinate system.

2. the time discrete method of binary wing model solves scheme：

By 6 α (τ) to be found a function in equation group (4), ξ (τ), ω_i(τ), i=1,2,3,4 assume initially that into Fourier progression forms.

The approximate solution (5) that it will be assumed is substituted into math equation (4), obtains residual error function：

Wherein R_jRepresentNow, residual error letter is forced Number R_j2N+1 equidistant time point τ over one period_iOn be that zero can be obtained by time discrete method Algebraic Equation set, this is System contains the individual equations 6 of 6 × (2N+1) × unknown numbers of (2N+1)+1.

3. the time discrete method Algebraic Equation set of compact is set up：

Time discrete treatment is carried out to time discrete method Algebraic Equation set each.To α (τ) in 2N+1 time point τ_iFrom Dissipate, can obtain：

Wherein θ_i=ω τ_i, above formula can be written as matrix form：

Also can be represented simply as,

Wherein

Q_α=[α₀,α₁,...,α_2N]^T

Similar hasWhereinQ_ξRespectively have with ξ Close with a harmonic coefficient,WithRespectively and ω_iIt is relevant with a harmonic coefficient.

Then, it is rightDiscrete obtaining is carried out at 2N+1 time point：

Above formula is easily written as matrix form：

Order matrix

Therefore have：

Wherein

Identical reason,WhereinWithIt is respectively on ξ and ω_jDerivative with point.

Then, it is rightIt is discrete at 2N+1 time point, have：

The matrix form of equation (9) is：

Square formation in above formula can be expressed as with FA2

Therefore,

SimilarWhereinWithIt is respectively on ξ and ω_jSecond dervative with point.

Above, we are processed equation item by item, have obtained the matrix form after time discrete.Time discrete method generation Following form is obtained after number equation is transformed：

Wherein D=FAF^-1, and

It is nonlinear due to there was only second equation, other five equations are linear equation, therefore equation group (10) Can further simplify.WillWithWithRepresent, then they are updated in nonlinear equation and are obtained：

Wherein A_1α, B_1ξ, A_2α, B_2ξFor：

A_1α=c₁ω²D²+c₃ωD+c₅I+c₆(ωD+∈₁I)^-1+c₇(ωD+∈₂I)^-1

B_1ξ=c₀ω²D²+c₂ωD+(c₄+c₁₀)I+c₈(ωD+∈₁I)^-1+c₉(ωD+∈₂I)^-1

A_2α=d₁ω²D²+d₃ωD+d₅I+d₆(ωD+∈₁I)^-1+d₇(ωD+∈₂I)^-1

B_2ξ=d₀ω²D²+d₂ωD+d₄I+d₈(ωD+∈₁I)^-1+d₉(ωD+∈₂I)^-1.

System (11) is the compact time discrete method Algebraic Equation set with time domain variable as variable.

When using-frequency conversion relational expressionWithSystem (11) is transformed to：

Equation (12) is the compact time discrete method Algebraic Equation set with frequency variable as variable, unknown containing 2N+2 Number, wherein, A₂=A_2α, A₁=A_1α, B₂=B_2ξ, B₁=B_1ξ。

4. the time discrete method Algebraic Equation set of above-mentioned compact is solved using scalar homotopy Method, and incorporating parametric is scanned Method solves the time discrete method Algebraic Equation set of above-mentioned compact：

Time discrete method is very sensitive to initial value, and arbitrarily given initial value typically cannot be solved truly, therefore, with initial value not Sensitive scalar Homotopy Method obtains semi analytic periodic solution to provide reasonable initial value.In order to carry out System Parameter Design, it is necessary to ask for The amplitude frequency curve of response.There is provided reasonable initial value to obtain amplitude-frequency response for time discrete method used here as parameter scanning method. Fig. 3 and Fig. 4 are the response curves on fundamental vibration frequency than flying speed that time discrete method and RK4 are calculated.Fig. 3 (a) shows Show, for forward scan,Before, the result of TDC9 and RK4 is coincide, and is started more than later TDC9 and RK4 Separate.Also, the response curve of TDC9 is consecutive variations, that is to say, that TDC9 can not predict binary wing flutter system Secondary bifurcation.For reverse scanning, it can be seen that TDC9 is consistent in whole lower branch curve with RK4.The bifurcation point predicted using RK4 It is 1.84, the response curve more than RK4 after the bifurcation point jumps to upper branch from lower branch.For TDC9, before bifurcation point TDC9 can obtain entirely descending branch response curve.But after surmounting bifurcation point, TDC9 response curves do not jump to branch response Curve, but cranky some other curve is obtained near lower branch.Because reaching parameter scanning method after bifurcation point Rational initial value, scanning method failure cannot be provided.Because TDC methods are to the extremely sensitive property of initial value, therefore TDC9 is in bifurcation point Result afterwards corresponds to non-physical solution.

Fig. 3 (b) is the frequency that TDC11 is calculated -- velocity-response curve figure.During forward scan, the result and mark of TDC11 Quasi- solution until 2.27 be all very close to, this is better than the result of TDC9.Deviateing occurs in the result of TDC11 of continuing to push the speed, This means non-physical solution is occurred in that.TDC11 can not detect secondary bifurcation point as TDC9,.During reverse scanning, TDC11 can obtain the lower branch of whole response curve, and it is more smaller than standard value 1.84 to have obtained secondary Hopf forks value.Due to Bifurcation point makes scanning method fail, therefore be given more than TDC11 after bifurcation point is non-physical solution.

TDC15 can accurately depict entirely lower branch response curve when Fig. 4 (a) shows reverse scanning, but forward direction is swept Whole upper branch can not be given when retouching.The result of TDC21 is shown in Fig. 4 (b).As illustrated, forward scan has obtained entirely going up branch curve, And the bifurcation point predicted is near the mark value 2.34 for 2.38.After surmounting bifurcation point, non-physical solution is inevitably occurred in that.It is inverse It is very good to be coincide with standard solution to the result of scanning, and entirely lower branch height of curve is consistent, and bifurcation point is also accurately obtained.Surmount fork After point, non-physical solution branch occurs, untilJust the upper of physical significance is jumped to from non-physical solution branch Branch curve.

In a word, use time discrete method can obtain Kind of Nonlinear Dynamical System with the method that parameter scanning method is combined Various response curves.Also, with the increase of harmonic wave number, computational accuracy is higher.

Claims (5)

1. the rapid simulation method that a kind of re-entry space vehicle wing flutter is responded, it is characterised in that comprise the following steps：

1) binary flutter of aerofoil kinetic model of the re-entry space vehicle under fluid structurecoupling is set up；

2) by carrying out solution setup time discrete method Algebraic Equation set to flutter of aerofoil kinetic model；

3) the time discrete method Algebraic Equation set of compact is set up；

4) above-mentioned compact time discrete method Algebraic Equation set is solved, the vibration response curve of aerofoil system is obtained.

2. the rapid simulation method that a kind of re-entry space vehicle wing flutter according to claim 1 is responded, it is characterised in that Step 1) it is specially：

Contain aeroelastic model during Hookean spring with reference to binary wing, and consider aerofoil system in pitching and sink-float two certainly By the structural nonlinear spent, aerofoil system equation is set up：

$\stackrel{\·\·}{\mathrm{\ξ}}+x_{\mathrm{\α}}\stackrel{\·\·}{\mathrm{\α}}+2{\mathrm{\ζ}}_{\mathrm{\ξ}}\frac{\stackrel{\‾}{\mathrm{\ω}}}{U*}\stackrel{\·}{\mathrm{\ξ}}+{\left(\frac{\stackrel{\‾}{\mathrm{\ω}}}{U*}\right)}^{2}G\left(\mathrm{\ξ}\right)=-\frac{1}{\mathrm{\π}\mathrm{\μ}}{C}_L\left(\mathrm{\τ}\right),---\left(1\right)$

$\frac{x_{\mathrm{\α}}}{r_{\mathrm{\α}}^2}\stackrel{\·\·}{\mathrm{\ξ}}+\stackrel{\·\·}{\mathrm{\α}}+2{\mathrm{\ζ}}_{\mathrm{\α}}\frac{1}{U*}\stackrel{\·}{\mathrm{\α}}+{\left(\frac{1}{U*}\right)}^{2}M\left(\mathrm{\α}\right)=\frac{2}{{\mathrm{\π\μr}}_{\mathrm{\α}}^2}{C}_M\left(\mathrm{\τ}\right),---\left(2\right)$

Wherein, x_αIt is dimensionless distance of the wing gas bullet coordinate origin to barycenter, ξ=h/b is dimensionless sink-float amount, () table Show the derivative to nondimensional time τ, τ=Ut/b, t are time, U^*It is nondimensional velocity, is defined as U^*=U/ (b ω_α), U is sky Gas speed of incoming flow,Wherein ω_ξAnd ω_αIt is respectively the intrinsic frequency of not coupled wave equation sink-float and pitch freedom, ζ_ξAnd ζ_αIt is damping ratio, r_αIt is the torque around elastic shaft, α is the angle of pitch, and h is deflection angle, μ=m/ π ρ b², m is wing quality, ρ It is atmospheric density, b is half chord length of wing；

M (α) and G (ξ) are respectively the nonlinear terms of pitching and the sink-float free degree, and expression formula is：

M (α)=alpha+beta α³, G (ξ)=ξ+γ ξ³, (3)

Wherein β and γ is nonlinear terms coefficient；

C_L(τ) and C_M(τ) is linear aerodynamic force and aerodynamic moment, and expression formula is：

$\begin{array}{c}{C}_L\left(\mathrm{\τ}\right)=\mathrm{\π}(\stackrel{\·\·}{\mathrm{\ξ}}-a_h\stackrel{\·\·}{\mathrm{\α}}+\stackrel{\·}{\mathrm{\α}})+2\mathrm{\π}\[\mathrm{\α}\left(0\right)+\stackrel{\·}{\mathrm{\ξ}}\left(0\right)+(\frac{1}{2}-a_h)\stackrel{\·}{\mathrm{\α}}\left(0\right)\]\mathrm{\φ}\left(\mathrm{\τ}\right)\\ +2\mathrm{\π}{\∫}_0^{\mathrm{\τ}}\mathrm{\φ}\left(\mathrm{\τ}-\mathrm{\σ}\right)\[\stackrel{\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)+\stackrel{\·\·}{\mathrm{\ξ}}\left(\mathrm{\σ}\right)+(\frac{1}{2}-a_h)\stackrel{\·\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)\]d\mathrm{\σ},\end{array}$

$\begin{array}{c}{C}_M\left(\mathrm{\τ}\right)=\mathrm{\π}(\frac{1}{2}+a_h)\[\mathrm{\α}\left(0\right)+\stackrel{\·}{\mathrm{\ξ}}\left(0\right)+(\frac{1}{2}-a_h)\stackrel{\·}{\mathrm{\α}}\left(0\right)\]\mathrm{\φ}\left(\mathrm{\τ}\right)+\frac{\mathrm{\π}}{2}(\stackrel{\·\·}{\mathrm{\ξ}}-a_h\stackrel{\·\·}{\mathrm{\α}})-\frac{\mathrm{\π}}{16}\stackrel{\·\·}{\mathrm{\α}}\\ -\left(\frac{1}{2}-a_h\right)\frac{\mathrm{\π}}{2}\stackrel{\·}{\mathrm{\α}}+\mathrm{\π}\left(\frac{1}{2}+a_h\right){\∫}_0^{\mathrm{\τ}}\mathrm{\φ}\left(\mathrm{\τ}-\mathrm{\σ}\right)\[\stackrel{\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)+\stackrel{\·\·}{\mathrm{\ξ}}\left(\mathrm{\σ}\right)+(\frac{1}{2}-a_h)\stackrel{\·\·}{\mathrm{\α}}\left(\mathrm{\σ}\right)\]d\mathrm{\σ},\end{array}$

Wherein Wagner function phis (τ) areψ₁, ψ₂, ε₁, ε₂It is Wagner normal Number, a_hIt is dimensionless distance of the wing axis to the gas bullet origin of coordinates；

Then, by introducing one group of integral transform, by C in above formula_LAnd C_MComprising integral term eliminate, so as to by integral differential Equations turned is differential equation group, is designated as following form：

$c_0\stackrel{\·\·}{\mathrm{\ξ}}+c_1\stackrel{\·\·}{\mathrm{\α}}+c_2\stackrel{\·}{\mathrm{\ξ}}+c_3\stackrel{\·}{\mathrm{\α}}+c_4\mathrm{\ξ}+c_5\mathrm{\α}+c_6{w}_1+c_7{w}_2+c_8{w}_3+c_9{w}_4+c_{10}G\left(\mathrm{\ξ}\right)=f\left(\mathrm{\τ}\right),$

$d_0\stackrel{\·\·}{\mathrm{\ξ}}+d_1\stackrel{\·\·}{\mathrm{\α}}+d_2\stackrel{\·}{\mathrm{\ξ}}+d_3\stackrel{\·}{\mathrm{\α}}+d_4\mathrm{\ξ}+d_5\mathrm{\α}+d_6{w}_1+d_7{w}_2+d_8{w}_3+d_9{w}_4+d_{10}M\left(\mathrm{\α}\right)=g\left(\mathrm{\τ}\right),$

${\stackrel{\·}{w}}_1=\mathrm{\α}-{\∈}_1{w}_1,$

${\stackrel{\·}{w}}_2=\mathrm{\α}-{\∈}_2{w}_2,$

${\stackrel{\·}{w}}_3=\mathrm{\ξ}-{\∈}_1{w}_3,$

${\stackrel{\·}{w}}_4=\mathrm{\ξ}-{\∈}_2{w}_4.---\left(4\right)$

Wherein, c₀=1+1/ μ, c₁=x_α-a_h/ μ,

c₃=[1+ (1-2a_h)(1-ψ₁-ψ₂)]/μ, c₄=2 (ε₁ψ₁+ε₂ψ₂)/μ,

c₅=2 [1- ψ₁-ψ₂+(1/2-a_h)(ε₁ψ₁+ε₂ψ₂)]/μ,

c₆=2 ε₁ψ₁[1-ε₁(1/2-a_h)]/μ, c₇=2 ε₂ψ₂[1-ε₂(1/2-a_h)]/μ,

$c_8=-2{\mathrm{\ϵ}}_1^2{\mathrm{\ψ}}_1,c_9=-2{\mathrm{\ϵ}}_2^2{\mathrm{\ψ}}_2/\mathrm{\μ},c_{10}={(\stackrel{\‾}{\mathrm{\ω}}/U^{*})}^2,$

$d_0=x_{\mathrm{\α}}/r_{\mathrm{\α}}^2-a_h/{\mathrm{\μr}}_{\mathrm{\α}}^2,d_1=1+(1+8a_h^2)/\left(8{\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_2=-(1+2a_h)(1-{\mathrm{\ψ}}_1-{\mathrm{\ψ}}_2)/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_4=-(1+2a_h)({\mathrm{\ϵ}}_1{\mathrm{\ψ}}_1+{\mathrm{\ϵ}}_2{\mathrm{\ψ}}_2)/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_5=-(1+2a_h)(1-{\mathrm{\ψ}}_1-{\mathrm{\ψ}}_2)/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right)-(1-4a_h^2)({\mathrm{\ϵ}}_1{\mathrm{\ψ}}_1+{\mathrm{\ϵ}}_2{\mathrm{\ψ}}_2)/\left(2{\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_6=-(1+2a_h){\mathrm{\ψ}}_1{\mathrm{\ϵ}}_1\[1-{\mathrm{\ϵ}}_1(1/2-a_h)\]/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_7=-(1+2a_h){\mathrm{\ψ}}_2{\mathrm{\ϵ}}_2\[1-{\mathrm{\ϵ}}_2(1/2-a_h)\]/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),$

$d_8=(1+2a_h){\mathrm{\ψ}}_1{\mathrm{\ϵ}}_1^2/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),d_9=(1+2a_h){\mathrm{\ψ}}_2{\mathrm{\ϵ}}_2^2/\left({\mathrm{\μr}}_{\mathrm{\α}}^2\right),d_{10}={(1/U^{*})}^2,$

${\mathrm{\ω}}_1={\∫}_0^{\mathrm{\τ}}{e}^{-c_1(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\α}\left(\mathrm{\σ}\right)d\mathrm{\σ},{\mathrm{\ω}}_2={\∫}_0^{\mathrm{\τ}}{e}^{-c_2(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\α}\left(\mathrm{\σ}\right)d\mathrm{\σ},$

$\begin{array}{cc}{\mathrm{\ω}}_3={\∫}_0^{\mathrm{\τ}}{e}^{-c_1(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\ξ}\left(\mathrm{\σ}\right)d\mathrm{\σ}& {\mathrm{\ω}}_4={\∫}_0^{\mathrm{\τ}}{e}^{-c_2(\mathrm{\τ}-\mathrm{\σ})}\mathrm{\ξ}\left(\mathrm{\σ}\right)d\mathrm{\σ},\end{array}$

Wherein, r_αIt is the radius of gyration of elastic coordinate system.

3. the rapid simulation method that a kind of re-entry space vehicle wing flutter according to claim 2 is responded, it is characterised in that Step 2) it is specially：

By step 1) 6 α (τ) to be found a function in the differential equation group (4) that obtains, ξ (τ), ω_i(τ), i=1,2,3,4, first Assume Fourier progression forms：

$\mathrm{\α}\left(\mathrm{\τ}\right)={\mathrm{\α}}_0+\underset{n=1}{\overset{N}{\Σ}}\[{\mathrm{\α}}_{2n-1}cos\left(n\mathrm{\ω}\mathrm{\τ}\right)+{\mathrm{\α}}_{2n}sin\left(n\mathrm{\ω}\mathrm{\τ}\right)\],$

$\mathrm{\ξ}\left(\mathrm{\τ}\right)={\mathrm{\ξ}}_0+\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[{\mathrm{\ξ}}_{2n-1}\cos \left(n\mathrm{\ω}\mathrm{\τ}\right)+{\mathrm{\ξ}}_{2n}\sin \left(n\mathrm{\ω}\mathrm{\τ}\right)\],$

$w_i\left(\mathrm{\τ}\right)=w_0^i+\underset{n=1}{\overset{N}{\Σ}}\[w_{2n-1}^{i}cos\left(n\mathrm{\ω}\mathrm{\τ}\right)+w_{2n}^{i}sin\left(n\mathrm{\ω}\mathrm{\τ}\right)\],i=1,2,3,4.---\left(5\right)$

Wherein, α_j, ξ_j, ω_j, j=0 ..., 2N is corresponding harmonic constant；

The approximate solution (5) that it will be assumed is substituted into differential equation group (4), obtains residual error function：

$R_1=c_0\stackrel{\·\·}{\mathrm{\ξ}}+c_1\stackrel{\·\·}{\mathrm{\α}}+c_2\stackrel{\·}{\mathrm{\ξ}}+c_3\stackrel{\·}{\mathrm{\α}}+(c_4+c_{10})\mathrm{\ξ}+c_5\mathrm{\α}+c_6{w}_1+c_7{w}_2+c_8{w}_3+c_9{w}_4\≠0,$

$R_2=d_0\stackrel{\·\·}{\mathrm{\ξ}}+d_1\stackrel{\·\·}{\mathrm{\α}}+d_2\stackrel{\·}{\mathrm{\ξ}}+d_3\stackrel{\·}{\mathrm{\α}}+d_4\mathrm{\ξ}+d_5\mathrm{\α}+d_6{w}_1+d_7{w}_2+d_8{w}_3+d_9{w}_4+d_{10}{\mathrm{\α}}^3\≠0,$

$R_3={\stackrel{\·}{w}}_1+{\∈}_1{w}_1-\mathrm{\α}\≠0,$

$R_4={\stackrel{\·}{w}}_2+{\∈}_2{w}_2-\mathrm{\α}\≠0,$

$R_5={\stackrel{\·}{w}}_3+{\∈}_1{w}_3-\mathrm{\ξ}\≠0,$

$R_6={\stackrel{\·}{w}}_4+{\∈}_2{w}_4-\mathrm{\ξ}\≠0.---\left(6\right)$

Wherein, R_jRepresent

Finally, residual error function R is forced_j2N+1 equidistant time point τ over one period_iOn be zero to obtain final product to time discrete method Algebraic Equation set, the aerofoil system contains 6 × (2N+1) individual equations ,+1 unknown number of 6 × (2N+1).

4. the rapid simulation method that a kind of re-entry space vehicle wing flutter according to claim 3 is responded, it is characterised in that Step 3) it is specially：

Time discrete treatment is carried out to time discrete method Algebraic Equation set each, to α (τ) in 2N+1 time point τ_iIt is discrete, can ：

$\mathrm{\α}\left({\mathrm{\τ}}_i\right)={\mathrm{\α}}_0+\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[{\mathrm{\α}}_{2n-1}cos\left({\mathrm{n\ω\τ}}_i\right)+{\mathrm{\α}}_{2n}sin\left({\mathrm{n\ω\τ}}_i\right)\].$

Wherein θ_i=ω τ_i, above formula is written as matrix form：

It is represented simply as：

${\tilde{Q}}_{\mathrm{\α}}\≡\left(\begin{array}{c}\mathrm{\α}\left({\mathrm{\τ}}_0\right)\\ \mathrm{\α}\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \mathrm{\α}\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)={\mathrm{FQ}}_{\mathrm{\α}},$

Wherein,

$F=\left(\begin{array}{cccccc}1& {\mathrm{cos\θ}}_0& {\mathrm{sin\θ}}_0& ...& \cos \left({\mathrm{N\θ}}_0\right)& \sin \left({\mathrm{N\θ}}_0\right)\\ 1& {\mathrm{cos\θ}}_1& {\mathrm{sin\θ}}_1& ...& \cos \left({\mathrm{N\θ}}_1\right)& \sin \left({\mathrm{N\θ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& {\mathrm{cos\θ}}_{2N}& {\mathrm{sin\θ}}_{2N}& ...& \cos \left({\mathrm{N\θ}}_{2N}\right)& \sin \left({\mathrm{N\θ}}_{2N}\right)\end{array}\right)$

Q_α=[α₀,α₁,...,α_2N]^T

SimilarWhereinQ_ξIt is respectively relevant with ξ with point Harmonic coefficient,WithRespectively and ω_iIt is relevant with a harmonic coefficient；

Then, it is rightDiscrete obtaining is carried out at 2N+1 time point：

$\stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_i\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[-{\mathrm{n\ω\α}}_{2n-1}sin\left({\mathrm{n\ω\τ}}_i\right)+{\mathrm{n\ω\α}}_{2n}cos\left({\mathrm{n\ω\τ}}_i\right)\]$

Above formula is written as matrix form：

Order matrix

Therefore have：

${\tilde{Q}}_{\stackrel{\·}{\mathrm{\α}}}={\mathrm{\ωFAQ}}_{\mathrm{\α}},---\left(8\right)$

Wherein

${\tilde{Q}}_{\stackrel{\·}{\mathrm{\α}}}\≡\left(\begin{array}{c}\stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_0\right)\\ \stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \stackrel{\·}{\mathrm{\α}}\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)$

Similarly have：WhereinWithIt is respectively On ξ and ω_jDerivative with point；

Then, it is rightIt is discrete at 2N+1 time point, have：

$\stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_i\right)=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\[-n^2{\mathrm{\ω}}^2{\mathrm{\α}}_{2n-1}cos\left({\mathrm{n\ω\τ}}_i\right)-n^2{\mathrm{\ω}}^2{\mathrm{\α}}_{2n}sin\left({\mathrm{n\ω\τ}}_i\right)\]---\left(9\right)$

Above formula (9) is written as matrix form：

Square formation FA in above formula²It is expressed as

${\mathrm{FA}}^2=\left(\begin{array}{cccccc}0& -{\mathrm{cos\θ}}_0& -{\mathrm{sin\θ}}_0& ...& -N^2\cos \left({\mathrm{N\θ}}_0\right)& -N^2\sin \left({\mathrm{N\θ}}_0\right)\\ 0& -{\mathrm{cos\θ}}_1& -{\mathrm{sin\θ}}_1& ...& -N^2\cos \left({\mathrm{N\θ}}_1\right)& -N^2\sin \left({\mathrm{N\θ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& -{\mathrm{cos\θ}}_{2N}& -{\mathrm{sin\θ}}_{2N}& ...& -N^2\cos \left({\mathrm{N\θ}}_{2N}\right)& -N^2\sin \left({\mathrm{N\θ}}_{2N}\right)\end{array}\right)$

Therefore,

${\tilde{Q}}_{\stackrel{\·\·}{\mathrm{\α}}}\≡\left(\begin{array}{c}\stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_0\right)\\ \stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ \stackrel{\·\·}{\mathrm{\α}}\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)={\mathrm{\ω}}^2{\mathrm{FA}}^2{Q}_{\mathrm{\α}}$

Similarly haveWhereinWithIt is respectively on ξ and ω_jSecond dervative with point；

Thus following form is obtained after time discrete method Algebraic Equation set is transformed：

$(c_0{\mathrm{\ω}}^2{D}^2+c_2\mathrm{\ω}D+c_{4}I+c_{10}I){\tilde{Q}}_{\mathrm{\ξ}}+(c_1{\mathrm{\ω}}^2{D}^2+c_3\mathrm{\ω}D+c_{5}I){\tilde{Q}}_{\mathrm{\α}}+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_{i+5}{\tilde{Q}}_{w_i}=0,$

$(d_0{\mathrm{\ω}}^2{D}^2+d_2\mathrm{\ω}D+d_{4}I){\tilde{Q}}_{\mathrm{\ξ}}+(d_1{\mathrm{\ω}}^2{D}^2+d_3\mathrm{\ω}D+d_{5}I){\tilde{Q}}_{\mathrm{\α}}+\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{d}_{i+5}{\tilde{Q}}_{w_i}+d_{10}{\tilde{M}}_{\mathrm{\α}}=0,$

$(\mathrm{\ω}D+{\∈}_{1}I){\tilde{Q}}_{{\mathrm{\ω}}_1}-{\tilde{Q}}_{\mathrm{\α}}=0,$

$(\mathrm{\ω}D+{\∈}_{2}I){\tilde{Q}}_{{\mathrm{\ω}}_2}-{\tilde{Q}}_{\mathrm{\α}}=0,$

$(\mathrm{\ω}D+{\∈}_{1}I){\tilde{Q}}_{{\mathrm{\ω}}_3}-{\tilde{Q}}_{\mathrm{\ξ}}=0,$

$(\mathrm{\ω}D+{\∈}_{2}I){\tilde{Q}}_{{\mathrm{\ω}}_4}-{\tilde{Q}}_{\mathrm{\ξ}}=0,---\left(10\right)$

Wherein D=FAF^-1, and

${\tilde{M}}_{\mathrm{\α}}=\left(\begin{array}{c}{\mathrm{\α}}^3\left({\mathrm{\τ}}_0\right)\\ {\mathrm{\α}}^3\left({\mathrm{\τ}}_1\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\α}}^3\left({\mathrm{\τ}}_{2N}\right)\end{array}\right)$

WillWithWithRepresent, then they are updated in nonlinear equation and are obtained：

$(A_{2\mathrm{\α}}-B_{2\mathrm{\ξ}}{B}_{1\mathrm{\ξ}}^{-1}{A}_{1\mathrm{\α}}){\tilde{Q}}_{\mathrm{\α}}+d_{10}{\tilde{M}}_{\mathrm{\α}}=0,---\left(11\right)$

Wherein A_1α, B_1ξ, A_2α, B_2ξRespectively：

A_1α=c₁ω²D²+c₃ωD+c₅I+c₆(ωD+∈₁I)^-1+c₇(ωD+∈₂I)^-1

B_1ξ=c₀ω²D²+c₂ωD+(c₄+c₁₀)I+c₈(ωD+∈₁I)^-1+c₉(ωD+∈₂I)^-1

A_2α=d₁ω²D²+d₃ωD+d₅I+d₆(ωD+∈₁I)^-1+d₇(ωD+∈₂I)^-1

B_2ξ=d₀ω²D²+d₂ωD+d₄I+d₈(ωD+∈₁I)^-1+d₉(ωD+∈₂I)^-1.

Formula (11) is the compact time discrete method Algebraic Equation set with time domain variable as variable, when using-frequency conversion relational expressionWithAbove formula is converted to the compact time discrete method generation with frequency variable as variable Number equation group (12)：

$(A_2-B_2{B}_1^{-1}{A}_1)Q_{\mathrm{\α}}+d_{10}{F}^{-1}{\left({\mathrm{FQ}}_{\mathrm{\α}}\right)}^3=0,---\left(12\right)$

Formula (12) containing 2N+2 unknown number, wherein, A₂=A_2α, A₁=A_1α, B₂=B_2ξ, B₁=B_1ξ。

5. the rapid simulation method that a kind of re-entry space vehicle wing flutter according to claim 1 is responded, it is characterised in that Step 4) the middle use scalar Homotopy above-mentioned compact time discrete method Algebraic Equation set of solution.