AE 454/777 Dynamics & Bifurcation
Tutorial 1
Deadline: 20th Feburary 2026
Instructions
Please adhere to the following guidelines for your submission:
• File Format: Submit your solutions as a single PDF report named T1_<roll_number>.pdf.
• Programming Languages: You may use MATLAB, Python, or Julia.
• Code Listing: Include your full code within the report. You can either paste the source code
directly or append a PDF export of your notebook to the end of your report.
• AI Policy: While there are no formal restrictions on the use of AI, you are encouraged
to develop the integrators independently. Please limit AI usage to syntax assistance and
debugging to ensure a deep understanding of the algorithms.
Prelude
Aircraft wings are subjected to aerodynamic moments that tend to twist the structure under load.
The inherent elasticity of the wing provides a restoring torsional stiffness, allowing the wing to be
modelled, in a simplified sense, as a torsional spring. As the wing deforms, the aerodynamic loads
acting on it change in response to the altered geometry, introducing a coupled interaction between
aerodynamics and structural elasticity. This interaction can give rise to aeroelastic instabilities. In
this tutorial, we examine a toy model to analyse the balance of moments on a wing and to study
the resulting static equilibrium behaviour.
Moments on the wing
In real planes, twisting and bending are coupled, but here we isolate the twisting. As illustrated in
figure 1, assume the wing to be a rod anchored to a torsional spring of stiffness k at its shear centre.
Let α be the average twist angle of the wing. The aerodynamic moment about shear centre is given
by M (α) = q∞ SCm (α) where q∞ = 0.5ρ∞ V∞ 2 , S is wing area and moment coefficient defined as:
Cm (α) = Cm0 + Cmα sin α cos α.
For a twist of α, the potential energy of the system is given by
1 q∞ SCmα
Π(α) = kα2 − sin2 (α) − q∞ SCm0 α (1)
2 2
Note: Assume that centre of gravity aligns with shear centre, hence we neglect moment due to
wing weight. Now:
1
Figure 1: Left: Twisting of wing due to aerodynamic moments, Right: Simplified 2D toy model
schematic
1. Show that equation for α at equilibrium can be written in dimensionless form as
µ0 + µ1 sin(2α) − α = 0, (2)
for appropriate choices of µ0 and µ1 . (Hint: The equilibrium equation is obtained from
−Π′ (α) = 0)
2. For µ0 = 0, give graphical analysis of the dimensionless equation for the cases µ1 > 0.5 and
µ1 < 0.5. How many equilibria can exist in each case?
3. Let r = 2µ1 − 1, and h = µ0 . Give a numerically accurate plot of the bifurcation curves in the
(r, h) plane. Interpret your results physically.
Remark: µ0 is directly related to the camber of the wing airfoil and µ1 is proportional to δ or the
separation between shear centre and centre of pressure. (see figure 1)
Dynamics of wing twisting
In the previous section, we analysed the static equilibriums and conservative moments acting on the
wing. Now the complete dynamic equation of evolution of α(t) is given by
I α̈ + ζ α̇ + Π′ (α) = 0,
where I is the wing moment of inertia and ζ is the damping coefficient. For small α and Cm0 = 0,
show that the above dynamic equation can be reduced to:
4µ1 2
I α̈ + ζ α̇ = kα r − α + O(α5 ), (3)
3
where r = 2µ1 − 1.
2
Amplitude dynamics
Moreover, solution for α(t) can be approximated as α(t) ≈ A(t) sin(ωt + ϕ(t)). Derivations of
dynamics of A and ϕ are beyond the scope of the current tutorial. Henceforth, taking the solution
form for granted and with appropriate variable change, we can obtain:
Ȧ = A(1 − A2 ) (4)
1. Solve the initial value problem for the ODE in (4) with A(0) = 0.5 analytically. What is the
value of A(1)?
2. Using the Euler method with step size ∆t = 1, estimate A(1) numerically - call this Â(1).
Then repeat using ∆t = 10−n where n = 1,2,3,4.
3. Plot the error E = |A(1) − Â(1)| as a function of ∆t. Then plot ln(E) against ln(∆t). Explain
the results.
4. Redo this exercise for
(a) Runge-Kutta - 2 method
(b) two-step Adams-Bashforth method
Phase dynamics
The evolution of ϕ(t) ∈ [0, 2π) often resemble that of a highly damped non-linear pendulum. The
dynamics equation is given by:
bϕ̇ + Γ1 sin(ϕ) = Γ0 (5)
where b is the damping coefficient and Γ0 is constant applied torque.
1. Reduce this equation to the form of a non-uniform oscillator ϕ′ = γ − sin ϕ with appropriate
values for t and γ.
2. Sketch sin ϕ(t) vs. t for different γ , including the limiting cases γ ≈ 1 and γ << 1.
3. Redo the exercise for ϕ′ (t) instead of sin ϕ(t).