Control Oriented Modeling, Experimentation, and Stability Analysis of an Autorotating Samara

Samara seed pods, a morphology that has evolved parallel in numerous plants throughout the world, effectively employ autorotation to slow descent speed. This allows organisms to produce heavier seeds/fruit that can be scattered across a larger area by prevailing winds and gusts, [1–6]. The study of these biological structures can provide insight into efficient design of aerodynamic systems, both small and large. Samaras can be classified by wing configuration and aerodynamic behavior in descent [7]. This paper will focus on samars of the single-winged variety that do not roll about their spanwise axis at steady-state.

A seminal work by Norberg [8], has presented a thorough qualitative and experimental analysis of samara stability. However, through dynamic modeling, a mathematical analysis can be performed to draw more quantitative insight in regards to samara descent behavior, e.g., ranges of physical and aerodynamic parameters that allow for stable autorotation. A simple mathematical analysis of samara stability was presented in Ref. [9] for the purpose of modeling and controlling a powered single-winged rotorcraft. This work was further expanded upon in Ref. [10]. Computational fluid dynamics was utilized in Refs. [11] and [12] to analyze the microscale effects of turbulence and leading-edge vortices. Discussion of the effects of leading-edge vortices and the robustness of samara stability to gusts has been given in Ref. [13]. An extensive and detailed model of samara dynamics and an analysis of stability have been presented by Rosen and Seter [14,15], which employ blade element method with special attention to the effects of low Reynolds number. Further experimental work can be found in Refs. [4], [12], [14], [16], and [17].

It is the goal of this paper to provide a more compact blade element theory-based model for a samara in vertical descent by neglecting lateral movement and assuming negligible roll angle. This model, presented in Sec. 2, is explored for analysis of both steady-state, and transient behavior. Expanding upon the work in Ref. [18], an experimental setup and results are discussed in Sec. 3 for tuning and validation of simulation results. Furthermore, an analytical derivation of conservative boundaries for aerodynamic and physical properties that allow steady autorotation is discussed in Sec. 4, and a stability analysis of the autorotational equilibrium is conducted in Sec. 5. This is followed by Concluding Remarks, Acknowledgments, References, and an Appendix, in that order. The paper provides a model that balances accuracy and ease of implementation for stability analysis as well as design and control of simple single-winged rotorcraft similar to that of Refs. [9], [10], [19], and [20].

In Eq. (9), note that the limits of integration r₀ and r_f are not equal to 0 and R, respectively. This is because, as evident from Fig. 1(b), the wing-span that contributes to the aerodynamic forces starts from $r0>0$ and ends at an $rf<R$⁠. The bottom boundary of r₀ is taken to account for losses due to the rounded seed geometry near the center of mass, whereas the top boundary of r_f is taken to account for tip losses. For practicality, constant values have been selected for r₀ and r_f as ratios of a given samara's radius, R (see Table 3). Also, as the forces and moments are resolved, the following two simplifying assumptions are applied:

1. The effect of $dFx3$⁠, the elemental force along the blade span, is neglected. Thus, $Fx3≈0$ is assumed and its effect on aerodynamics is neglected. This simplifies the dynamic model. The assumption can be removed by using the net relative wind velocity on the (x₃, y₃) plane and accordingly considering blade elements to be tilted from the y₃ axis instead of being parallel to it, as shown in Fig. 2(b).

2. The rolling moment is assumed $Mx3≈0$⁠. This is based on the observation that the samara's motion is dominated by the yaw rate $ϕ˙$⁠, the pitching motion $θ,θ˙$ (coning), and the vertical motion v₀.

where f is a tunable factor to account for the nonuniform mass distribution which differs from samara to samara. The resulting reduced moment of inertia helps achieve coning angles that correlate well with data published in the literature, [8]. For the purposes of producing a simplified model, the center of mass is taken to be at O.

As mentioned in Sec. 2.1, the integration of independent blade elements perpendicular to the span ignores the spanwise airflow of the leading edge vortex which has been documented for single-winged samara structures [11–13]. This simplification is emphasized by assumption 1 (above) which states that $Fx3≈0$⁠. It will be shown in Sec. 3 that a strong agreement with steady-state values can be achieved with corrective adjustment of two tunable parameters: f from Eq. (10) and later $CD0$ from Eq. (19). The goal of this model is to provide a platform for reaching accurate steady-state results with realistic dynamic behavior. It is expected that a direct implementation of leading-edge vortex aerodynamics would improve the predictive accuracy of the model for the extreme dynamics before settling into autorotation. This addition would, however, significantly increase model complexity which can be detrimental to control applications.

Equation (20) must be solved numerically to determine the steady value of λ. Once solved, θ can be determined from either Eqs. (16) or (17). Here, θ can be solved explicitly. Next $ϕ˙$ can be determined explicitly from Eq. (18). Finally, knowing λ and $ϕ˙$⁠, v₀ can also be explicitly determined from the definition of λ in Eq. (14).

For observation of samara performance and collection of validation data, an experimental setup was created in the following manner:

1. Samaras were collected from local maple trees (specifically Red Maples—Acer rubrum).

2. A measuring tape was hung plum in camera view of a Samsung Z Fold3.

3. Samaras were dropped from sufficient height in plane with the measuring tape and recorded at 960 frames per second (visualized in Fig. 3).

4. Samaras were also recorded from a top-down perspective at the same frame rate.

5. The collected footage was analyzed frame-by-frame to extract θ, $ϕ˙$⁠, and v₀.

It is important to note that for accurate falling speed estimation, samaras should be recorded falling in the same plane as the measuring tape. Analyzing recordings of samaras with significant offset from the distance reference requires correction for parallax effects and can increase error. Furthermore, angle measurements should be taken when the samara is nearest to the center of view, such that minimal perspective distortion is present. A visualization of angle measurement is seen in Fig. 4. A sample video of samara descent from both front and top views is available at the following link.¹

Performance measurements were conducted for five sample samaras. The physical parameters of this sample group are presented in Table 1, along with the physical parameters of three larger samaras from previous works [8,12]. The parameters f_opt and $CD0opt$ are optimized values and will be discussed in Sec. 3.2. It should be noted that the Norberg and Holden 2 specimens are not Red Maple samaras, but Norway maple (Acer plantanoides) and silver maple (Acer saccharinum), respectively. Experimental values for steady-state θ, $ϕ˙$⁠, and v₀ are presented in Table 2 for the five tested samaras as well as the three samaras from the literature.

It has been observed that samaras, even within one species, vary significantly in size, mass, curvature, etc. For this reason, selecting an appropriate value for $Iy3y3$ and C_D is a nontrivial effort. Direct measurement of these values requires microscale testing and in many cases, destruction of a specimen for analysis of mass distribution. For this reason, two tunable parameters have been added, namely, f and $CD0$ in Eqs. (10) and (19), respectively, to adjust approximations of $Iy3y3$ and C_D for variation in morphology and higher-order aerodynamic response. Given the model above and collected experimental data, as well as other important modeling parameters listed in Table 3, a two-dimensional steepest descent optimization can be employed to select values for f and $CD0$⁠. A successful optimization should not only produce modeled steady-state values of θ, $ϕ˙$⁠, and v₀ which show good agreement with experimental values but also produce realistic values for f and $CD0$⁠.

Optimized values for the eight studied samaras are presented in Table 1. It can be see in Table 2 that strong agreement is achieved for θ, $ϕ˙$⁠, and v₀. For a study of appropriate $Iy3y3$ values, samara 2 was separated into two pieces: seed and wing. These separate sections were measured for mass and physical dimensions, and an approximation of $Iy3y3$ was performed by considering the seed to be a cylinder and the wing a rectangular flat plate, each of uniform density. From this analysis, it was found that $Iy3y3≈2×10−9kg m2$ for samara 2. Due to the rounded geometry of both the seed and wing, the above approximation is an overestimate. Given the measured $Iy3y3$⁠, Eq. (10) can be used to show $f≈0.47$⁠. This value is larger, yet of similar magnitude to the value optimized through steepest descent. The optimized values of $CD0$ in Table 1 are also realistic, as they are comparable to the low Reynolds Number performance of insect wings [14]. As will be discussed further in Sec. 4, the optimized $CD0$ values also fall within the necessary range of stability.

As previously stated, the above steady-state model was produced for an arbitrary roll angle, ψ. This allows for simulation of a range of equilibrium results for a span of ψ values. To explore the effects of varying ψ, the physical parameters of the Norberg specimen (see Table 1) were inputted, and the steady-state values of $ϕ˙$⁠, θ, v₀, and λ were modeled for a range of $−15deg<ψ<15deg$⁠. The results of this analysis are presented in Figs. 5 and 6.

It can be observed in Fig. 5(a) that a significant error in steady-state $ϕ˙$ is associated with negative ψ. Furthermore, the steep, positive correlation of θ and ψ seen in Fig. 5(b) suggests that adjusting ψ by a few degrees to either side of ψ = 0 results in a significant deviation of steady-state θ agreement. These trends, however, are produced using the f and $CD0$ values of Table 1, which were optimized under the assumption of ψ = 0. For this reason, an analogous three-dimensional steepest descent algorithm has been employed to not only optimize f and $CD0$⁠, but to also allow ψ to vary. For this optimization, the values of Table 1 were taken as initial guesses to see if significant deviation would occur when the assumption of ψ = 0 was removed. The results of this process are presented in Table 4. The resulting optimized ψ values do not deviate further than 0.6 deg from 0 deg and produce only minimal improvements in steady-state agreement of θ, $ϕ˙$⁠, and v₀. Furthermore, the variability in magnitude and sign of the optimized ψ suggests that such a correction is specific to individual morphological differences among samaras rather than suggesting a innate, nonzero roll for all samaras.

The transient behavior of a falling samara is next studied through simulations. Steady-state analysis has shown that it is reasonable to assume ψ = 0 near equilibrium. It has been observed that a falling samara establishes autorotation in a matter of a few seconds. This coupled with the minuscule rolling moment of inertia, $Ix3x3$⁠, suggests that achieving negligible roll angle will be nearly instantaneous. It is therefore reasonable for this control-oriented model to extend the assumption of ψ = 0 to transients. For consistency with the steady-state analysis, ψ has been presented as a constant input. The resulting dynamic model is given in Eqs. (11) and (12). Transient simulation results for θ and $θ˙$ as well as $ϕ˙$ and v₀ are presented in Figs. 7 and 8, respectively. Case 1 (bold) shows the model response with all initial conditions set to 0 (i.e., ${θ θ˙ ϕ˙ v0}={0 0 0 0}$⁠). Case 2 (dashed) has initial conditions of $θ=45 deg, θ˙=0.175$ rad/s, $ϕ˙=4$ rev/s, and $v0=−0.4$ m/s.

For the dynamic responses of θ, $θ˙, ϕ˙$⁠, and v₀, an aggressive stabilization is observed before smoothing out at approximately 0.5 s. The following smooth dynamic settles into steady autorotation within the first second. Special attention should be given to the response of v₀ which at approximately 0.25 s achieves maximum descent speed. Past this point, the samara slows to its autorotative terminal velocity. The steady values of the dynamic simulation show complete agreement with those of the steady-state model, as expected.

The equilibrium of autorotation can be viewed as a balancing of counteracting yawing moments along the span of a rotor-craft blade. This effect is visualized in Fig. 9 which displays the moment provided by each blade element of a samara. For the presented case (case 1), the entire blade provides a moment to increase the rotation of the samara. This effect reduces over time until steady-state is reached. It can be seen that at 1 s an approximately equal region of positive and negative moment is present, suggesting equilibrium. This result agrees with the predicted performance of helicopter blades in autorotation [22].

Observation of samaras found in nature and simulation of the presented model from various initial conditions has suggested the autorotational equilibrium of a single-winged samara is characterized by a large region of attraction. An investigation of samara stability will be presented in Secs. 4 and 5.

From the above transformed necessary condition, which is valid for the case when $ψ=0 deg$⁠, a solution for $(λ/χ)$ can exist in Eq. (27) only if $γ>0$⁠, implying $CD0>0$⁠. The term $CD0$ in the drag coefficient of Eq. (19) is thus a critical component in determining autorotational equilibrium. Section 5 presents a stability analysis of the autorotational equilibrium.

In $Ae*$⁠, the dynamics of θ and $θ˙$ are decoupled from that of $ϕ˙$ and v₀, thus resulting in two $(2×2)$ linear systems. In $A¯e$⁠, the dynamics of θ, $θ˙$ and v₀ are decoupled from that of $ϕ˙$⁠, resulting in one $(3×3)$ and a scalar dynamical system. In both cases, the stability analysis is simplified. To evaluate which simplification is better suited for analysis, a numerical root-locus study is performed. The four main parameters of the samara model, namely, the mass m, the blade radius R, the mass distribution factor f, and the drag factor $CD0$ are varied by $±30%$ and the corresponding eigenvalues are plotted for the original linearized system $Ae$ and the simplified ones, $Ae*$ and $A¯e$⁠. The eigenvalue variations and comparisons are shown in Fig. 10. It is noted that for all parameter perturbations in Fig. 10, the eigenvalues of $Ae$⁠, which is the original linearized system, are well-approximated by those of $Ae*$⁠. Furthermore, they are better approximated by the eigenvalues of $Ae*$ than those of $A¯e$⁠. The same observation is made for decoupling the dynamics of the descent velocity, v₀, from that of θ, $θ˙$⁠, and $ϕ˙$⁠. The decoupling in the form of two $(2×2)$ linear systems, as in $Ae*$⁠, thus provides a way to do a tractable stability analysis of the original linearized system while better retaining the characteristics of $Ae$ than $A¯e$⁠. It is also evident from Fig. 10 that the off-block-diagonal entries of $Ae$⁠, with respect to the block-diagonal form of $Ae*$⁠, do not produce enough effect on its eigenvalues to induce instability.