Method for identifying acoustic modal attenuation rate and uncertainty thereof under pulse excitation

Abstract

The invention relates to a method for identifying the acoustic modal attenuation rate and the uncertainty thereof under pulse excitation, belongs to the technical field of acoustic modal identification, and solves the technical problems that only the acoustic modal attenuation rate of a combustion chamber can be identified and the uncertainty thereof cannot be obtained under the action of pulse excitation, and the identification method comprises the following steps: establishing a regression equation of pressure data; determining a prior probability distribution for each parameter in the acoustic modal parameter vector; writing posterior probability distribution based on a Bayesian statistical inference theory; a Markov chain Monte Carlo method is adopted, posterior probability distribution is solved, and the acoustic modal attenuation rate and uncertainty are recognized. The identification method is used for rocket engine combustion stability margin test evaluation.

Description

Technical Field

[0001] The invention belongs to the technical field of acoustic mode recognition, and in particular relates to a method for identifying an acoustic mode attenuation rate and its uncertainty under pulse excitation. Background Art

[0002] Pulse excitation is a commonly used method in rocket engine combustion stability test evaluation, which mainly includes two methods: 1) installing explosive bombs at specific positions in the combustion chamber, and 2) installing pulse guns on the side walls of the combustion chamber. Foreign large-thrust engines have all undergone combustion stability evaluations, such as SSME and its upgraded version, RS-68, F1, Vulcain, LE-7, RD-0120, RD-120, and RD-170. In China, the pulse gun perturbation method (YF-X1~X5) is used in a series of conventional propellant engines. The pulse gun method is also used in solid rocket engine combustion response measurements (T-type burners) and full-scale engine stability tests.

[0003] Pulse excitation will produce strong pressure disturbances inside the combustion chamber, thereby exciting the pressure oscillations of the inherent acoustic mode of the combustion chamber. If the acoustic mode decay rate is small, the acoustic mode may show continuous large-amplitude pressure oscillations after excitation, that is, the combustion process and the acoustic mode undergo thermoacoustic coupling, which is called combustion instability. If the acoustic mode decay rate is large, the acoustic mode will show an oscillatory decay process after excitation and finally return to the noise level. For the latter, it is necessary to further identify the decay time or decay rate (characterizing the speed of oscillation decay) to evaluate the combustion stability margin of the combustion chamber.

[0004] The US CPIA 655 states that the main frequency component f of the combustion chamber is When the internal shock amplitude decays to below the allowable amplitude, the combustion stability margin of the engine is considered to be sufficient. The decay time threshold value obviously depends on the maximum shock amplitude and the allowable amplitude (which it eventually decays to). Russia requires that the maximum shock amplitude in the combustion chamber or generator be 15 to 25 times the root mean square value of the combustion noise, and the recommended decay time (the time taken for the maximum shock amplitude to decay to 1 / e times) should be less than 15ms. Some people take the time when the maximum shock amplitude decays to ±5% of the chamber pressure as the decay time. Some people have compared the decay time when the maximum shock amplitude decays to 10% and 1 / e times, and believe that the latter is more reasonable, and the former does not consider the influence of the maximum shock amplitude on the decay time.

[0005] The decay time is an overall evaluation of the decay speed of multiple acoustic modes, and cannot characterize the decay speed of each individual acoustic mode. When the decay time is slow, that is, the stability margin of the combustion chamber is insufficient and a passive damping device needs to be added, it is usually necessary to design for a single acoustic mode, so it is necessary to identify the decay time or decay rate of each individual acoustic mode.

[0006] Time-frequency analysis methods, such as short-time Fourier transform (STFT), wavelet transform (WT), Hilbert-Huang transform (HHT), digital filter + Hilbert transform (filter + HT), etc., are often used in the identification of structural vibration modal parameters. These methods can theoretically be applied to the identification of combustion chamber acoustic modes under pulse excitation, but the uncertainty of the identification results cannot be obtained. Summary of the invention

[0007] In order to solve the problem that the existing method can only identify the acoustic mode decay rate under pulse excitation but cannot obtain its uncertainty, the present invention proposes a method for identifying the combustion chamber acoustic mode decay rate and its uncertainty under pulse excitation.

[0008] The technical solution adopted by the present invention to solve the technical problem is:

[0009] A method for identifying the attenuation rate of the acoustic mode of a combustion chamber under pulse excitation and its uncertainty comprises the following steps:

[0010] Step S1, establishing a regression equation for pressure data;

[0011] Step S2, determine the prior probability distribution: for the acoustic modal parameter vector θ = (ν i ,ω i ,θ i ,A i ) determines the prior probability distribution.

[0012] Step S3, write out the posterior probability distribution: Based on the Bayesian statistical inference theory, write out the posterior probability distribution p(θy).

[0013] Step S4, identifying the acoustic mode decay rate and uncertainty: using the Markov Chain Monte Carlo (MCMC) method to solve the posterior probability distribution p(θy) and identify the acoustic mode decay rate and uncertainty.

[0014] In the above identification method, the step S1 of establishing a regression equation for pressure data further comprises:

[0015] The combustion chamber pressure data y(x, t) after pulse excitation contains several acoustic modal parameters, that is, the linear superposition of several acoustic modal parameters. The calculation formula of the combustion chamber pressure data y(x, t) is as follows:

[0016]

[0017] In formula (1), y(x,t) represents the pressure at position x at time t, η i (t) and ψ i (x) are the amplitude and vibration shape of the i-th order acoustic mode, and ξ(t) is the noise.

[0018] When the measuring point is fixed, the vibration mode function ψ i (x) degenerates to a constant B i , the pressure data y(x,t) is abbreviated as y(t), and the expression of the pressure data y(t) at a fixed measuring point is:

[0019]

[0020] When there is no nonlinear coupling between the acoustic modes, the amplitude of each acoustic mode is η in the free decay stage of the combustion chamber acoustic mode after pulse excitation. i (t) satisfies the following second-order linear vibration equation:

[0021]

[0022] In formula (3), ν i is the attenuation rate, ω i is the oscillation angular frequency.

[0023] Solving equation (3), we can obtain the amplitude of the i-th order acoustic mode of the combustion chamber as:

[0024]

[0025] In formula (4), C i is the initial amplitude, θ i For phase.

[0026] Combining equations (2) and (4), the pressure data y(t) at the fixed measuring point of the combustion chamber after pulse excitation is obtained. The regression equation of the pressure data y(t) at the fixed measuring point is:

[0027]

[0028] In formula (5), A i =B i C i .

[0029] In the above identification method, the step S2 of determining the prior probability distribution further comprises:

[0030] Assume that the acoustic modal parameter vector θ = (ν i ,ω i ,θ i ,A i ) where each component is independently distributed, the prior probability distribution p(θ) of the parameter vector θ is:

[0031] p(θ)=p(ν i )p(ω i )p(θ i )p(A i ) (7)

[0032] According to the characteristics of pulse excitation data and the regression equation of pressure data, the attenuation rate ν is determined i , oscillation angular frequency ω i , initial amplitude A i are all positive real values, and the phase θ i Take the value [-π,π]; further determine the attenuation rate ν i , oscillation angular frequency ω i , initial amplitude A i , phase θ i All satisfy uniform distribution, and ξ(t) satisfies normal distribution, that is:

[0033]

[0034] In formula (8), A max is the maximum initial amplitude, σ c is the standard deviation of the normal distribution, and the root mean square value of all data points of the pulsating pressure y(t) is taken; the acoustic mode number i is taken as the number of sudden frequencies after FFT analysis of the pressure data y(t) at a fixed measuring point, ω min and ω max are the minimum and maximum values ​​of the oscillation angular frequency respectively.

[0035] In the above identification method, the initial amplitude maximum value A max Take twice the maximum amplitude of the pressure data y(t) at the fixed measuring point and the minimum value of the oscillation angular frequency ω min and the maximum value ω max Take the pressure data y(t) at the fixed measuring point and perform FFT analysis on 0.7 times the minimum burst frequency and 1.2 times the maximum burst frequency.

[0036] In the above recognition method, in step S3, the expression of the posterior probability distribution p(θy) is as follows:

[0037] p(θy)∝p(yθ)p(θ)(6)

[0038] In formula (6), p(yθ) is the likelihood function of the pressure data y(t) at a fixed measuring point when the parameter vector θ is known.

[0039] In the above identification method, in step S4, sampling is performed by a random walk algorithm, and the sampling sequence {θ1, θ2, θ3, ...} constitutes a Markov chain, whose stationary distribution is the same as the probability distribution of the model parameters, and the mathematical expectation of the acoustic modal parameter vector is The calculation formula is as follows:

[0040]

[0041] The relationship between the confidence level and the posterior probability distribution p(θy) in the credible interval [δ1δ2] is as follows:

[0042]

[0043] In formula (10), the symbol β represents the confidence level, and the credible interval [δ1δ2] is the uncertainty of the mathematical expectation θ of the acoustic modal parameter vector.

[0044] In the above identification method, the random walk algorithm may be a Metropolis-Hastings algorithm or a Gibbs algorithm.

[0045] In the above identification method, the confidence level β=0.95.

[0046] The above identification method uses a dynamic pressure sensor to monitor the pressure data y(t) at a fixed measuring point in the combustion chamber after pulse excitation, wherein the sampling frequency is 10kHz and the sampling time is 15ms; A max Take twice the maximum value of the pressure data at the fixed measuring point, which is 40MPa; σ c The root mean square value of all data points of the pressure data y(t) at the fixed measuring point is 4.9MPa.

[0047] The above identification method performs fast Fourier transform FFT analysis on the pressure data y(t) at a fixed measuring point. The number of acoustic modes i is 2, ω min and ω max They are 0.70kHz and 2.64kHz respectively.

[0048] The beneficial effects of the present invention are:

[0049] A method for identifying the decay rate of acoustic modal parameters in a combustion chamber under pulse excitation and its uncertainty uses a Bayesian nonlinear regression method and empirical knowledge to establish the prior probability distribution and likelihood function (statistical model) of acoustic modal parameters. Based on experimental data and Bayesian theorem, the Markov chain Monte Carlo (MCMC) method is used to solve the posterior probability distribution. The marginal probability distribution of each acoustic modal parameter is obtained from the joint distribution, and its mathematical expectation is calculated as the estimated value of the parameter. The credible interval is calculated to quantify the uncertainty.

[0050] A method for identifying acoustic modal decay rate and its uncertainty under pulse excitation is proposed. The parameter estimation based on Bayesian statistical inference can not only identify the acoustic modal parameters, but also evaluate the uncertainty of the modal parameters. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 is the pulsating pressure time series of the second embodiment of the present invention;

[0052] Figure 2 This is the FFT analysis of the pulsating pressure time series in the second embodiment of the present invention;

[0053] Figure 3 is the probability density of the attenuation rate ν1 of the second embodiment of the present invention;

[0054] Figure 4 is the probability density of the attenuation rate ν2 in the second embodiment of the present invention. DETAILED DESCRIPTION

[0055] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0056] Embodiment 1

[0057] A method for identifying acoustic mode attenuation rate and its uncertainty under pulse excitation comprises the following steps:

[0058] (1) Establishing the regression equation of pulsating pressure data

[0059] The combustion chamber pressure data y(x,t) after pulse excitation is non-stationary data, containing multiple acoustic modal parameters, which can be written as the linear superposition of each acoustic modal parameter.

[0060]

[0061] Where y(x,t) represents the pressure time series at position x, η i (t) and ψ i (x) are the amplitude and vibration mode of the ith order acoustic mode, and ξ(t) is the noise. When the measuring point is fixed, the vibration mode function ψ i (x) degenerates to a constant B i , so y(x,t) can be shortened to y(t)

[0062]

[0063] Assuming that there is no nonlinear coupling between the acoustic modes, the amplitude of each acoustic mode in the combustion chamber during the free decay stage after pulse excitation is η i (t) satisfies the following second-order linear vibration equation:

[0064]

[0065] The general solution of the equation is

[0066]

[0067] Among them C i is the initial amplitude, ν i is the attenuation rate, ω i is the oscillation angular frequency, θ iCombining equations (2) and (4), the pulsating pressure data at a certain measuring point in the combustion chamber after pulse excitation is formed by the following equation:

[0068]

[0069] Among them A i =B i C i The purpose of this application is to use experimental data y(t) to identify the acoustic modal parameter vector θ = (ν i ,ω i ,θ i ,A i ) and quantify its uncertainty.

[0070] (2) Bayesian statistical model

[0071] Based on the Bayesian statistical inference theory, the acoustic modal parameter vector θ = (ν i ,ω i ,θ i ,A i ) can be written as conditional probability form p(θy). From Bayes’ theorem, we know

[0072] p(θy)∝p(yθ)p(θ)(6)where p(yθ) is the likelihood function of the data y(t) when the parameter vector θ is known, and p(θ) is the prior probability distribution of the parameter vector θ.

[0073] (3) Determine the prior probability distribution

[0074] The acoustic modal parameter vector θ = (ν i ,ω i ,θ i ,A i ) determines the prior probability distribution. Assume that the parameter vector θ = (ν i ,ω i ,θ i ,A i ) are independently distributed, that is,

[0075] p(θ)=p(ν i )p(ω i )p(θ i )p(A i ) (7)

[0076] According to the characteristics of pulse excitation data and formula (5), the attenuation rate ν can be determined i , oscillation angular frequency ω i , initial amplitude A i are all positive real values, and the phase θ i Take the value [-π,π]. Further determine the attenuation rate νi , oscillation angular frequency ω i , initial amplitude A i , phase θ i All satisfy uniform distribution, and ξ(t) satisfies normal distribution, that is,

[0077]

[0078] A max Take twice the maximum amplitude of the pulsating pressure data y(t), σ c is the standard deviation of the normal distribution, and the root mean square value of all data points of the pulsating pressure y(t) is taken. The acoustic mode number i is taken as the number of sudden frequencies after FFT analysis of the data y(t), ω min and ω max They are taken as 0.7 times the minimum burst frequency and 1.2 times the maximum burst frequency after FFT analysis of the data y(t).

[0079] (4) Calculate the posterior probability distribution

[0080] The Markov chain Monte Carlo (MCMC) method is used to solve the above posterior probability distribution p(θy). Specifically, sampling is performed through a random walk algorithm, such as the mature Metropolis-Hastings (MH) algorithm or the Gibbs algorithm. The sampling sequence {θ1, θ2, θ3, ...} constitutes a Markov chain, whose stationary distribution is the same as the probability distribution of the model parameters. The mathematical expectation is calculated using formula (9): As an estimate of the parameter θ

[0081]

[0082] Using formula (10), we can calculate the confidence interval [δ1δ2] with confidence level β for the parameter value quantify the uncertainty.

[0083] The confidence level is also called confidence degree or probability, with the symbol β and is usually 0.95.

[0084]

[0085] Embodiment 2

[0086] The dynamic pressure sensor is used to monitor the time series y(t) of the pulsating pressure in the combustion chamber after pulse excitation. In this example, the sampling frequency is 10 kHz and the sampling time is 15 ms. Figure 1 As shown. max Take twice the maximum amplitude of the pulsating pressure data, that is, 40MPa. c Take the root mean square value of all data points of the pulsating pressure y(t), which is 4.9MPa.

[0087] The fast Fourier transform (FFT) analysis of the pulsating pressure time series y(t) is performed, and the results are as follows: Figure 2 As shown. Determine the number of acoustic modes i to be 2, ω min and ω max They are 0.70kHz and 2.64kHz respectively.

[0088] The Markov chain Monte Carlo (MCMC) method is used to solve the posterior probability distribution and the probability density of the decay rates ν1 and ν2 is obtained as follows: Figure 3 , Figure 4 The mathematical expectation is calculated to obtain the estimated values ​​of the decay rates ν1 and ν2, which are 1.65 and 3.33 respectively, and the credible intervals with a probability of 0.95 are [1.55, 1.74] and [3.16, 3.50] respectively.

Claims

1. A method for identifying the decay rate of combustion chamber acoustic modes and its uncertainty under pulse excitation, characterized in that: The following steps are involved: Step S1, establishing a regression equation for pressure data; Step S2, determine the prior probability distribution: for the acoustic modal parameter vector θ=(ν i ,ω i ,θ i ,A i ) determines the prior probability distribution for each parameter; Step S3, write the posterior probability distribution: Based on the Bayesian statistical inference theory, write the posterior probability distribution p(θy); Step S4, identifying the acoustic mode decay rate and uncertainty: using the Markov Chain Monte Carlo (MCMC) method, Solve the posterior probability distribution p(θy) and identify the acoustic mode decay rate and uncertainty.

2. The identification method according to claim 1, characterized in that: The step S1 of establishing a regression equation for pressure data further includes: The combustion chamber pressure data y(x, t) after pulse excitation contains several orders of acoustic modes, that is, the linear superposition of several orders of acoustic modes. The calculation formula of the combustion chamber pressure data y(x, t) is as follows: In formula (1), y(x,t) represents the pressure at position x at time t, η i (t) and ψ i (x) are the amplitude and vibration shape of the i-th order acoustic mode, ξ(t) is the noise; When the measuring point is fixed, the vibration mode function ψ i (x) degenerates into a constant B i , the pressure data y(x,t) is abbreviated as y(t), and the expression of the pressure data y(t) at a fixed measuring point is: When there is no nonlinear coupling between the acoustic modes of each order, the acoustic mode of the combustion chamber is in the free decay stage after the pulse excitation. Amplitude of each acoustic mode η i (t) satisfies the following second-order linear vibration equation: In formula (3), ν i is the decay rate, ω i is the oscillation angular frequency; Solving equation (3), we can obtain the amplitude of the i-th order acoustic mode of the combustion chamber as follows: In formula (4), C i is the initial amplitude, θ i is the phase; Combining equations (2) and (4), the pressure data y(t) at the fixed measuring point in the combustion chamber after pulse excitation is obtained. The regression equation of the pressure data y(t) at the fixed measuring point is: In formula (5), A i =B i C i .

3. The identification method according to claim 1, characterized in that: The step S2 of determining the prior probability distribution further comprises: Assume that the acoustic modal parameter vector θ = (ν i ,ω i ,θ i ,A i ) where each component is independently distributed, and the prior probability distribution p(θ) of the parameter vector θ is: p(θ)=p(ν i )p(ω i )p(θ i )p(A i ) (7) According to the characteristics of pulse excitation data and the regression equation of pressure data, the attenuation rate ν is determined i , oscillation angular frequency ω i 、 Initial amplitude A i are all positive real values, and the phase θ i Take the value [-π,π]; further determine the attenuation rate ν i , oscillation angular frequency ω i , initial amplitude A i , phase θ i All satisfy uniform distribution, and ξ(t) satisfies normal distribution, that is: In formula (8), A max is the maximum initial amplitude, σ c is the standard deviation of the normal distribution, and the root mean square value of all data points of the pulsating pressure y(t) is taken; the acoustic mode number i is taken as the number of sudden frequencies after FFT analysis of the pressure data y(t) at a fixed measuring point, ω min and ω max are the minimum and maximum values ​​of the oscillation angular frequency respectively.

4. The identification method according to claim 3, characterized in that: The initial amplitude maximum value A max Take 2 times the maximum amplitude of the pressure data y(t) at the fixed measuring point and the minimum value of the oscillation angular frequency ω min and the maximum value ω max The pressure data y(t) at the fixed measuring point is taken as 0.7 times of the minimum burst frequency and 1.2 times of the maximum burst frequency after FFT analysis.

5. The identification method according to claim 1, characterized in that: In step S3, the posterior probability distribution p(θy) is expressed as follows: p(θy)∝p(yθ)p(θ)(6)In formula (6), p(yθ) is the likelihood function of the pressure data y(t) at a fixed measuring point when the parameter vector θ is known.

6. The identification method according to claim 1, characterized in that: In step S4, sampling is performed by a random walk algorithm, and the sampling sequence {θ1, θ2, θ3, ...} constitutes a Markov chain, whose stationary distribution is the same as the probability distribution of the model parameters, and the mathematical expectation of the acoustic modal parameter vector is The calculation formula is as follows: The calculation relationship of the credible interval [δ1δ2] with a confidence level of β for the posterior probability distribution p(θy) is as follows: In formula (10), the credible interval [δ1δ2] is the mathematical expectation of the acoustic modal parameter vector uncertainty.

7. The identification method according to claim 6, characterized in that: The random walk algorithm may be a Metropolis-Hastings algorithm or a Gibbs algorithm.

8. The identification method according to claim 6, characterized in that: The confidence level β=0.

95.

9. The identification method according to claim 1, characterized in that: The pressure data y(t) at a fixed measuring point in the combustion chamber after pulse excitation is monitored by a dynamic pressure sensor, where the sampling frequency is 10 kHz and the sampling time is 15 ms; A max Take twice the maximum value of the pressure data at the fixed measuring point, which is 40MPa; c The root mean square value of all data points of the pressure data y(t) at a fixed measuring point is 4.9 MPa.

10. The identification method according to claim 9, characterized in that: Perform fast Fourier transform (FFT) analysis on the pressure data y(t) at a fixed measuring point. The number of acoustic modes i is 2, ω min and ω max They are 0.70kHz and 2.64kHz respectively.

Images