Links

Obtaining stability derivatives from forced oscillations

This article explains how to obtain the stability derivatives of a design from forced oscillations
Consider a simple wing-tail design. The design has a wing span
bref=1.7321 mb_{ref}=1.7321\text{ }m
, wing area
Sref=0.3 m2S_{ref}=0.3\text{ }m^2
, and wing chord of
cref=0.1732 mc_{ref}=0.1732\text{ }m
. The image below shows the mesh of the design.
Mesh of a wing-tail design
A body-fixed coordinate system is placed at the leading edge of the wing. The reference speed is
Vref=25 m/s|V_{ref} |=25\text{ }m/s
and the reference density, pressure and viscosity are evaluated at an altitude of
h=0 mh=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 for more information on how to determine the neutral point.
The table below shows the
CZC_{Z}
and
CmC_{m}
values used for determining the location of the neutral point,
xnpx_{np}
.
Simulation type
α\alpha^\circ
CZC_{Z}
CmC_{m}
Fixed-wake
0
−0.3215
-0.0223
Fixed-wake
5
-0.7913
-0.1973
The location of the neutral point,
xnpx_{np}
, can be obtained from the following equation:
xnp=crefCmαCZαΔCmΔCZCmα=5Cmα=0CZα=5CZα=0.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
CZC_{Z}
and
CmC_{m}
values from the table into the above equation,
xnpx_{np}
evaluates to
0.0645 m-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,
SMSM
:
SM=xcgxnpcrefSM=\frac{x_{cg}-x_{np}}{c_{ref}}
For
SM=0.1 (10%)SM=0.1\text{ }(10\%)
,
xcg=0.0385 mx_{cg}=-0.0385\text{ }m
.
For a statically stable design the centre of gravity is ahead of the neutral point.
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 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 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:
θ=Asin(ωt),\theta=Asin(\omega t),
where
θ\theta
is the pitch angle,
AA
is the amplitude of the pitch oscillation,
ω\omega
is the angular frequency of the oscillation (
ω=2πf\omega=2\pi f
), and
tt
is the time. In this example
A=5 degA=5\text{ }deg
and
f=5Hzf=5 Hz
. The pitch rate
qq
can be obtained by differentiating the
θ\theta
equation with respect to time:
q=θ˙=ωAcos(ωt).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×103dt=4\times10^{-3}
and
Ntimesteps=100N_{timesteps}=100
. The image below shows
θ\theta
and
qq
. Note that
qq
is divided by 100 so the difference between
θ\theta
and
qq
is easy to see. In the actual motion file
qq
is not divided by 100.
Generated pitch angle and pitch rate

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.
Wing-tail design forced pitch oscillation
After running the simulation the aerodynamic coefficients can be plotted against the angle of attack
α\alpha
. The image below shows the
CXC_{X}
,
CZC_{Z}
, and
CmC_{m}
aerodynamic coefficients.
Aerodynamic coefficients versus angle of attack
The locations of the minimum and maximum
α\alpha
and
qq
are shown in the images. The
CXC_{X}
coefficient exhibits a quadratic behavior, which is expected.
The
CX0C_{X0}
,
CX,αC_{X,\alpha}
,
CX,qˉC_{X,\bar{q}}
,
CZ0C_{Z0}
,
CZ,αC_{Z,\alpha}
,
CZ,qˉC_{Z,\bar{q}}
,
Cm0C_{m0}
,
Cm,αC_{m,\alpha}
, and
Cm,qˉC_{m,\bar{q}}
stability derivatives can be obtained in multiple ways. Here the
CZC_{Z}
and
CmC_{m}
stability derivatives are obtained with the least squares method. To obtain the
CmC_{m}
stability derivatives the following system of equations is defined:
[1θ(t)q(t)cref2Vref]{Cm0Cm,αCm,qˉ}=Cm(t). \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
1\mathbf{1}
is a vector of ones which length is equal to the length of the
θ\theta
and the
qq
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
Cm0C_{m0}
,
Cm,αC_{m,\alpha}
, and
Cm,qˉC_{m,\bar{q}}
. The same process is repeated for the
CZC_{Z}
aerodynamic coefficient.
A polynomial fitting method is used for
CXC_{X}
in order to capture its quadratic behavior. A third order polynomial in
α\alpha
is selected. The polynomial is given by the following equation:
CX(α)=CX0+CX,αα+CX,α2α2.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
CX,qˉC_{X,\bar{q}}
. The
CX,qˉC_{X,\bar{q}}
can be determined from the following equation:
CX,qˉ=CXqmaxCXqmin2kA,C_{X,\bar{q}}=\frac{C_{X}\rvert_{q_{max}}-C_{X}\rvert_{q_{min}}}{2kA},
where
kk
is the reduced frequency, given by
ωcref2Vref\omega\frac{c{ref}}{2|V_{ref}|}
. The table below shows the obtained stability derivatives:
0_{0}
α\frac{\partial}{\partial\alpha}
α\frac{\partial}{\partial\alpha}
qˉ\frac{\partial}{\partial\bar{q}}
CXC_{X}
-0.0219
0.2595
3.1367
-0.2831
CZC_{Z}
-0.3149
-4.9830
5.9714
CmC_{m}
0.0458
-1.3909
-19.2330
The wing-tail design is both statically and dynamical stable as
Cm,αC_{m,\alpha}
and
Cm,qˉ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
CXC_{X}
,
CZC_{Z}
, and
CmC_{m}
:
CX(t)=CX0+CX,αα(t)+CX,α2α(t)2+CX,qˉq(t)cref2Vref,CZ(t)=CZ0+CZ,αα(t)+CZ,qˉq(t)cref2Vref,Cm(t)=Cm0+Cm,αα(t)+Cm,qˉq(t)cref2Vref.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
CX(t)C_{X}(t)
,
CZ(t)C_{Z}(t)
, and
Cm(t)C_{m}(t)
obtained from the simulation. The variation of angle of attach
α\alpha
and the pitch rate
qq
with respect to time is known and is given by:
θ=α=Asin(ωt),q=θ˙=ωAcos(ωt).\theta=\alpha=Asin(\omega t),\\ q=\dot{\theta}=\omega A cos(\omega t).
The image below shows the comparison of
CXC_{X}
,
CZC_{Z}
, and
CmC_{m}
:
Comparison of the aerodynamic coefficients from the simulation and from the obtained aerodynamic derivatives
As seen from the image the aerodynamic coefficients agree well. The minor differences from
t=0 st=0\text{ }s
to
t=0.03 st=0.03\text{ }s
are due to the simulation transients (the time it takes for the flow to develop).