# Obtaining stability derivatives from forced oscillations Consider a simple wing-tail design. The design has a wing span $$b\_{ref}=1.7321\text{ }m$$, wing area $$S\_{ref}=0.3\text{ }m^2$$, and wing chord of $$c\_{ref}=0.1732\text{ }m$$. The image below shows the mesh of the design. ![Mesh of a wing-tail design](https://2298927434-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-MWUM1iuDAQw73diDpYu%2F-MfJs7x1PRknD27-Jrg6%2F-MfJtASExNbeyRXfCyRw%2Fwing_tail_configuration.png?alt=media\&token=d1f53b45-6bca-4b27-b9ab-fb97b76da82a) A body-fixed coordinate system is placed at the leading edge of the wing. The reference speed is $$|V\_{ref} |=25\text{ }m/s$$ and the reference density, pressure and viscosity are evaluated at an altitude of $$h=0\text{ }m$$(ISA). Before performing forced oscillation simulations the location of the neutral point of the wing-tail configuration must be determined. Two simulations at different angles of attack are required to determine the neutral point location. Please refer to the [Stability article](https://blog.aviumtechnologies.com/topics/stability) for more information on how to determine the neutral point. The table below shows the $$C\_{Z}$$and $$C\_{m}$$values used for determining the location of the neutral point, $$x\_{np}$$. | Simulation type | $$\alpha^\circ$$ | $$C\_{Z}$$ | $$C\_{m}$$ | | --------------- | ---------------- | ---------- | ---------- | | Fixed-wake | 0 | −0.3215 | -0.0223 | | Fixed-wake | 5 | -0.7913 | -0.1973 | The location of the neutral point, $$x\_{np}$$, can be obtained from the following equation: $$ x\_{np}=-c\_{ref}\frac{\frac{\partial C\_{m}}{\partial \alpha}}{\frac{\partial C\_{Z}}{\partial \alpha}}\approx\frac{\Delta C\_{m}}{\Delta C\_{Z}}\approx\frac{C\_{m}\rvert\_{\alpha=5}-C\_{m}\rvert\_{\alpha=0}}{C\_{Z}\rvert\_{\alpha=5}-C\_{Z}\rvert\_{\alpha=0}}\text{.} $$ After substituting the $$C\_{Z}$$and $$C\_{m}$$ values from the table into the above equation, $$x\_{np}$$ evaluates to $$-0.0645\text{ }m$$. The minus sign indicates that the neutral point is aft of the origin of the body-fixed coordinate system. The following equation can be used to determine the centre of gravity for a specific static margin, $$SM$$: $$ SM=\frac{x\_{cg}-x\_{np}}{c\_{ref}} $$ For $$SM=0.1\text{ }(10%)$$, $$x\_{cg}=-0.0385\text{ }m$$. {% hint style="info" %} For a statically stable design the centre of gravity is ahead of the neutral point. {% endhint %} Once the centre of gravity location is obtained, forced oscillation simulations can be performed. Remember to change the x\_cg value in the **.conf** file and to re-run the preprocessor before proceeding. The preprocessor will move the origin of the body-fixed coordinate system to the new location specified by x\_cg, y\_cg, and z\_cg. Please refer to the [Configuration file section](https://docs.aviumtechnologies.com/user-guide/the-configuration-file) of the documentation for more information. ### Creating a forced oscillation motion file A motion file describing the translation and rotation of the body-fixed coordinate system is required to perform a forced oscillation simulation. Please refer to the [Custom motion section](https://docs.aviumtechnologies.com/user-guide/solver/unsteady-simulations/custom-motion) of the documentation for additional information. To obtain the longitudinal stability derivatives of the design, a forced pitch oscillation motion is required. The pitch angle during the oscillation is defined by the following equation: $$ \theta=Asin(\omega t), $$ where $$\theta$$ is the pitch angle, $$A$$is the amplitude of the pitch oscillation, $$\omega$$is the angular frequency of the oscillation ($$\omega=2\pi f$$), and $$t$$is the time. In this example $$A=5\text{ }deg$$and $$f=5 Hz$$. The pitch rate $$q$$can be obtained by differentiating the $$\theta$$equation with respect to time: $$ q=\dot{\theta}=\omega A cos(\omega t). $$ Appropriate time step size and number of time steps must be selected to complete at least 2 cycles. In this example $$dt=4\times10^{-3}$$ and $$N\_{timesteps}=100$$. The image below shows $$\theta$$and $$q$$. Note that $$q$$is divided by 100 so the difference between $$\theta$$ and $$q$$is easy to see. In the actual motion file $$q$$is not divided by 100. ![Generated pitch angle and pitch rate](https://2298927434-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-MWUM1iuDAQw73diDpYu%2F-MfKI0XI_pwaBk9mOKsm%2F-MfKJ_8SXasQNriEGjx5%2Fmotion.png?alt=media\&token=ee23e223-5243-4d8f-b3da-76d9ec22d5d4) ### Obtaining the stability derivatives from the simulation results The video below shows the development of the wing-tail configuration wake during the forced oscillation simulation. {% embed url="" %} Wing-tail design forced pitch oscillation {% endembed %} After running the simulation the aerodynamic coefficients can be plotted against the angle of attack $$\alpha$$. The image below shows the $$C\_{X}$$,$$C\_{Z}$$, and $$C\_{m}$$aerodynamic coefficients. ![Aerodynamic coefficients versus angle of attack](https://2298927434-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-MWUM1iuDAQw73diDpYu%2F-MfKI0XI_pwaBk9mOKsm%2F-MfKJ_8TQLhtiMNOfLJk%2Fcoefficients.png?alt=media\&token=e81375d9-236a-48c5-8c55-808ee43a6676) The locations of the minimum and maximum $$\alpha$$ and $$q$$are shown in the images. The $$C\_{X}$$coefficient exhibits a quadratic behavior, which is expected. The$$C\_{X0}$$, $$C\_{X,\alpha}$$, $$C\_{X,\bar{q}}$$, $$C\_{Z0}$$, $$C\_{Z,\alpha}$$, $$C\_{Z,\bar{q}}$$,$$C\_{m0}$$, $$C\_{m,\alpha}$$, and $$C\_{m,\bar{q}}$$stability derivatives can be obtained in multiple ways. Here the $$C\_{Z}$$and $$C\_{m}$$stability derivatives are obtained with the least squares method. To obtain the$$C\_{m}$$stability derivatives the following system of equations is defined: $$ \begin{bmatrix} \mathbf{1} & \theta(t) & q(t)\frac{c\_{ref}}{2|V\_{ref}|} \end{bmatrix}\begin{Bmatrix}C\_{m0}\C\_{m,\alpha}\C\_{m,\bar{q}}\end{Bmatrix}=C\_{m}(t). $$ Note that in the above system of equations$$\mathbf{1}$$is a vector of ones which length is equal to the length of the $$\theta$$and the $$q$$vectors. The system is overdetermined - it has more equations than unknowns. However the least squares method can be used to solve the overdetermined system. Solving the system gives the values for $$C\_{m0}$$, $$C\_{m,\alpha}$$, and $$C\_{m,\bar{q}}$$. The same process is repeated for the $$C\_{Z}$$aerodynamic coefficient. A polynomial fitting method is used for $$C\_{X}$$in order to capture its quadratic behavior. A third order polynomial in $$\alpha$$is selected. The polynomial is given by the following equation: $$ C\_{X}(\alpha)=C\_{X0}+C\_{X,\alpha}\alpha+C\_{X,\alpha^2}\alpha^2. $$ Unlike the least squares method the polynomial fitting method cannot be used to determine $$C\_{X,\bar{q}}$$. The $$C\_{X,\bar{q}}$$can be determined from the following equation: $$ C\_{X,\bar{q}}=\frac{C\_{X}\rvert\_{q\_{max}}-C\_{X}\rvert\_{q\_{min}}}{2kA}, $$ where $$k$$is the reduced frequency, given by $$\omega\frac{c{ref}}{2|V\_{ref}|}$$. The table below shows the obtained stability derivatives: | | $$\_{0}$$ | $$\frac{\partial}{\partial\alpha}$$ | $$\frac{\partial}{\partial\alpha}$$ | $$\frac{\partial}{\partial\bar{q}}$$ | | ---------- | --------- | ----------------------------------- | ----------------------------------- | ------------------------------------ | | $$C\_{X}$$ | -0.0219 | 0.2595 | 3.1367 | -0.2831 | | $$C\_{Z}$$ | -0.3149 | -4.9830 | | 5.9714 | | $$C\_{m}$$ | 0.0458 | -1.3909 | | -19.2330 | The wing-tail design is both statically and dynamical stable as$$C\_{m,\alpha}$$and $$C\_{m,\bar{q}}$$are negative. ### Verifying the obtained stability derivatives The time history of the aerodynamic coefficients can be used to verify the obtained stability derivatives. The stability derivatives are used to describe the time history of $$C\_{X}$$, $$C\_{Z}$$, and $$C\_{m}$$: $$ C\_{X}(t)=C\_{X0}+C\_{X,\alpha}\alpha(t)+C\_{X,\alpha^2}\alpha(t)^2+C\_{X,\bar{q}}q(t)\frac{c\_{ref}}{2|V\_{ref}|},\\ C\_{Z}(t)=C\_{Z0}+C\_{Z,\alpha}\alpha(t)+C\_{Z,\bar{q}}q(t)\frac{c\_{ref}}{2|V\_{ref}|},\\ C\_{m}(t)=C\_{m0}+C\_{m,\alpha}\alpha(t)+C\_{m,\bar{q}}q(t)\frac{c\_{ref}}{2|V\_{ref}|}. $$ The above equations can be compared against $$C\_{X}(t)$$, $$C\_{Z}(t)$$, and $$C\_{m}(t)$$obtained from the simulation. The variation of angle of attach $$\alpha$$and the pitch rate $$q$$with respect to time is known and is given by: $$ \theta=\alpha=Asin(\omega t),\\ q=\dot{\theta}=\omega A cos(\omega t). $$ The image below shows the comparison of $$C\_{X}$$, $$C\_{Z}$$, and $$C\_{m}$$: ![Comparison of the aerodynamic coefficients from the simulation and from the obtained aerodynamic derivatives](https://2298927434-files.gitbook.io/~/files/v0/b/gitbook-legacy-files/o/assets%2F-MWUM1iuDAQw73diDpYu%2F-MfK95BUVa6Zf94xQzjt%2F-MfKHYJDmCupxw2tB3E9%2Ffit.png?alt=media\&token=71b2af44-fc91-4bad-9f4a-c04e2433283c) As seen from the image the aerodynamic coefficients agree well. The minor differences from $$t=0\text{ }s$$to $$t=0.03\text{ }s$$are due to the simulation transients (the time it takes for the flow to develop).