Stick Model Lumping: Method

Stick Model Lumping: Lumped Masses, Harmonic-Mean Stiffness, and the Myklestad and Holzer Frequency Passes

The problem

Before a wing, a tail or a fuselage is modelled with shell and beam elements, and long before a flutter analysis, its dynamics are first captured by a stick model: a line of point masses at stations along the elastic axis joined by massless elastic segments. The stick model is what Myklestad (1944) and Holzer (1921) computed by hand, what NACA TN 1522 (Houbolt and Anderson 1948) worked through example by example for nonuniform beams, and what every beam-element dynamic model still is underneath. Building it from the distributed data of a loads report, running mass and stiffness against span, raises three questions that this page answers in a fixed, checkable way: how the running mass is turned into point masses without losing the mass or moving the centre of gravity, what constant stiffness a tapering segment should be given, and what the resulting uncoupled bending and torsion frequencies are, found by the two classical transfer passes and confirmed by an eigen-solution of the same model.

Conventions and units

The span coordinate y runs from the root (station 0, y = 0) to the tip (station N, y = L) along the elastic axis; the data rows give EI (bending about the chordwise axis), GJ, the running weight w, the chordwise offset e of the running mass centroid from the elastic axis (aft positive) and its radius of gyration r about that centroid, and concentrated weights W_c with their own offsets e_c and radii r_c, all varying linearly between rows. Weights are converted to mass with g = 9.80665 m/s² (386.088 in/s²). Lengths, running weights and weights go through the unit system (in, lb/in, lb or m, N/m, N); EI and GJ are plain numbers in lb-in² or N-m² because the dictionary has no stiffness unit, and they are not changed by the unit switch. Masses are reported in the consistent dynamic units lb·s²/in (the "snail", weight in lb divided by 386.088) or kg, inertias in lb·s²·in or kg·m², frequencies in Hz and rad/s. Deflection w is positive up, the bending moment is M = EI d²w/dy² (positive bends the tip up), the shear is V = dM/dy, so that a tip load P in the +w direction gives V = −P along the beam; twist φ is positive nose up about the elastic axis and the torque is T = GJ dφ/dy.

Lumping the running mass

Two rules are offered. The tributary-length rule, the one of the handbooks and of Myklestad's own examples, gives station k the running mass between the mid-points of its two neighbouring segments, m_k = ∫ μ dy from ½(y_k−1 + y_k) to ½(y_k + y_k+1), the root station from y = 0 to the first mid-point and the tip station from the last mid-point to the tip; a concentrated mass goes to the nearest station. The total mass is preserved exactly; the first moment is not, so the centre of gravity moves by a fraction of a segment length, which the page prints. The lever rule (the statically equivalent or "consistent" lumping) splits every piece of running mass between the two stations of its segment in inverse proportion to the distances from its centroid, m_k−1 += m_p(y_k − ȳ_p)/l_s and m_k += m_p(ȳ_p − y_k−1)/l_s, and a concentrated mass between stations likewise; mass and first moment are both exact. The chordwise static moment m e and the mass moment of inertia about the elastic axis m (e² + r²) of each lump are integrated over the same lengths (or split with the same fractions), so that the station carries, for a later coupled or torsional analysis, the inertia of the mass it represents: I_θ,k = ∫ μ (e² + r²) dy + Σ m_c(e_c² + r_c²). The integrals are evaluated by three-point Gauss quadrature on every linear piece, which is exact for the cubic integrand μ (e² + r²). The root lump is listed for the mass check but does not move: in a cantilever it carries no kinetic energy, and in a finite-element model it is simply grounded.

Both rules converge to the continuous beam as the stations are refined; the tributary rule gives slightly higher first-mode frequencies than the lever rule on the same stations because it keeps the mass nearer the root for a tapering distribution. For the uniform cantilever lumped by the tributary rule the first three frequencies are 1.8 %, 5.9 % and 9.3 % low with 5 stations, 0.46 %, 1.6 % and 2.6 % low with 10, 0.11 %, 0.40 % and 0.65 % low with 20 and 0.03 %, 0.10 % and 0.16 % low with 40, converging like 1/N² from below; the torsion modes of the uniform shaft converge like 1/N² too but from a lower start (10 stations: 0.10 %, 0.92 % and 2.6 % low). A rule of thumb is ten stations per mode wanted to one percent.

Segment stiffness

Between two stations the beam becomes a massless segment of constant EI and GJ. The constant that preserves the segment's flexibility, the relative rotation of its ends under a constant moment ∫ M dy/EI, is the harmonic mean EI_s = l_s / ∫ dy/EI(y). For EI varying linearly from EI_a to EI_b over the segment the integral is exact, ∫ dy/EI = l_s ln(EI_b/EI_a)/(EI_b − EI_a), so EI_s is the logarithmic mean (EI_b − EI_a)/ln(EI_b/EI_a), always below the arithmetic mean; a segment that spans several data rows is integrated piece by piece. The arithmetic mean would make a tapering segment too stiff (by 1.5 % for a 2:1 taper within one segment, 12 % for 10:1). The same mean serves GJ for the torsional flexibility ∫ T dy/GJ. The segment EI is exact for the end rotation but not for the deflection under a varying moment within the segment, which is one more reason the stations should be close enough that EI does not change by more than a factor of about two across any one of them.

Myklestad's bending pass

For harmonic motion at circular frequency ω the Euler–Bernoulli beam equation (EI w″)″ = μ ω² w is written as four first-order relations for the state (w, θ, M, V): θ = w′, M = EI θ′, V = M′, V′ = μ ω² w. Along a massless segment the inertia term vanishes, so the state at the outboard end follows from the static beam formulas: w⁺ = w + θ l + M l²/(2EI) + V l³/(6EI), θ⁺ = θ + M l/EI + V l²/(2EI), M⁺ = M + V l, V⁺ = V (Myklestad's field matrix). Across a point mass the deflection, slope and moment are continuous and the shear jumps by the inertia force, V⁺ = V + m ω² w (the point matrix). At a clamped root w = θ = 0 and the moment and shear are unknown, so the root state is a combination a (0, 0, 1, 0) + b (0, 0, 0, 1); at a root restrained by a rotational spring K_θ the slope is free but M = K_θ θ, and the basis is (0, 1, K_θ, 0) and (0, 0, 0, 1). Both basis vectors are carried station by station to the tip, where the free end requires M = 0 and V = 0 after the tip mass; the condition that a non-trivial (a, b) exists is the vanishing of the 2 × 2 determinant Δ(ω) = M_A V_B − M_B V_A. The pass is Myklestad's 1944 tabular method written as a transfer-matrix recursion (Pestel and Leckie 1963); Houbolt and Anderson's TN 1522 reaches the same frequencies by Stodola iteration of the deflection shape and is the public-domain record of the method's accuracy on nonuniform beams.

Two numerical points. The components of the state grow by many orders of magnitude along a long beam at a high trial frequency, so both basis vectors are divided by the same scale factor at every station (the ratio of the two is all the determinant needs) and the determinant is normalised by the lengths of its two rows, Δ̂ = Δ / (|(M_A, M_B)| |(V_A, V_B)|), which lies between −1 and 1 whatever the units and is the quantity drawn in the page's figure. The sweep then runs over a uniform grid of 800 trial frequencies up to 15 % above the highest mode wanted, locates the sign changes of Δ̂ and refines each by bisection to twelve digits; the determinant of a lumped chain is a polynomial in ω², so it has no poles to be mistaken for roots (the Holzer residual torque is likewise polynomial). The mode shape at a root is the null vector (a, b) = (M_B, −M_A) applied to the stored states, normalised to a unit tip deflection, and the page tabulates w, θ, M and V at every station from it; the shear tabulated is the value just inboard of the station mass.

Cross-check by the flexibility matrix. The same lumped model has the cantilever influence coefficients δ_ij = ∫₀^{min(y_i, y_j)} (y_i − y)(y_j − y)/EI dy, integrated exactly over the piecewise-constant segments (plus y_i y_j/K_θ for the root spring), and the eigenvalues λ of the symmetric matrix M^½ δ M^½ give every natural frequency as ω = 1/√λ (Jacobi rotations). The page prints these beside the Myklestad roots: the two must agree to the digits shown, since they are two solutions of one discrete model, and the comparison is the guarantee that the sweep found every root in its range. The test script filed with the tool also checks both against an independent numpy transfer-matrix pass and a Hermite beam finite-element model with the rotations condensed, on the default case and on the TN 1522 examples.

Holzer's torsion pass

The torsional equation (GJ φ′)′ = −I_θ ω² φ gives the state (φ, T) with T = GJ φ′: along a massless segment φ⁺ = φ + T l/GJ and T⁺ = T; at a station of inertia I_θ the torque jumps by the inertia torque, T⁺ = T − I_θ ω² φ. Holzer's method starts at the fixed root with φ = 0 and a unit torque, carries the pair to the tip and reads the residual torque there, which must vanish at a natural frequency; the page normalises it by the largest torque along the span, sweeps and bisects as for bending, and prints the eigenvalues of the torsional flexibility matrix δ_ij = Σ l_s/GJ_s (segments inboard of the nearer station) beside the roots. The station inertias are the lumped m (e² + r²) about the elastic axis: the chordwise offsets of the masses enter here and nowhere else, because the two passes are uncoupled.

Limits

The passes are uncoupled: a mass centroid off the elastic axis couples bending and torsion in reality, and the coupled frequencies of a wing with a heavy forward engine differ from the uncoupled ones; the page reports the static unbalance Σ m e and the inertias so that a coupled analysis (a two-state-vector Myklestad pass, or the beam-element model written from the paste-ready table) can follow. Shear deformation and rotary inertia are neglected (Euler–Bernoulli), which overestimates the higher modes of deep, short beams; chordwise bending, sweep, dihedral and a free–free airplane are not represented, and the root is either clamped or restrained by a single rotational spring (the Stodola examples of TN 1522 for beams on translational springs and free–free beams are outside this tool). The frequencies are those of the lumped model: their distance from the continuous beam's is the lumping error, which the N-stations option lets the user measure by refining the stations.

Reading the results

The distributed-input section lists the rows as read and the exact integrals of the total weight, the mass, the spanwise and chordwise centre of gravity, the static moment and the inertia about the elastic axis. The lumped-stations section draws the stick, tabulates every station's tributary length, running and concentrated weight, mass, chordwise centroid, static moment and inertia, and checks the lumps against the input (the c.g. shift of the tributary rule is printed); Figure 2 lays the lumped running weight over the entered curve. The segment-stiffness section prints the harmonic-mean formula with the first segment's numbers and the table of all segments, with Figure 3. The Myklestad section tabulates the frequencies against the flexibility eigenvalues with the generalised masses, draws the normalised determinant and the mode shapes, and tabulates w, θ, M and V per station; the Holzer section does the same for torsion. The finite-element section gives the GRID, CONM2 and CBAR / PBAR (or segment table) text in the current units, ready to paste.