# A transfer function relating jet engine thrust to pitch angle during a phugoid motion

Crossposted on Mathematics SE

I am doing a school project that requires me to find a transfer function that relates the thrust of a jet engine (which could change with time in one way or another) and the aircraft's pitch angle during a phugoid motion. Phugoid motion is when an aircraft pitches and flies up along this pitch angle while maintaining constant thrust. The drag and the weight would slow the aircraft down, causing the pitch angle to decrease as lift--which depends on aircraft velocity along the flight path. Due to the lift being no longer able to counter the weight, the aircraft droops down, accelerated by gravity and its own thrust while speeding up. This increases lift and allows it to climb back up again. Then the cycle is repeated because all of the flight controls are down. My project's goal is to use fuzzy logic to increase and decrease the thrust so the pitch angle and rate of pitch angle change would be damped to 0.

A set of equations describing the phugoid motion: https://www.math.stonybrook.edu/~scott/Book331/Seeing_flight_path.html

The $$\theta$$ is the pitch angle: https://www.math.stonybrook.edu/~scott/Book331/Phugoid_model.html

The way I see it, the acceleration/deceleration of the aircraft along the flight path can be found through: $$a_{fp}(t)=T(t)-D(t)-mg\sin(\theta(t))$$ m is the mass of the aircraft,g is the gravity, T is thrust, and D is drag.

The velocity along the flight path can be found through $$v_{fp}=v_{fp_0}-\int_{t_0}^{t_1}a_{fp}(t)dt$$ $$v_{fp_0}$$ is the initial aircraft speed along flight path (a given value). The thrust is a time-varied input, the drag can be found through $$D(t)=C_d\frac{\rho(v_{fp}(t))^2}{2}S_{wet}$$ Everything except vfp is known. The lift of the aircraft is just like the drag: $$L(t)=C_L\frac{\rho(v_{fp}(t))^2}{2}S_{wet}$$ The lift and drag equation here are all with respect to the inertia reference frame instead of with respect to the flight path. The pitch angle should be $$\theta(t)=\theta_0-\int_{t_0}^{t_1}\frac{\frac{L(t)-W}{m}}{\int_{t_0}^{t_1}v_{fp}(t)dt}dt$$ Where $$\frac{L(t)-W}{m}$$ is the vertical acceleration of the aircraft, $$\int_{t_0}^{t_1}v_{fp}(t)dt$$ is the displacement of the aircraft from $$t_0$$ to $$t_1$$ So this $$\frac{\frac{L(t)-W}{m}}{\int_{t_0}^{t_1}v_{fp}(t)dt}$$ is actually an angular acceleration. So we integrate this angular acceleration and find the change in pitch angle from initial $$t_0$$ to "final" $$t_1$$.

But I can't figure out what to do next, given that there are so many integrals.

• The lift of the aircraft is just like the drag - no, I don't think so. There should be a phase shift between the two. Next, for thrust variation with speed it is important to know the bypass ratio. Did you try to model the airplane using finite differences? That avoids those pesky integrals ... Commented Nov 23, 2023 at 16:37
• @PeterKämpf It is a school project, so I am going to simplify it and make thrust an input value of my Simulink system. I mean, it may be good if I can get a theta(t+1)=f(theta(t)) Commented Nov 23, 2023 at 18:26
• To perform such analysis the steps normally are: 1) write the equation of motion 2) linearize it via the small-perturbations hypothesis 3) do some simplifications according to the relative magnitude of each term in the equations 4) apply Laplace to change in the s-domain 5) from the transfer function get the searched values. Didn't your professor teach that? If not, any book about flight mechanics can guide you. Commented Nov 23, 2023 at 19:14
• Why not hold vfp constant, since that's your goal. Thrust damping could be based on Drag (now more or less constant too), so your fuzzy is only looking for +/- changes in velocity due to the gravity component of climb and descent angle. It would then add or subtract engine thrust accordingly to help damp the phugoid. Commented Nov 24, 2023 at 10:29
• For these kinds of control analysis and control synthesis, one typically uses a set of ordinary differential equations (ODE's). Phugoid motion (and Short term oscillation) are textbook examples for control of aircraft as these show the oscillation frequency and damping of these modes very well. I will seek out the according equations as soon as I have time. From that point on you can then easily design an according control system. Commented Nov 24, 2023 at 10:57

For control design, it is much more convenient to characterize the motion of a dynamic system via a set of Ordinary Differential Equations (ODE's) instead of directly trying to solve the equations... I will try to outline the "classical textbook" way of deriving and characterizing the phugoid motion of an aircraft which then can be used to design a control system. You can find similar derivations in most textbooks on aircraft dynamics and control, a free example is this one (although this might not be the best textbook on the matter).

For the analysis of the dynamic responses of aircraft, we first start with the well known (longitudinal) equations of motion of an aircraft. These are then linearized (a step which I am not showing). The resulting state space system in well known linear matrix equation form $$\dot{x}=\mathbf{A}\cdot x + \mathbf{B} \cdot u$$ is given as follows:

$$\begin{bmatrix} \dot{u}\\ \dot{w} \\ \dot{q} \\ \dot{\theta} \end{bmatrix} = \begin{bmatrix} X_u & X_w & 0 &-g \cdot \cos(\theta_0) \\ Z_u & Z_w & u_0 & -g \cdot \sin(\theta_0) \\ M_u+M_{\dot{w}} Z_u & M_w+M_{\dot{w}} Z_w & M_q+u_0 M_{\dot{w}} & M_{\dot{w}} \\ 0 & 0 & 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} \dot{u}\\ \dot{w} \\ \dot{q} \\ \dot{\theta} \end{bmatrix} + \begin{bmatrix} X_{\delta_e} & X_{\delta_T} \\ Z_{\delta_e} & Z_{\delta_T} \\ M_{\delta_e}+M_{\dot{w}} Z_{\delta_e} & M_{\delta_T}+M_{\dot{w}} Z_{\delta_T} \\ 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} \delta_e & \delta_T \end{bmatrix}$$ The nomenclature is as follows: $$u$$ and $$w$$ are the forward and downward velocities, $$q$$ the rotational velocities and $$\theta$$ the pitch angle. Thrust and elevator deflection are given by $$\delta_T$$ and $$\delta_e$$. The derivatives like $$X_u$$ and $$M_w$$ are the aerodynamic/thrust terms derivated in respect to the different variables (therefore $$u$$ and $$w$$ in this example). The source (p. 83) also states how to relate these to more commonly known aircraft parameters here.

This fourth order system now mirrors two distinct aircraft motions:

• The short period oscillation, therefore involving mainly $$w$$ and $$q$$
• The phugoid oscillation, involving mainly $$u$$ and $$\theta$$

As you specified, you want to control the phugoid motion only with engine thrust (I pressume in reference to flight incidents like this one). Therefore you can neglect $$\delta_e$$ and also the coefficients $$X_{\delta_e}$$, $$Z_{\delta_e}$$, and $$M_{\delta_e}$$. To calculate the transfer function of the fourth order system, you can now either directly form the transfer function by calculating $$H(s)=\mathbf{C}(s\mathbf{I} - \mathbf{A})^{-1}\mathbf{B} + \mathbf{D}$$ , or we split the system according to the two different modes into two systems. Once into a system reflecting the short period oscillation

$$\begin{bmatrix} \dot{w} \\ \dot{q} \end{bmatrix} = \begin{bmatrix} Z_w & u_0 \\ M_w+M_{\dot{w}} Z_w & M_q+u_0 M_{\dot{w}} \end{bmatrix} \begin{bmatrix} w \\ q \end{bmatrix} + \begin{bmatrix} Z_{\delta_T} \\ M_{\delta_T}+M_{\dot{w}} Z_{\delta_T} \end{bmatrix} \cdot \begin{bmatrix} \delta_T \end{bmatrix}$$

and once into a system reflecting the phugoid motion (after some reformulation and setting $$w=\dot{w}=0$$, p 88)

$$\begin{bmatrix} \dot{u}\\ \dot{\theta} \end{bmatrix} = \begin{bmatrix} X_u &-g \cdot \cos(\theta_0) \\ -\frac{Z_u}{u_0 + Z_q} & \frac{g \cdot \sin{\theta_0}}{u_0 + Z_q} \end{bmatrix} \cdot \begin{bmatrix} u \\ \theta \end{bmatrix} + \begin{bmatrix} X_{\delta_T} \\ 0 \end{bmatrix} \cdot \begin{bmatrix} \delta_T \end{bmatrix}$$

by neglecting all cross-coupling terms. It should be noted that by neglecting the cross-coupling terms, the accuracy of the system also diminishes. Therefore, if you can, you should avoid this, although it is easier to handle 2x2 systems directly.

Either way, I would strongly reccomend doing all of this numerically with tools such as matlab, which is very well suited for analyzing such systems, and then to also synthesize a controller.

Footnote: A much more convenient way of choosing the state variables is to choose the following state vector $$x=\begin{bmatrix} V & \alpha & q & \theta \end{bmatrix}^T$$, however sources of this choice are unfortunately closed source.