Lateral modes: roll, spiral and Dutch roll

Sideways motion has its own family of modes, and they are more tangled than the longitudinal ones because roll and yaw are coupled through the dihedral effect, adverse yaw and the yawing moment due to roll rate. In the simulator Sam probes each one. A sharp aileron input followed by release shows how quickly the roll rate builds and stops: the roll mode. Banking 10 degrees and letting go shows whether the bank slowly grows or shrinks over half a minute: the spiral mode. A rudder doublet sets the nose wagging from side to side while the wings rock: the Dutch roll. This lesson builds the lateral-directional linear model of the Tern Block 2 in cruise and finds all three.

The lateral state-space model

The lateral equations come from the side-force, rolling and yawing equations of Lesson 1, linearized in the same way as the longitudinal ones. Using sideslip \(\beta\) instead of side velocity, with state \([\beta, p, r, \phi]\), the free motion obeys

Here the primed notation often used in textbooks, which folds in the product of inertia \(I_{xz}\), is dropped because the Tern's stability-axis \(I_{xz}\) is small. The new derivatives are \(L_p\), roll damping, the rolling moment opposing roll rate because the down-going wing meets the air at a higher angle of attack; \(N_r\), yaw damping, mainly from the fin; \(L_r\), the rolling moment due to yaw rate, because the outer wing in a yaw moves faster and lifts more; and \(N_p\), the yawing moment due to roll rate, usually adverse. Each is converted from its non-dimensional coefficient using \(b/2u_0\) for the rate terms.

# Lateral-directional modes of the Tern Block 2 in cruise: roll subsidence, spiral and
# Dutch roll from the full model, compared with the classic approximations.
import math
def charpoly(A):
    """Coefficients of det(sI - A), highest power first (Faddeev-LeVerrier)."""
    n = len(A); M = [[0.0] * n for _ in range(n)]; c = [1.0]
    I = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
    for k in range(1, n + 1):
        AM = [[sum(A[i][l] * M[l][j] for l in range(n)) for j in range(n)] for i in range(n)]
        M = [[AM[i][j] + c[-1] * I[i][j] for j in range(n)] for i in range(n)]
        AM = [[sum(A[i][l] * M[l][j] for l in range(n)) for j in range(n)] for i in range(n)]
        c.append(-sum(AM[i][i] for i in range(n)) / k)
    return c

def roots(c, it=500):
    """All roots of a polynomial by the Durand-Kerner iteration."""
    n = len(c) - 1; a = [x / c[0] for x in c]
    z = [complex(0.4, 0.9) ** k for k in range(n)]
    for _ in range(it):
        new = []
        for i in range(n):
            p = sum(a[k] * z[i] ** (n - k) for k in range(n + 1))
            d = 1
            for j in range(n):
                if j != i: d *= (z[i] - z[j])
            new.append(z[i] - p / d)
        z = new
    return sorted(z, key=lambda r: -abs(r))

def describe(s):
    import math
    if abs(s.imag) < 1e-6:
        t = -1 / s.real if s.real < 0 else 1 / s.real
        kind = 'time constant' if s.real < 0 else 'doubles every'
        val = t if s.real < 0 else math.log(2) / s.real
        return f"real root {s.real:+.4f}: {kind} {val:.2f} s"
    wn = abs(s); z = -s.real / wn
    return f"pair {s.real:+.4f} +/- {abs(s.imag):.4f}j: wn {wn:.3f} rad/s, zeta {z:.3f}, period {2*math.pi/abs(s.imag):.1f} s"

m, Ixx, Izz, S, b, g, rho, u0 = 6500.0, 30000.0, 52000.0, 28.0, 17.0, 9.81, 0.909, 105.0
QS = 0.5 * rho * u0 ** 2 * S; k = b / (2 * u0)
d = dict(Cyb=-0.75, Clb=-0.099, Cnb=0.109, Clp=-0.50, Cnp=-0.040, Clr=0.120, Cnr=-0.130)
def model(d):
    Yb = d['Cyb'] * QS / m
    Lb, Lp, Lr = (d[x] * QS * b / Ixx * f for x, f in (('Clb', 1), ('Clp', k), ('Clr', k)))
    Nb, Np, Nr = (d[x] * QS * b / Izz * f for x, f in (('Cnb', 1), ('Cnp', k), ('Cnr', k)))
    A = [[Yb / u0, 0.0, -1.0, g / u0],
         [Lb, Lp, Lr, 0.0],
         [Nb, Np, Nr, 0.0],
         [0.0, 1.0, 0.0, 0.0]]
    return A, Yb, Lb, Lp, Lr, Nb, Np, Nr
A, Yb, Lb, Lp, Lr, Nb, Np, Nr = model(d)
print("state x = [beta, p, r, phi]; A =")
for row in A: print("  " + "  ".join(f"{v:+8.4f}" for v in row))
r = roots(charpoly(A))
for s in r:
    if s.imag >= -1e-9: print("  " + describe(s))
# approximations
tau_roll = -1 / Lp
s_spiral = (Lb * Nr - Lr * Nb) / Lb               # positive means a slowly divergent spiral
wn_dr = math.sqrt(Nb + Yb * Nr / u0)
z_dr = -(Yb / u0 + Nr) / (2 * wn_dr)
print(f"\napproximations: roll time constant {tau_roll:.2f} s, spiral root {s_spiral:+.4f} per s")
print(f"Dutch roll wn {wn_dr:.2f} rad/s, zeta {z_dr:.3f}")
print(f"spiral condition Lb*Nr - Lr*Nb = {Lb*Nr - Lr*Nb:+.3f} (positive gives a stable spiral)")
# what would make the spiral stable? more dihedral (Cl_beta) or more yaw damping (Cn_r)
for name, chg in (('Cl_beta -0.120', dict(Clb=-0.120)), ('Cn_r -0.160', dict(Cnr=-0.160)), ('Cn_beta 0.140', dict(Cnb=0.140))):
    dd = dict(d); dd.update(chg)
    rr = [s for s in roots(charpoly(model(dd)[0])) if abs(s.imag) < 1e-6]
    sp = max(rr, key=lambda s: s.real)
    dr = [s for s in roots(charpoly(model(dd)[0])) if s.imag > 1e-6][0]
    print(f"  with {name}: spiral root {sp.real:+.4f}, Dutch roll zeta {-dr.real/abs(dr):.3f}")

state x = [beta, p, r, phi]; A =
   -0.1542   +0.0000   -1.0000   +0.0934
   -7.8711   -3.2181   +0.7723   +0.0000
   +4.9997   -0.1485   -0.4827   +0.0000
   +0.0000   +1.0000   +0.0000   +0.0000
  real root -3.3280: time constant 0.30 s
  pair -0.2636 +/- 2.3271j: wn 2.342 rad/s, zeta 0.113, period 2.7 s
  real root +0.0003: doubles every 2184.37 s

approximations: roll time constant 0.31 s, spiral root +0.0079 per s
Dutch roll wn 2.25 rad/s, zeta 0.141
spiral condition Lb*Nr - Lr*Nb = -0.062 (positive gives a stable spiral)
  with Cl_beta -0.120: spiral root -0.0037, Dutch roll zeta 0.106
  with Cn_r -0.160: spiral root -0.0042, Dutch roll zeta 0.135
  with Cn_beta 0.140: spiral root +0.0047, Dutch roll zeta 0.103

The roll mode

The roll mode, also called roll subsidence, is a fast, heavily damped real root. After an aileron input, roll rate builds exponentially toward a steady value set by the balance between aileron moment and roll damping. Because roll damping is so strong, a one-degree-of-freedom model is enough:

For the Tern the program gives a time constant of 0.31 seconds from the approximation and 0.30 seconds from the full model. Roll rate reaches 63 percent of its final value within a third of a second: the aircraft feels crisp in roll. A long wing raises roll damping strongly, which is why long-span gliders, and to a lesser extent the high-aspect-ratio Tern, reach only modest steady roll rates even though they respond without lag.

The spiral mode

The spiral is a slow real root, stable or unstable, describing what happens to a small bank angle if the pilot does nothing. Bank tilts the lift and the aircraft starts to sideslip toward the low wing. Two effects then compete. The dihedral effect rolls the aircraft back toward wings level. Directional stability yaws the nose toward the sideslip, into the turn, and the faster outer wing then lifts more through \(L_r\), steepening the bank. If the second effect wins, the bank and turn rate grow slowly and the nose drops: a spiral dive. The approximate root is

Strong dihedral effect and strong yaw damping stabilize the spiral; strong directional stability destabilizes it. For the Tern the condition is slightly negative, −0.062, and the approximation predicts a mildly divergent spiral doubling its bank every 88 seconds. The full model, which includes the gravity coupling more exactly, puts the root almost on the origin: it doubles only every 2,200 seconds. Either way, the spiral is so slow that a pilot corrects it without thinking, and handling-qualities criteria allow a mildly unstable spiral provided its time to double is long, at least 20 seconds for level 1 in the categories used in Lesson 9.

The Dutch roll

The Dutch roll is a lightly damped oscillation combining yaw, sideslip and roll. The nose swings left and right under the weathercock stiffness, and through the dihedral effect each sideslip rolls the aircraft, so the wingtips trace an ellipse as seen from behind. The name is usually attributed to the motion of a skater on Dutch canals. If directional stability were the only restoring effect and yaw damping the only damping, the Dutch roll would be a simple yaw oscillation:

This approximation gives 2.25 rad/s and a damping ratio of 0.14, against 2.34 rad/s and 0.11 from the full model. It is noticeably worse than the short-period approximation because the roll coupling, which it ignores, is significant: rolling motion in the Dutch roll interacts with the dihedral effect and reduces the damping. The ratio of roll to sideslip in the mode grows with \(C_{l_\beta}/C_{n_\beta}\), and when that ratio is large the motion becomes roll-dominated and uncomfortable. For the Tern, with a ratio close to one, the wings rock about as much as the nose wags.

The last part of the program shows the classic design tension. More dihedral effect, \(C_{l_\beta}\) of −0.120, stabilizes the spiral but slightly reduces Dutch roll damping. More yaw damping, \(C_{n_r}\) of −0.160, helps both. More directional stability, \(C_{n_\beta}\) of 0.140, makes the spiral divergent and reduces the Dutch roll damping ratio. Only yaw damping improves everything, and the most effective way to add yaw damping without a larger fin is a yaw damper: an automatic rudder input proportional to yaw rate, the subject of Lesson 10.

Measuring the lateral modes in flight

Each mode has its own flight test input. For the Dutch roll, the pilot applies a rudder doublet, a short pulse one way and then the other, and releases the controls; the period is read from the yaw rate trace and the damping from the ratio of successive peaks, using the logarithmic decrement. For the roll mode, a step aileron input is held and the rise of roll rate fitted with an exponential. For the spiral, the aircraft is trimmed wings level, banked 10 or 20 degrees with aileron, and the controls released; the time for the bank to double or halve is timed with a stopwatch or read from recorded data. Sam will fly each test at forward and aft CG, light and heavy, clean and with flaps, and the results will be compared with the linear model's predictions. Where they disagree, the derivatives will be adjusted by system identification, a fitting process that finds the derivative values that best reproduce the measured responses.