## Abstract

This paper presents a control-oriented model for describing the steady-state and dynamic behavior of a single-winged samara seed-pod in autorotative descent. A negligible lateral center of mass motion and constant, prescribed roll-angle to develop a simplified and compact model. Spanwise aerodynamic dependence is exchanged for an independent blade element representation with two tuned parameters to account for the effects of leading-edge vortex phenomena. The resulting model is a fourth-order nonlinear dynamical system. The accuracy of this model is established by validating it against our own experimental data as well as against those reported in the literature by other researchers. The validation exercise reveals that zero roll-angle is a viable assumption that significantly reduces model complexity while retaining accuracy. A necessary condition is derived for the existence of steady autorotation of the samara under free descent. Furthermore, a stability analysis is conducted suggesting that the eigenvalues of the fourth-order system, linearized about the autorotational equilibrium, can be well-represented by those of two decoupled two-dimensional systems. The analysis reveals the critical parameters that determine stability of sustained autorotation. Such stability analysis provides a platform for similar analytical exploration of future model improvements. The validity of this compact model suggests the plausibility of designing and controlling simple autorotative mechanisms based on these dynamics.

## 1 Introduction

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].

## 2 Modeling

### 2.1 Dynamic Model.

*x*-axis, and $Fx3$ refers to the net force in the direction of the body-fixed x-axis. The same naming convention follows for the body-fixed

*y*- and

*z*-axes. The angular velocity of the samara in the body-fixed frame $(x3,y3,z3)$ can be expressed using Euler angles as follows:

*P*shown in Fig. 1(a) can be expressed relative to

*O*as given below:

*O*of the samara has motion predominantly in the vertical direction with negligible motion in lateral directions. From Eqs. (3) and (4)

*P*is assumed $vw/P=\u2212vP$. This approximation neglects the effects of the leading edge vortex of a single-winged samara which couples blade elements. Efforts were made in Ref. [14] to account for spanwise flow with skewed blade elements. Instead, a corrective tuning adjusts the coefficient of drag and moment of inertia for the effects of the coupled aerodynamics. This correction will be discussed further in Secs. 2.2 and 3.2. From Fig. 2 for an element along the samara blade, the local drag and lift forces,

*dD*and

*dL*, respectively, are given by

*α*is the spanwise local angle of attack and from Eq. (5)

*y*

_{3}and

*z*

_{3}directions are

### 2.2 Simplifying Assumptions.

In Eq. (9), note that the limits of integration *r*_{0} 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*

_{0}is taken to account for losses due to the rounded seed geometry near the center of mass, whereas the top boundary of

*r*is taken to account for tip losses. For practicality, constant values have been selected for

_{f}*r*

_{0}and

*r*as ratios of a given samara's radius,

_{f}*R*(see Table 3). Also, as the forces and moments are resolved, the following two simplifying assumptions are applied:

The effect of $dFx3$, the elemental force along the blade span, is neglected. Thus, $Fx3\u22480$ 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*_{3},*y*_{3}) plane and accordingly considering blade elements to be tilted from the*y*_{3}axis instead of being parallel to it, as shown in Fig. 2(b).The rolling moment is assumed $Mx3\u22480$. This is based on the observation that the samara's motion is dominated by the yaw rate $\varphi \u02d9$, the pitching motion $\theta ,\theta \u02d9$ (coning), and the vertical motion

*v*_{0}.

*O*(see Figs. 1(a) and 1(b)). It is imperative that this be taken into consideration when approximating the value for $Iy3y3$. For this paper, Eq. (10) was employed for moment of inertia approximation

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\u22480$. 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.

*ψ*, remains small. With this observation, and consideration of the symmetry of the samara's airfoil profile, it is reasonable to assume that the angle

*ψ*will be nearly if not exactly zero at steady-state. Further, the assumption of $Mx3\u22480$ along with $Ix3x3=0$ and $Iz3z3=Iy3y3$ implies that the first equation of Eq. (1) is identically zero. The dynamics and statics of the system can be studied for different constant values of

*ψ*. Statics analysis reveals that only a small range of

*ψ*around zero is allowable. This simplification reduces the statics problem to that of determining three unknowns, namely, $\varphi \u02d9$,

*θ*, and

*v*

_{0}. The dynamics reduce to four states, namely, $\varphi \u02d9,\u2009\theta \u02d9$,

*θ*, and

*v*

_{0}. The simplified dynamics of an autorotating samara are next provided. From the second and third equations of Eq. (1), from Eq. (2), and imposing $\psi \u02d9=0$ on $\omega x3,\u2009\omega y3,\u2009\omega z3$ in Eq. (3)

### 2.3 Conditions For Steady Autorotation.

*ψ*, is allowed to assume nonzero constant values. The value,

*ψ*, is treated as an input in determining the possible set of autorotational equilibria. The exercise confirms that feasible steady values of

*ψ*lie only a few degrees around

*ψ*= 0. The equations for steady autorotation, obtained from Eqs. (1)–(3) are

*χ*represents the span-ratio which is the position of a blade element with respect to the length of the blade. The ratio

*λ*represents the tip-speed-ratio which is the relation of the vertical descent speed to the speed of the tip of the blade. The ratio $\lambda /\chi $ produces a local speed ratio which describes the relation of the vertical descent speed of the entire samara to the local velocity of a blade element. Inclusion of the local speed ratio provides insight into spanwise characteristics. Applying the ratios in Eq. (14) produce the following expressions for the local $||U\u221e||$ and

*α*under steady autorotation:

*ψ*, the right-hand side of Eq. (15) reduces to the local speed ratio of the blade element. Expanding Eq. (13) using Eqs. (8), (9), and (15) produces the following system of equations:

*C*and

_{L}*C*, the coefficients of lift and drag, respectively [21]. Specifically

_{D}*θ*,

*v*

_{0}, and

*λ*for a range of values of

*ψ*. The process of determining the steady conditions is to first solve for

*λ*from Eqs. (16) and (17). This is done by noting that

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 $\varphi \u02d9$ can be determined explicitly from Eq. (18). Finally, knowing *λ* and $\varphi \u02d9$, *v*_{0} can also be explicitly determined from the definition of *λ* in Eq. (14).

## 3 Results and Validation

### 3.1 Experimental Setup and Results.

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

Samaras were collected from local maple trees (specifically Red Maples—

*Acer rubrum*).A measuring tape was hung plum in camera view of a Samsung Z Fold3.

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

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

The collected footage was analyzed frame-by-frame to extract

*θ*, $\varphi \u02d9$, and*v*_{0}.

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.^{1}

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 *θ*, $\varphi \u02d9$, and *v*_{0} are presented in Table 2 for the five tested samaras as well as the three samaras from the literature.

### 3.2 Parameter Optimization.

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*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

_{D}*f*and $CD0$. A successful optimization should not only produce modeled steady-state values of

*θ*, $\varphi \u02d9$, and

*v*

_{0}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 *θ*, $\varphi \u02d9$, and *v*_{0}. 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\u22482\xd710\u22129kg\u2009m2$ 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\u22480.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.

### 3.3 Analysis of Variable Roll.

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 $\varphi \u02d9$, *θ*, *v*_{0}, and *λ* were modeled for a range of $\u221215deg<\psi <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 $\varphi \u02d9$ 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 *θ*, $\varphi \u02d9$, and *v*_{0}. 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.

### 3.4 Dynamics Simulations.

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 $\theta \u02d9$ as well as $\varphi \u02d9$ and *v*_{0} are presented in Figs. 7 and 8, respectively. Case 1 (bold) shows the model response with all initial conditions set to 0 (i.e., ${\theta \u2009\theta \u02d9\u2009\varphi \u02d9\u2009v0}={0\u20090\u20090\u20090}$). Case 2 (dashed) has initial conditions of $\theta =45\u2009deg,\u2009\theta \u02d9=0.175$ rad/s, $\varphi \u02d9=4$ rev/s, and $v0=\u22120.4$ m/s.

For the dynamic responses of *θ*, $\theta \u02d9,\u2009\varphi \u02d9$, and *v*_{0}, 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*_{0} 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.

## 4 A Necessary Condition for the Existence of Autorotational Equilibrium

*θ*, and

*v*

_{0}are constant. Note that the solution of steady autorotation can only exist if there exists a solution

*λ*for Eq. (20), as

*r*varies over the span $r\u2208[r0,rf]$. Substituting the expressions of $My3$ and $Mz3$ from Eq. (16) and (17), respectively, Eq. (20) can be rewritten as

From the above transformed necessary condition, which is valid for the case when $\psi =0\u2009deg$, a solution for $(\lambda /\chi )$ can exist in Eq. (27) only if $\gamma >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.

## 5 Stability Analysis

*i*=

*2, 3, 4 and $j=1,2,3,4$ are given in the Appendix. The characteristic equation of $Ae$ is*

In $Ae*$, the dynamics of *θ* and $\theta \u02d9$ are decoupled from that of $\varphi \u02d9$ and *v*_{0}, thus resulting in two $(2\xd72)$ linear systems. In $A\xafe$, the dynamics of *θ*, $\theta \u02d9$ and *v*_{0} are decoupled from that of $\varphi \u02d9$, resulting in one $(3\xd73)$ 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 $\xb130%$ and the corresponding eigenvalues are plotted for the original linearized system $Ae$ and the simplified ones, $Ae*$ and $A\xafe$. 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\xafe$. The same observation is made for decoupling the dynamics of the descent velocity, *v*_{0}, from that of *θ*, $\theta \u02d9$, and $\varphi \u02d9$. The decoupling in the form of two $(2\xd72)$ 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\xafe$. 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.

*A*

_{21}given in Eq. (A1). All terms of Eq. (A1) will be negative if $v0,e<0,\u2009\varphi \u02d9e>0,\u20090<\theta e<\pi /4,\u2009CD0>0,\u2009Iy3y3>0$ and if all integral terms in Eq. (A1) are positive. Note that the samara geometry, aerodynamic properties, and equilibrium conditions ensure that the above requirements are satisfied and thereby $A21<0$ provided $0<\theta e<\pi /4$. Next consider the expression of

*A*

_{22}given in Eq. (A2). It can be expressed as

*A*

_{22}is also guaranteed by the samara geometry, aerodynamic properties, and equilibrium conditions. This establishes the stability of the subsystem $Ae,1*$. Next, consider the subsystem $Ae,2*$, given by

*A*

_{33}can be expressed as

*A*

_{44}can be expressed as

*A*

_{43}can be expressed as

*A*

_{34}can be expressed as

The condition in Eq. (44) is not restrictive. This can be verified by noting that since a typical value of $\lambda e\u22480.25$. Substituting this value into Eq. (44) yields $CD0\u22640.93$. For a steady-state $\lambda \u22480.28$, as obtained in Fig. 6(b) for $\psi =0$ deg, the requirement becomes $CD0\u22641.14$. Thus, stability is also dependent on the parameter $CD0$ in our model and the stability condition on $CD0$ is more stringent in comparison to the necessary condition of autorotation, as derived in Sec. 4. From the zoomed-in views of Figs. 10(a) and 10(d), it is evident that the dominant eigenvalue of $Ae$ is around –9. This would imply a rise-time [23], of about $1.8/9=0.2$ s. This matches well with the approximate rise-times in the transient simulations of Figs. 8(a) and 8(b).

It should be noted that, while there are alternate methods to study the stability of the aforementioned system, such as using the *Gerschgorin circle theorem* [24] or its block analog [25], or even the *small gain theorem* for interconnected systems [26], these analysis yield significantly more conservative stability criteria for this study. It is also arguable that the partial derivatives in Eq. (29) could be numerically determined instead of the analytical approach adopted in this paper. Indeed, the numerical approach was attempted and it was observed these calculations were lacking in convergence. A similar convergence-related observation was made in Ref. [15], which ultimately relied on numerical values based on additional investigations on the accuracy of the numerical process.

## 6 Conclusions

A simplified and compact model for the steady-state behavior of a single-winged samara has been presented. It has been shown that an assumption of negligible roll angle, *ψ* = 0, is reasonable, and can simulate accurate steady-state and realistic transient behavior. With small, nonzero *ψ* values, improvements can be made to simulation agreement with data, but these deviations appear to be more correlated with individual morphological differences between samaras as opposed to a general trend of *ψ* deviation. Experimental results have been shown for five Red Maple samaras to tune and validate the model. The use of two tunable parameters has adequately encapsulated the effects of complex geometry and higher-order aerodynamic phenomena on drag and moment of inertia. Necessary conditions for the existence of an autorotation equilibrium have been analytically derived with further analysis on stability, suggesting ranges of parameters that can be employed in the development of biomimicking aerodynamic mechanisms. The stability analysis reveals that in the neighborhood of the autorotational equilibrium, the eigenvalues of the higher-order coupled system can be well-approximated by two decoupled lower-order systems. This lends itself to a tractable stability analysis which shows the reliance of stability on key parameters of the samara. Furthermore, the compact form of the presented model lends itself well to implementation in aerodynamic control of single-winged crafts. A topic of future interest is the exploration of lateral stability and the effects of wind gusts as well as employing the presented model in a control system.

## Funding Data

National Science Foundation CMMI, USA (Grant No. 1762986; Funder ID: 10.13039/100000001).

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*C*=_{D}coefficient of drag

- $CD0$ =
additional drag effects

*C*=_{L}coefficient of lift

*f*=adjustment for mass distribution

*F*=net force

*g*=acceleration due to gravity

*I*=moment of inertia

*m*=mass

*M*=net moment

*r*=radial position

*R*=blade radius

*v*_{0}=vertical velocity of samara

*w*=blade width

*α*=local angle of attack

*θ*=pitch angle

*λ*=tip speed ratio

*ρ*=air density

- $\varphi $ =
yaw angle

*χ*=span ratio

*ψ*=roll angle

*ω*=angular velocity