CN112487574B - Combustion stability margin assessment method - Google Patents
Combustion stability margin assessment method Download PDFInfo
- Publication number
- CN112487574B CN112487574B CN202011332465.9A CN202011332465A CN112487574B CN 112487574 B CN112487574 B CN 112487574B CN 202011332465 A CN202011332465 A CN 202011332465A CN 112487574 B CN112487574 B CN 112487574B
- Authority
- CN
- China
- Prior art keywords
- time sequence
- pressure
- combustion
- power spectrum
- combustion chamber
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/10—Noise analysis or noise optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/12—Timing analysis or timing optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Geometry (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Combined Controls Of Internal Combustion Engines (AREA)
- Measuring Fluid Pressure (AREA)
Abstract
The invention provides a combustion stability margin assessment method, which solves the problem that the prior method can not accurately acquire the natural mode frequency and the non-stability judgment criterion when a system stably works, and comprises the following steps: comprises the steps of collecting a combustion chamber pressure pulsation time sequence; calculating a pressure time sequence power spectrum; obtaining peak frequency of the power spectrum through parameter identification; band-pass filtering is carried out by taking peak frequency as a center; calculating an autocorrelation function of the filtered pressure time sequence; obtaining the envelope of the autocorrelation function through Hilbert transformation; fitting the envelope curve by using an exponential function to obtain an attenuation coefficient; the combustion process stability in the combustion chamber is evaluated based on the damping coefficient. The method can evaluate the margin of combustion stability in the combustion chamber in the combustion noise stage, and has guiding effect on evaluating the working reliability of the engine and evaluating the effectiveness of the improved scheme.
Description
Technical Field
The invention belongs to the field of aerospace propulsion systems, and particularly relates to a combustion stability margin assessment method.
Background
During rocket engine development, external pulse disturbances are typically introduced to evaluate the combustion stability of the engine. If the disturbance triggers unstable combustion, the system is in a bistable region, and the stability margin of the system is insufficient. If the oscillation amplitude generated by disturbance gradually decays with time, the distance between the system and the bifurcation point (theoretically, the attenuation coefficient of the system at the bifurcation point is zero) can be judged by calculating the attenuation coefficient, namely, the stability margin of the system. However, the introduction of external pulse disturbance to evaluate system stability requires additional excitation devices with the risk of triggering instability and thus damaging the engine, and thus methods to quantify system stability margin without external excitation devices are needed.
Methods available in the prior literature, such as Lieuwen, T, online Combustor Stability Margin Assessment Using Dynamic Pressure Data, ASME J.Eng.gas Turbines Power,2005,127 (3): 478-482.Lieuwen, based on nonlinear coupled vibrator models, calculate the system attenuation coefficient using the autocorrelation function of the steady stage pressure time series; such as Stadlmair, n.v., hummel, t., and Sattelmayer, t.thermoacousic damping rate determination from combustion noise using bayesian statistics, ASME Turbo expose 2017:Turbomachinery Technical Conference and Exposition (50848), v04at04a025.stadlmair is based on the Lieuwen method, and a bayesian statistical recognition method is adopted to simultaneously recognize attenuation coefficients of a plurality of acoustic modes; i.e., yi, t., gutmark, e.j., online prediction of the onset of combustion instabilities based on the computation of damping ratios.j.sound vib, 2008,310 (1-5): 442-447, yi converts the method of Lieuwen to the frequency domain and obtains the attenuation coefficient by fitting the pressure time series power spectrum.
Of the above methods, the Lieuwen method requires bandpass filtering at the center frequency to eliminate interactions between different acoustic modes. Because the pressure time sequence in the noise stage is broadband oscillation near the acoustic mode of the combustion chamber, the judgment of the existence of larger arbitrary center frequency can have larger influence on the calculation result. Lieuwen carries out band-pass filtering by taking the resonant frequency after the oscillation combustion of the system as the center frequency, and the method has two problems, namely, when the system parameters are changed (such as mixing ratio), the resonant frequency is changed; in practical applications, it is desirable to determine the stability margin of the system during the stable combustion (combustion noise) phase, and the resonant frequency of the system after the oscillating combustion cannot be known "accurately" at this time.
Theoretically, the system would be at zero decay factor, indicating that combustion instability of the system would occur, i.e., the threshold for decay factor is zero. In practical applications, combustion instability occurs in the combustion chamber when the damping coefficient does not reach zero, so that the threshold value of the damping coefficient needs to be corrected according to the existing experimental result.
Disclosure of Invention
The invention aims to solve the problem that the existing method cannot accurately acquire the natural mode frequency and the non-stability judgment criterion when a system works stably, and provides a combustion stability margin assessment method.
In order to achieve the above purpose, the technical scheme of the invention is as follows:
a combustion stability margin assessment method comprising the steps of:
step one, acquiring a combustion chamber pressure pulsation time sequence;
acquiring pressure pulsation time sequence of combustion chamber or combustion chamber before spraying by dynamic pressure sensor
Step two, calculating a pressure time sequence power spectrum;
calculating the power spectrum of the pressure time sequence by adopting a periodical graph method, namely, for the pressure pulsation time sequencePerforming discrete Fourier transform, dividing the square of the modulus by the number of data points to obtain a pressure time sequence power spectrum +.>
Step three, obtaining the peak frequency of the pressure time sequence power spectrum through parameter identification;
fitting a pressure time series power spectrum usingFitting by nonlinear least square methodObtaining peak frequency omega i ；
Wherein S is pp (ω；ω i ,Γ,ν i ) Power spectrum for theoretical pressure time series; omega is the angular frequency, Γ is the intensity of the noise, omega i Is the peak frequency, v i Is the attenuation coefficient;
step four, using peak frequency omega i Band-pass filtering is carried out for the center;
determining a filter bandwidth omega H At peak frequency omega i Band-pass filtering is carried out for the center, and the time sequence of the filtered pressure is recorded as p t ；
Step five, calculating an autocorrelation function rho of the filtered pressure time sequence τ ；
Autocorrelation function ρ τ The calculation is performed by the following method,
wherein E is calculated as the mathematical expectation, p t+τ Is p t Time series with a delay of τ, μ being p t Is a function of the average value or mathematical expectation of (c),is p t Is a variance of (2);
step six, obtaining an envelope curve H of the autocorrelation function through Hilbert transformation t ；
Obtaining an envelope curve H of an autocorrelation function by Hilbert transform t The calculation is as follows:
wherein t is the sampling time of the pressure pulsation time sequence, and τ is the delay time of the autocorrelation function;
step seven, fitting an envelope curve by using an exponential function to obtain an attenuation coefficient;
using exponential functionsFitting an envelope curve H t Obtaining attenuation coefficient v i ；
And step eight, comparing the distance between the attenuation coefficient and the threshold value, and evaluating the stability of the combustion process in the combustion chamber.
Further, in step five, the varianceThe following calculation is performed
Where n is the pressure time series length.
Further, in the eighth step, the threshold value is 0.016.
Compared with the prior art, the technical scheme of the invention has the following advantages:
the invention provides a combustion stability margin assessment method of a combustion chamber, which utilizes a system identification method to identify peak frequency from a pressure time sequence frequency domain power spectrum as the center frequency of band-pass filtering, reduces the randomness when the center frequency is selected in band-pass filtering calculation, and can calculate the attenuation coefficient of the natural mode frequency of the system more accurately; meanwhile, the invention establishes a conversion method of the attenuation coefficient and the empirical formula, obtains the attenuation coefficient threshold value when instability occurs, determines the stability judgment criterion, and improves the engineering applicability and reliability of the combustion stability margin assessment method.
Drawings
FIG. 1 is a schematic diagram of an oscillation curve of an autocorrelation function in an embodiment of the method of the present invention.
Detailed Description
The present invention will be described in further detail with reference to the accompanying drawings and specific embodiments.
Because the quantitative method for determining the center frequency of the band-pass filter is not available in the prior art, the invention provides a method for quantitatively determining the center frequency, namely a combustion stability margin assessment method is provided, and the margin of combustion stability in a combustion chamber can be assessed in a combustion noise stage by adopting the method, so that the method has a guiding effect on assessing the working reliability of an engine and assessing the effectiveness of an improvement scheme. The method comprises the following steps: comprises the steps of collecting a combustion chamber pressure pulsation time sequence; calculating a pressure time sequence power spectrum; obtaining peak frequency of the power spectrum through parameter identification; band-pass filtering is carried out by taking peak frequency as a center; calculating an autocorrelation function of the filtered pressure time sequence; obtaining the envelope of the autocorrelation function through Hilbert transformation; fitting the envelope curve by using an exponential function to obtain an attenuation coefficient; the combustion process stability in the combustion chamber is evaluated based on the damping coefficient.
The principle of the method is as follows: the random broadband signal (white noise) excited oscillation system (combustion chamber) can be described by the following equation,
wherein p is i Amplitude, omega of pressure oscillation of ith order mode i Is the peak frequency of the ith order mode, v i For the i-th order modal attenuation coefficient, xi (t) is a noise excitation term, and solving the power spectral density S of the upper pressure oscillation amplitude pp (ω) to give the following formula,
wherein ω is angular frequency, Γ is the intensity of noise ζ (t), the pressure time series obtained in the test satisfies formula (1), and the frequency domain power spectrum form satisfies formula (2), so the peak frequency can be obtained by fitting the above formula.
The combustion stability margin assessment method provided by the invention specifically comprises the following steps:
step one, acquiring a combustion chamber pressure pulsation time sequence;
monitoring of combustion chamber or combustion chamber pre-injection pressure pulse time sequence by dynamic pressure sensorAnd saving the measurement result;
step two, calculating a pressure time sequence power spectrum;
calculating the power spectrum of the pressure time sequence by adopting a classical periodogram method, namely, the pressure pulsation time sequencePerforming Discrete Fourier Transform (DFT), and dividing the square of modulus (amplitude function) by the number of data points to obtain pressure time sequence power spectrumIn addition, the pressure time series power spectrum can also be calculated by adopting an autocorrelation method or other parameter and non-parameter model spectrum estimation methods>
Step three, obtaining peak frequency omega of pressure time sequence power spectrum through parameter identification i ；
Fitting a pressure time series power spectrum usingFitting uses a nonlinear least square method, in particular by finding a set of undetermined parameters (ω i ,Γ,ν i ) The sum of squares of the residual errors is minimized, and the calculation can be carried out by adopting Gauss-Newton, levenberg-Marquardt algorithm and the like;
wherein S is pp (ω；ω i ,Γ,ν i ) Power spectrum for theoretical pressure time series; omega is the angular frequency, Γ is the intensity of the noise, omega i Is the peak frequency, v i Is the attenuation coefficient;
step four, band-pass filtering is carried out by taking peak frequency as a center;
determining a filter bandwidth omega H At peak frequency omega i Band-pass filtering is carried out for the center, and the time sequence of the filtered pressure is recorded as p t ；
The bandpass filtering algorithm is mature, taking a Butterworth filter as an example, and the general flow is as follows: firstly, designing a corresponding analog filter according to parameter requirements, then converting the analog filter into a digital filter by an impulse response invariant method or a bilinear conversion method, and recording a filtered pressure time sequence as p t ；
Step five, calculating an autocorrelation function of the filtered pressure time sequence;
autocorrelation function ρ τ The calculation is performed by the following method,
wherein E is calculated as the mathematical expectation, p t+τ Is p t Time series with a delay of τ, μ being p t Is a function of the average value or mathematical expectation of (c),is p t Variance, variance>The following calculation is performed
Wherein n is the length of the pressure time series;
step six, obtaining an envelope curve H of the autocorrelation function through Hilbert transformation t ；
Obtaining an envelope curve H of an autocorrelation function by Hilbert transform t The calculation is as follows:
wherein t is the sampling time of the pressure pulsation time sequence, and τ is the delay time of the autocorrelation function;
step seven, fitting an envelope curve by using an exponential function to obtain an attenuation coefficient;
using exponential functionsFitting an envelope curve H t Obtaining attenuation coefficient v i ；
In this step, it is equivalent to using the linear expression y= - ω i ν i t+b fitting curve lnH t A linear least square method is adopted, in particular to find a group of undetermined parameters (v i B), enabling the square sum of the residual errors to be minimum, wherein the linear least square method is an existing maturation algorithm;
step 8) comparing the attenuation coefficient with the distance of the threshold value to evaluate the stability of the combustion process in the combustion chamber.
The farther the decay factor is from the threshold value, the more stable the combustion state in the combustion chamber, and correspondingly, the closer the decay factor is from the threshold value, the more unstable the combustion state in the combustion chamber.
Document M.L.Dranovsky, V. (Ed.) Yang, F. (Ed.) Culick, and D. (Ed.) Talley.Combition Instabilities in Liquid Rocket Engines: testing and Development Practices in Russia.AIAA Progress in Astronautics and Aeronautics, volume 221,2007. Pulse pressure decay rate δT is defined, and a number of experiments have shown that the combustion chamber can operate stably within a range δT≡0.1 … 0.3.3. By calculation, the attenuation coefficient v i Is related with the attenuation rate delta T as v i =δt/2pi, so the attenuation coefficient threshold can be determined to be 0.016. The farther the decay factor is from the threshold value, the more stable the combustion state within the combustion chamber is indicated. Accordingly, the closer the decay factor is to the threshold value, the more unstable the combustion state in the combustion chamber.
The process according to the invention is described in further detail below by way of example.
Monitoring of combustion chamber or combustion chamber pre-injection pressure pulsation time sequence by dynamic pressure sensorIn the example, the sampling frequency is 12.8kHz, the sampling time is 5s, and the measurement result is saved;
for pressure pulsation time seriesPerforming Discrete Fourier Transform (DFT), dividing the square of modulus (amplitude function) by the number of data points n, and calculating pressure time series power spectrum +.>
Obtaining a power spectrum by fittingPeak frequency omega of (2) i In this example, the nonlinear least square method is used to perform parameter identification calculation, and in this example, the peak frequency ω i 2900Hz;
determining a filter bandwidth omega H At peak frequency omega i Band-pass filtering is carried out for the center, and the time sequence of the filtered pressure is recorded as p t Filter bandwidth extraction10% -20% of peak frequency, filtering bandwidth omega in this example H For peak frequency omega i 20% of (a), i.e. 580Hz;
calculating an autocorrelation function ρ of a filtered pressure time series τ The calculation formula is shown in step five, in order to effectively identify the envelope curve in the next step, the autocorrelation curve is calculated to have a length of 3-10 oscillation cycles (as shown in fig. 1, the calculation result of the step five autocorrelation function is an oscillation curve, but the length is controlled by a delay time τ, which is given in advance, and in order to effectively identify the envelope curve in the next step, the length of the autocorrelation function is controlled to have 3-10 oscillation cycles by giving a suitable value of the delay time τ), and the autocorrelation function ρ is calculated in this example τ About 8 oscillation cycles;
the envelope of the autocorrelation function is calculated by Hilbert transform, the calculation formula is shown in step six, and the exponential function is utilizedFitting an envelope curve H t Obtaining attenuation coefficient v i Attenuation Rate v obtained in this example i A damping coefficient threshold of 0.016 greater than 0.11 indicates a stable combustion state in the combustion chamber.
Claims (3)
1. A combustion stability margin evaluation method, characterized by comprising the steps of:
step one, acquiring a combustion chamber pressure pulsation time sequence;
acquiring pressure pulsation time sequence of combustion chamber or combustion chamber before spraying by dynamic pressure sensor
Step two, calculating a pressure time sequence power spectrum;
calculating the power spectrum of the pressure time sequence by adopting a periodical graph method, namely, for the pressure pulsation time sequencePerforming discrete Fourier transform, dividing the square of the modulus by the number of data points to obtain a pressure time sequence power spectrum +.>
Step three, obtaining the peak frequency of the pressure time sequence power spectrum through parameter identification;
fitting a pressure time series power spectrum usingFitting by using a nonlinear least square method to obtain peak frequency omega i ；
Wherein S is pp (ω；ω i ,Γ,ν i ) Power spectrum for theoretical pressure time series; omega is the angular frequency, Γ is the intensity of the noise, omega i Is the peak frequency, v i Is the attenuation coefficient;
step four, using peak frequency omega i Band-pass filtering is carried out for the center;
determining a filter bandwidth omega H At peak frequency omega i Band-pass filtering is carried out for the center, and the time sequence of the filtered pressure is recorded as p t ；
Step five, calculating an autocorrelation function rho of the filtered pressure time sequence τ ；
Autocorrelation function ρ τ The calculation is performed by the following method,
wherein E is calculated as the mathematical expectation, p t+τ Is p t Time series with a delay of τ, μ being p t Is a function of the average value or mathematical expectation of (c),is p t Is a variance of (2);
step six, obtaining an envelope curve H of the autocorrelation function through Hilbert transformation t ；
Obtaining an envelope curve H of an autocorrelation function by Hilbert transform t The calculation is as follows:
wherein t is the sampling time of the pressure pulsation time sequence, and τ is the delay time of the autocorrelation function;
step seven, fitting an envelope curve by using an exponential function to obtain an attenuation coefficient;
using exponential functionsFitting an envelope curve H t Obtaining attenuation coefficient v i ；
And step eight, comparing the distance between the attenuation coefficient and the threshold value, and evaluating the stability of the combustion process in the combustion chamber.
2. The combustion stability margin evaluation method according to claim 1, characterized in that: in step five, the varianceThe following calculation is performed
Where n is the pressure time series length.
3. The combustion stability margin evaluation method according to claim 1 or 2, characterized in that: in step eight, the threshold is 0.016.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011332465.9A CN112487574B (en) | 2020-11-24 | 2020-11-24 | Combustion stability margin assessment method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011332465.9A CN112487574B (en) | 2020-11-24 | 2020-11-24 | Combustion stability margin assessment method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112487574A CN112487574A (en) | 2021-03-12 |
| CN112487574B true CN112487574B (en) | 2023-08-04 |
Family
ID=74934251
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011332465.9A Active CN112487574B (en) | 2020-11-24 | 2020-11-24 | Combustion stability margin assessment method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112487574B (en) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106709144A (en) * | 2016-11-22 | 2017-05-24 | 中国人民解放军装备学院 | Autocorrelation theory-based engine instability prediction and assessment method |
| CN109184913A (en) * | 2018-10-08 | 2019-01-11 | 南京航空航天大学 | The aero-engine aerodynamic stability active composite control method with prediction is estimated based on stability |
| CN109344510A (en) * | 2018-10-08 | 2019-02-15 | 南京航空航天大学 | A kind of active stability control method based on the estimation of aero-engine stability margin |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7194382B2 (en) * | 2004-02-06 | 2007-03-20 | Georgia Tech Research Corporation | Systems and methods for detection of combustor stability margin |
-
2020
- 2020-11-24 CN CN202011332465.9A patent/CN112487574B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106709144A (en) * | 2016-11-22 | 2017-05-24 | 中国人民解放军装备学院 | Autocorrelation theory-based engine instability prediction and assessment method |
| CN109184913A (en) * | 2018-10-08 | 2019-01-11 | 南京航空航天大学 | The aero-engine aerodynamic stability active composite control method with prediction is estimated based on stability |
| CN109344510A (en) * | 2018-10-08 | 2019-02-15 | 南京航空航天大学 | A kind of active stability control method based on the estimation of aero-engine stability margin |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112487574A (en) | 2021-03-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102436515B (en) | For estimating the method and system of the residual life of turbine airfoil | |
| George et al. | Experimental comparison of axial turbine performance under steady and pulsating flows | |
| CN106197807B (en) | A kind of measurement method for dynamic force | |
| EP2446234B1 (en) | Determining the resonance parameters for mechanical oscillators | |
| CN106709144B (en) | Engine instability prediction and evaluation method based on autocorrelation theory | |
| CN109684898A (en) | Aero-engine and its vibration signal blind separating method and device | |
| Stadlmair et al. | Thermoacoustic damping rate determination from combustion noise using bayesian statistics | |
| US10830659B2 (en) | Method of controlling a positive-ignition internal combustion engine by means of a knock estimator | |
| CN112487574B (en) | Combustion stability margin assessment method | |
| Feng et al. | Experimental investigation on characteristics of cavitation-induced vibration on the runner of a bulb turbine | |
| Costall et al. | Fundamental characterization of turbocharger turbine unsteady flow behavior | |
| RU2542162C1 (en) | Method of diagnostics of pre-emergency modes of operation of dry rocket engines (dre) in hold down tests | |
| CN116306359A (en) | Random wave simulation method for multidimensional space non-stationary non-uniform wind field | |
| Torres et al. | Identification of the fuel-thrust dynamics of a gas turbo engine | |
| CN105973516B (en) | A method of the pulsation thrust of solid propellant rocket for identification | |
| CN114838924A (en) | Structural damping ratio identification method based on wind-induced vibration non-stationary response | |
| EP1147378A1 (en) | Analysis method and analyzer | |
| CN112487573A (en) | Online prediction method for combustion instability of combustion chamber | |
| Hernley et al. | Aerodynamic Forcing Models for Compressor Aeromechanics | |
| CN112149291B (en) | Weak harmonic signal detection system and method | |
| JP2002188955A (en) | Strength deterioration detection method of structure by using ambiguous external force | |
| Balasubramanian et al. | Parametric estimation of dynamical system data using autoregressive modeling | |
| CN114358072B (en) | Leakage flow structure vibration source identification method based on prototype actual measurement vibration response | |
| Fernelius | Determining the Pressure Pulse Severity of a Pressure Gain Combustion Device | |
| Tukmakov et al. | Definition of the Generation Instants in Time of the Running Shaft Resonant Flexural Vibrations by Means of Indicator Function of the Number of States of Dynamic System |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |