# Research on Real-Time Model of Turboshaft Engine with Surge Process
## 3 Real-Time Surge Model of Turboshaft Engine
### Stability Analysis of MG Model  
For a specific MG pure surge model of a compression system, there is a \(B_{cr}\) which represents the critical value of \(B\) and when the \(B\) parameter is less than \(B_{cr}\), the instability type of the MG model is the surge. Additionally, the value of \(B_{cr}\) can be obtained with Hopf bifurcation analysis [22; 23; 24]. Express Equations (14) and (15) in the form of \(\dot{x}=f(x,B)\) where \(x=\left[ \begin{array}{cc}\Phi&\Psi\end{array} \right]^{T}\) and suppose the equilibrium point of the system is \((\Phi_{0},\Psi_{0})\), then the Jacobian matrix is obtained:  
\[L=D_{x}f(\Phi_{0},\Psi_{0})=\left[ \begin{array}{cc}\frac{\partial f_{1}}{ \partial\Phi_{0}^{2}}&\frac{\partial f_{1}}{\partial\Psi_{0}^{2}}\\ \frac{\partial f_{2}}{\partial\Phi}&\frac{\partial\Psi}{\partial\Psi}\end{array}  \right]_{\Phi=\Phi_{0},\Psi=\Psi_{0}}=\left[ \begin{array}{cc}\frac{\dot{ \Psi}_{C}(\Phi_{0})}{l_{C}}&\frac{1}{l_{C}}\\ \frac{1}{4B^{2}l_{c}}&\frac{-\partial_{T}(\Psi_{0})}{4B^{2}l_{c}}\end{array} \right] \tag{16}\]  
The characteristic equation is obtained according to \(\det(\lambda I-L)=0\):  
\[\lambda^{2}+\frac{\left[l_{c}\dot{\Phi}_{T}(\Psi_{0})-4B^{2}l_{c}\dot{\Psi}c( \Phi_{0})\right]}{4B^{2}l_{c}^{2}}\lambda+\frac{1-\dot{\Psi}c(\Phi_{0})\dot{ \Phi}_{T}(\Psi_{0})}{4B^{2}l_{c}^{2}}=0 \tag{17}\]  
where \(\dot{\Phi}_{T}(\Psi_{0})\) and \(\dot{\Psi}c(\Phi_{0})\) are respectively the characteristic derivatives of throttle valve and compressor at the equilibrium point:  
\[\dot{\Phi}_{T}(\Psi_{0})=\frac{\gamma_{T}}{2\sqrt{\Psi_{0}}} \tag{18}\]  
\[\dot{\Psi}c(\Phi_{0})=\frac{3H}{2W}\Bigg{(}1-\bigg{(}\frac{\Phi_{0}}{W}-1 \bigg{)}^{2}\Bigg{)} \tag{19}\]  
The change of the real part of the eigenvalue of the matrix determines the number of special structures in the phase space and the properties of each of them. Considering MG model with \(H=0.18\), \(W=0.25\), \(\gamma_{T}=0.55\), the compressor characteristic parameters \(H\), \(W\) and \(\gamma_{T}\) keep unchanged to analyze the distribution of eigenvalue of the Jacobian matrix with different \(B\) parameters. The changing trajectory of the two eigenvalues with the parameter is shown in Figure 4. With the increase of the \(B\) parameter, the distribution of the two eigenvalues goes through three stages: two negative real eigenvalues, a pair of conjugate complex eigenvalues and two positive real eigenvalues. In the second stage, the pair of common complex eigenvalues cross the imaginary axis from the left half complex plane to the right half complex plane. When crossing the imaginary axis, the real part of the common eigenvalue becomes zero, which means the stability of the equilibrium point changes and a limit cycle is generated at the equilibrium point. When \(B=0.4404\), the real part of the two eigenvalues is 0 and Hopf bifurcation occurs. Therefore, the \(B_{cr}\) of MG model, in this case, is 0.4404. The classical MG model has certain compressor characteristics, that is, \(H\) and \(W\) are determined, so the throttle opening directly affects the \(B_{cr}\). For the characteristic equation as Equation (17), \(B_{cr}\) can be accurately solved by using the Routh Hurwitz criterion [25]. Figure 5 shows the value of \(B_{cr}\) under \(\gamma_{T}\) and they are positively correlated. For a given \(B\) parameter, the \(B_{cr}\) can be reduced by reducing \(\gamma_{T}\), so that the \(B\) parameter that originally made the system stable becomes greater than \(B_{cr}\), making the model enter or exit surge. Figure 4: The Variation of eigenvalue trajectories with \(B\) parameter: (**a**) Variation of two eigenvalues; (**b**) Variation of real part of two eigenvalues; (**c**) Variation of imaginary part of two eigenvalues. Use Matcont toolbox [26] to further analyze the change of limit cycle with \(B\) parameter during the surge, and the results are shown in Figure 6. When the \(B\) parameter is greater than the \(B_{cr}\), the limit cycle occurs and its size is related to the size of the \(B\) parameter.