Investigation of the Free-Fall Dynamic Behavior of a Rectangular Wing with Variable Center of Mass Location and Variable Moment of Inertia

Abstract

In recent years, the air-drop launch technology of near-space UAVs has attracted much attention. Between downfall from the carrier and the flight control system’s initiation, the UAV presents free-fall movement. This free-fall process is very important for the control effect of the flight control system and is also crucial for the safety of the UAV and the carrier. Focus is required on two important dynamic parameters of the UAV: the moment of inertia and the center of mass position. In this paper, we used a quasi-steady model proposed by predecessors to address the flat-plate falling problem with modifications to describe the freely falling motion of the wing. Computational fluid dynamics (CFD) were used to simulate the free-fall movement of the wing with various parameters, and the wing release behavior was analyzed to check the quasi-steady model. Research shows that the movement characteristics of the falling wing are mostly reflected in the longitudinal plane, and the developed quasi-steady analytical model can more accurately describe the dynamic behavior of free-fall to some extent. By using CFD methods, we further investigated the aerodynamic performance of the free-fall wing. The results show that the wing mainly presents tumbling and fluttering motion. Changing the moment of inertia around the tumbling axis changes the tumbling frequency and the time point as the wing enters tumbling. In contrast, changing the position of the center of mass significantly changes the form of falling and makes the free-fall motion more complex. Therefore, it is necessary to carefully configure the center of mass in the UAV design process.

1. Introduction

In recent years, the air-drop launch technology of unmanned aerial vehicles (UAVs) has become a new hot topic [1,2,3,4,5]. It refers to the use of balloons, airships, rockets, re-entry spacecraft, and other vehicles to release the UAVs at a specified altitude and speed. Compared to the traditional take-off mode, air-drop launch is independent of the runway and ground environment. Moreover, it can save time and costs for some special UAVs to reach the designated airspace and improve safety in their approach. For example, a near-space solar-powered UAV [6,7,8] requires a lot of time and energy to climb to its cruising altitude, and it always encounters safety problems caused by atmospheric turbulence in the climb process. If air-drop launch technology is adopted, the problems associated with near-space solar-powered UAV takeoff sites can be effectively solved, flight safety can be improved, and limited energy can be used for the mission flight rather than for climbing. Hence, air-drop launch technology offers broad application opportunities in the fields of high-altitude aircraft, fast-deploying UAVs, and alien exploration [9]. Although much research on air-drop launch technology has been conducted in recent decades [3,10,11,12], there are still open questions that need to be understood if further improvements are to be realized.

For near-space solar-powered UAVs, the ideal air-drop launch process is undoubtedly free-fall motion from balloons [13,14] without any energy injection. However, during the uncontrolled falling motion between being dropped from the carrier to the start of flight control, the UAV presents uncontrolled falling. It is well known that the free-fall motion exists widely in nature and industry, such as the falling motion of leaves, paper, feathers, dust particles, diffusion of plant seeds [15,16,17], the uncontrolled re-entry of spacecraft [18], etc. Lightweight objects such as falling leaves, pieces of paper, thin disks [19], and plant seeds exhibit complex and abundant dynamic behaviors in the falling process. This is mainly due to the strong coupling between aerodynamics and kinematics, which makes it difficult to predict fall trajectory and to analyze the unsteady aerodynamics of the object. As early as 1854, Maxwell studied the irregular tumbling phenomenon of thin plates [20]. Willmarth et al. [21] generated a phase diagram of dropping disks with steady descent, tumbling, and fluttering. EH Smith [22] and Field [23] observed the chaotic motion of thin disks in experiments and provided phase diagrams of the motion form with respect to the dimensionless moment of inertia and Reynolds number. With the development of numerical simulation methods and experimental technology, Andersen [24] used high-speed photography to record the the falling track of a flat plate and compared experimental results with the numerical solution of the two-dimensional Navier–Stokes equation. Lam et al. [25] studied the free-fall of a disk in a stably stratified fluid using a similar experimental method. Esteban [26] used high-speed imaging technology to investigate how roundness influences the fall trajectory. Toupoint [27] experimentally studied the falling motion of a finite-length cylinder in a low-viscosity fluid. Chern et al. [28] used the direct-forcing immersed boundary (DFIB) method to numerically study freely dropped objects and compared the falling behavior of discs, equilateral triangles, and oval objects. Kim [29] used a three-dimensional particle tracking velocimeter (PTV) and a particle image velocimeter (PIV) to reveal the induced flow for a cone with different centers of gravity and density ratios during free fall. In addition, Tam [30] developed a theoretical model for a tumbling flexible wing and observed the dynamics of a passive flexible wing. Similar studies were conducted by Esteban [31,32], Toby Howison [33], Lee [34], and Zhou [35]. Several dynamic models have been proposed for freely falling plates and disks, such as those proposed by Kozlov [36]; Tanabe and Kaneko [37]; Belmonte et al. [38]; and Andersen et al. [39].

Although the aforementioned studies provided valuable understanding of the free-fall dynamic behavior of different kinds of objects, the aerodynamic mechanisms have not been fully discovered. The majority of the literature deals with small, light objects falling in a laboratory environment with a rather low Reynolds number, and very few studies in the literature have provided any explanations for the free-fall motion of large objects, such as wings or UAVs, into large-scale space. Hence, it is necessary to carry out more extensive investigations on the free-fall dynamic behavior of large wings and to deeply study the aerodynamic influences of several important parameters on the free-fall motion of air-drop-launched UAVs. This is the topic of the present work.

Considering that the position of the moment of inertia (MOI) and the center of mass (COM) play important roles in determining the dynamic behavior of falling disks and cones [25,29], we studied the free-fall dynamic behavior of a rectangular wing with a variable center of mass location and a variable moment of inertia, mainly focusing on the longitudinal movement characteristics of the falling wing. The following sections describe the development of a quasi-steady analytical model based on the Andersen–Pesavento–Wang model with some corrections. Furthermore, in order to verify the feasibility of the analytical model and to improve the understanding of a freely falling wing, a two-dimensional approximate computational fluid dynamics (CFD) method and a set of wing free-fall experiments were also developed and are introduced in detail. Then, the dynamic process of a freely falling wing was simulated and analyzed using the analytical model and experiments. Simulations were conducted under different conditions by adjusting the parameters and repeating the experiment several times in different states, and the effects of the wing’s MOI and COM position on free-fall were studied.

2. Research Methods and Validations

The COM method has the advantage of being simple in form and clear in physical implications and has been applied to the analysis of the dynamics of the free-fall of an object, but has some limitations, such as not taking into account the uneven mass distribution. This paper therefore developed an analytical model focused on this problem. In order to verify the feasibility of the quasi-steady analytical model and to gain a clearer understanding of the free-fall motion, a quasi-steady approximate 2D CFD numerical simulation method based on overset meshes for airfoils was proposed to enhance the understanding and knowledge of the development of aerodynamic characteristics through numerical simulation of the free-fall of elliptical airfoils. CFD free-fall simulations were carried out using a wing with elliptical airfoil and the chord length was 0.25 m as an example. Additionally, a 3D realistic wing free-fall experimental test system was developed to realistically represent the dynamic behavior of the freely falling wing and to validate the results with the analytical model of this paper further to enhance the understanding based on the theoretical analysis.

2.1. Quasi-Steady Analytical Model

The scenario was based on the Northwestern Polytechnical University’s “MY” full-wing solar-powered UAV [7,40]. The simplified wing is shown in Figure 1a. Eccentricity of the airfoil: $e=b/a=1/10$, where b and a are the semiminor axis and the semimajor axis of the ellipse respectively. Past research about freely falling plates has shown that the movement characteristics in the $XOZ$ plane are the most obvious, where the $XOZ$ plane is an aircraft-carried normal earth-fixed system. The aspect ratio of the wing in Figure 1 is 4, so this study simplified the three-dimensional falling problem of the wing to the motion in the $XOZ$ plane. The force model is based on the body coordinate system $X^{\prime }O{Z}^{\prime }$. $X_C$ is the chordal distance between centroid and geometric center, and the pitch angle is $\theta$, see Figure 1b.

Anderson, Pesavento, and Wang [39] established a finite dimensional pattern for the free-fall behavior of a flat plate and elliptical disk. In the model, the quasi-steady assumption is adopted for aerodynamic force, and the model is written as follows:

$$\left(m+m_{x^{\prime }}\right){\dot{v}}_{x^{\prime }}=(m+m_z)v_{z^{\prime }}\dot{\theta }-{\rho }_f\mathsf{\Gamma }{v}_{z^{\prime }}-mg\left(1-{\rho }_f{{\rho }_s}^{-1}\right)sin\theta -F_{x^{\prime }}^v$$

(1)

$$\left(m+m_{z^{\prime }}\right){\dot{v}}_{z^{\prime }}=-(m+m_{x^{\prime }})v_{x^{\prime }}\dot{\theta }+{\rho }_f\mathsf{\Gamma }{v}_{z^{\prime }}-mg\left(1-{\rho }_f{{\rho }_s}^{-1}\right)cos\theta -F_{z^{\prime }}^v$$

(2)

$$\left(I+J\right)\ddot{\theta }=(m_{x^{\prime }}-m_{z^{\prime }})v_{z^{\prime }}{v}_{x^{\prime }}-{\tau }_{\theta }+l_{\tau }{\rho }_f\mathsf{\Gamma }\sqrt{v_{x^{\prime }}^2+v_{z^{\prime }}^2}.$$

(3)

The velocity components $v_{x^{\prime }}$ and $v_{z^{\prime }}$, respectively, represent the speed in the $X^{\prime }$ and $Z^{\prime }$ directions of the body coordinate system. The variables m, I, ${\rho }_f$, and ${\rho }_s$ represent the mass of the wing, MOI about the $Y^{\prime }$-axis, density of the air and the wing, respectively. The added mass coefficients $m_{x^{\prime }}$ and $m_z^{\prime }$, and the added moment of inertia J in inviscid theory are coefficients given by Sedov [41]:

$$m_{x^{\prime }}=\pi {\rho }_{\theta }{b}^2,m_{z^{\prime }}=\pi {\rho }_{\theta }{a}^2,J=\frac{1}{8}\pi {\rho }_f{\left(a^2-b^2\right)}^2.$$

(4)

In the analytical model, the aerodynamic force of the freely falling wing is composed of an additional mass term, circulation term, and translational drag term. The circulation includes the translation term and rotation term, which can be written as

$$\begin{array}{ccc}\mathsf{\Gamma }=-\frac{1}{2}{C}_{T}V{c}_{A}lsin\left(2\alpha \right)+\frac{1}{2}{C}_R{c_A}^{2}l\dot{\theta },& c_A=2a,& V=\sqrt{v_{x^{\prime }}^2+v_{z^{\prime }}^2},\end{array}$$

(5)

where $\alpha$ is the angle of attack and the translation terms $F_{x^{\prime }}^v$, $F_{z^{\prime }}^v$ are given by:

$$\left(\begin{array}{c}{F}_{{}_{x^{\prime }}}^v\\ F_{z^{\prime }}^v\end{array}\right)=\frac{1}{2}{\rho }_f{V}^2{c}_{A}l{C}_D\left(\begin{array}{c}cos\alpha \\ sin\alpha \end{array}\right).$$

(6)

In the study of Wang et al. [42], the drag coefficient is given by:

$$C_D\approx C_A-C_{B}cos2\alpha .$$

(7)

Thus, the translational drag term is written as:

$$\left(\begin{array}{c}{F}_{x^{\prime }}^v\\ F_{z^{\prime }}^v\end{array}\right)={\rho }_{f}a\left(C_A-C_B\frac{v_{x^{\prime }}^2-v_{z^{\prime }}^2}{v_{x^{\prime }}^2+v_{z^{\prime }}^2}\right)\sqrt{v_{x^{\prime }}^2+v_{z^{\prime }}^2}\left(\begin{array}{c}{v}_{x^{\prime }}\\ v_{z^{\prime }}\end{array}\right).$$

(8)

In the model, dimensionless coefficients $C_T$, $C_R$, $C_A$, $C_B$ are related to the thickness and geometry of the airfoil [43]. Combining the finite dimensional model given by other researchers, the coefficient is determined as:

$$C_T=1.2,C_R=\pi ,C_A=1.92,C_B=1.55.$$

(9)

In angular movement, ${\tau }_{\theta }$ represents dissipative torque, caused by the wing’s angular movement, given by

$${\tau }_{\theta }=\frac{1}{2}{\rho }_f{C}_D\left(\frac{\pi }{2}\right)\underset{-\frac{c_A}{2}}{\overset{\frac{c_A}{2}}{\int }}\left|v_n\right|v_{n}rdr.$$

(10)

In Equation (10), $C_D\left(\frac{\pi }{2}\right)$ represents the drag coefficient at an angle of attack of $\pi /2$, r is the length between local position and the COM, and $v_n$ represents local normal velocity, given by $v_n=v_{{}_{z^{\prime }}}+r\dot{\theta }$. Considering the influence of the position of the COM, substitute Equation (10) as:

$$\begin{array}{c}{\tau }_{\theta }=\frac{1}{2}{\rho }_f{C}_D\left(\frac{\pi }{2}\right)\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|v_{z^{\prime }}+r\dot{\theta }\right|(v_{z^{\prime }}+r\dot{\theta })rdr\hfill \\ \le \frac{1}{2}{\rho }_f{C}_D\left(\frac{\pi }{2}\right)\left[\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|v_{z^{\prime }}\right|(v_{z^{\prime }}+r\dot{\theta })rdr+\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|r\dot{\theta }\right|(v_{z^{\prime }}+r\dot{\theta })rdr\right].\hfill \end{array}$$

(11)

Let

$$A=\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|v_{z^{\prime }}\right|(v_{z^{\prime }}+r\dot{\theta })rdr,B=\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|r\dot{\theta }\right|(v_{z^{\prime }}+r\dot{\theta })rdr.$$

(12)

Further, B can be expressed as

$$B=\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|r\dot{\theta }\right|(v_{z^{\prime }}+r\dot{\theta })rdr=\underset{0}{\overset{a-X_C}{\int }}\left|\dot{\theta }\right|(v_{z^{\prime }}+r\dot{\theta })r^{2}dr-\underset{-a-X_C}{\overset{0}{\int }}\left|\dot{\theta }\right|(v_{z^{\prime }}+r\dot{\theta })r^{2}dr.$$

(13)

The integral of A can be obtained

$$A=\underset{-a-X_C}{\overset{a-X_C}{\int }}\left|v_{z^{\prime }}\right|(v_{z^{\prime }}+r\dot{\theta })rdr=-\left|v_{z^{\prime }}\right|v_{z^{\prime }}{X}_{C}a+\frac{2}{3}\left|v_{z^{\prime }}\right|\dot{\theta }\left(2a^3+6{X_C}^{2}a\right).$$

(14)

Similarly, the integral of B is

$$B=-\frac{1}{3}\left|\dot{\theta }\right|v_{z^{\prime }}\left(6a^2{X}_C+3{X_C}^3\right)+\frac{1}{4}\left|\dot{\theta }\right|\dot{\theta }\left(a^4+12a^2{X_C}^2+2{X_C}^4\right).$$

(15)

Since $X_C$ is a small quantity, the higher order terms about $X_C$ are ignored and only the primary term is retained. Considering that Equality (11) seriously overestimates the dissipative moment term, it is assumed that inequality is always valid after ignoring a small number of higher order terms. Thus, ${\tau }_{\theta }$ can be approximated as:

$${\tau }_{\theta }\approx \frac{1}{2}{\rho }_f{C}_D\left(\frac{\pi }{2}\right)C_S\left[-\left|v_{z^{\prime }}\right|v_{z^{\prime }}{X}_{C}a+\frac{4}{3}\left|v_{z^{\prime }}\right|\dot{\theta }{a}^3-2\left|\dot{\theta }\right|v_{z^{\prime }}{a}^2{X}_C+\frac{1}{4}\left|\dot{\theta }\right|\dot{\theta }{a}^4\right],$$

(16)

where $C_S$ is a scale factor to reduce the overestimation of ${\tau }_{\theta }$. Otherwise, since the airfoil is an ellipse, the parameter $l_{\tau }$ in Equation (3) equals zero [42]. Then, let the total aerodynamic force in the $Z^{\prime }$ direction be $F_{z^{\prime }}^v$. Finally, the angular dynamic equation can be written as:

$$\left(I+J\right)\ddot{\theta }=(m_x-m_{z^{\prime }})v_{z^{\prime }}{v}_{x^{\prime }}-{\tau }_{\theta }-F_{z^{\prime }}^S{X}_C.$$

(17)

So far, we have deduced the quasi-steady three-degree-of-freedom dynamic model for the freely falling wing. Given the initial state, the motion trajectory of the wing in the falling process can be obtained by using the Runge–Kutta method (refer to Appendix A.1 for details).

In order to verify the feasibility of the analytical model and to gain a clearer understanding of the free-fall motion, CFD free-fall simulations were carried out using a wing with elliptical airfoil and the chord length was 0.25 m as an example, while the MOI of the wing was $4.5\times {10}^{-4}$ kg·m${}^2$. The wing drop experiment was designed to compare with the CFD result and the analytical model.

2.2. CFD Numerical Method

Having derived the theoretical model, this paper used a CFD method to further investigate the free-fall motion of the wing, aiming to gain a deeper understanding of the aerodynamic phenomena involved in the process. Due to the significance of studying motion in the longitudinal plane, this paper utilized a two-dimensional CFD method to investigate the freely falling process in the longitudinal plane.

The commercially available CFD package Fluent version 19.0 was used to simulate the unsteady flow-field around the wing configuration in this study; the two-dimensional RANS governing equations were discretized by the standard cell-centered finite volume scheme. Air was considered as the ideal gas, and the Sutherland viscosity law was used for better accuracy in terms of viscous effects. The turbulence effects were considered by using the k-$\mathsf{\Omega }$ shear stress transport (SST) turbulence model. Moreover, the second-order accurate Roe flux-difference upwind scheme was used for the convective fluxes of turbulence equations, and the lower-upper symmetric Gauss–Seidel (LU-SGS) implicit method was adopted for the time integration. Figure 2 displays the computational mesh around the baseline airfoil. The O-topology grids with 361 cells in the normal-to-wall direction and 150 cells around the airfoil were used. The boundary layer thickness was applied based on the Reynolds number of the flow, and the height of the first layer was set to be $5.76\times {10}^{-5}$ m so that the non-dimensional y-plus fell in a range less than 0.5. An overset grid approach was used for calculation; the background grid had an x-direction range of −12 to 12 m and a z-direction range of −15 to 2 m. The overset component was fused with the background grid using a cell-based grid merging strategy. The simulation of free-fall was realized using the SDOF macro in UDF, with a time step of $1\times {10}^{-4}$ s. The grid of the calculation domain is shown in Figure 2 and Figure 3. Through a two-dimensional CFD calculation, we can obtain the flow field information of the freely falling wing.

2.3. Experiment Method of Freely Falling Wing

In this section, an experimental method was designed for the freely falling wing. The experiments on a real wing verify the accuracy of the theoretical quasi-steady analytical model, and demonstrate the conclusions drawn from the analytical model and CFD numerical simulation. The experiment needs to measure the attitude and position of the wing’s movement. In previous studies on small objects such as paper and discs, high-speed photography technology has been used to record the falling trajectory of objects [44], and the experiments on freely falling objects are mostly carried out in narrow spaces such as water tanks. This is because the size of objects such as paper and discs is small, and the space required for the full development of the falling trajectory is also small. Related to the vision and lighting conditions, in this context, high-speed photography technology is not reliable. In prior research, the dynamic characteristics of plates can be measured by sensors including accelerometers, rategyros and magnetometers [45]. Therefore, in this study, the inertial measurement unit (IMU), shown in Figure 4, installed at the center of mass was used to measure the acceleration, attitude angle, and angular velocity of the wing. After filtering and integration, it was possible to obtain the falling path.

In order to easily change the MOI of the wing about the $Y^{\prime }$ axis, a light thin stalk extends from the center of mass of the wing parallel to the chord. Metal clump weights with the same quality and shape are installed along the stalk, and the metal clump weights can be fixed at any position on the stalk. By making the clump weights and stalk symmetrical about the COM and by synchronously changing the position of the weights on the stalk, the MOI of the wing can be changed while keeping the mass and COM unchanged: see Figure 5.

However, it is very difficult to change the position of center of mass while keeping the moment of inertia unchanged. In order to change the position of COM, we adjusted the position of the weights on the stalk, but this will bring an obvious change of the inertia. Therefore, this work first researched the influence of the MOI on the freely falling wing, and then researched the effect of the COM position on falling motion. Finally, the total mass of the system was 253 g. The experimental procedure involved lifting the wing and release mechanism to a height, and releasing the wing free, as shown in Figure 6.

2.4. Validation and Discussion

To validate the quasi-steady analytical model and to improve the knowledge of the fluid field characteristics during the falling of the wing, this paper used CFD to numerically simulate the falling of a wing, with an initial attitude angle ${\theta }_0=-{75}^{\circ }$. The results of the analytical model and CFD simulation as well as the experimental results are shown in Figure 7.

It can be seen that the analytical model can predict the angular motion well during the falling process, and also gives an accurate indication of the aerodynamic forces on the wing.

For further investigation, this paper used the $\mathsf{\Omega }$-method (a new generation vorticity identification method; refer to Appendix A.2 for details) to visualize vorticity around the falling wing simulated by CFD numerical calculations shown in Figure 8.

The process of the tumbling of the wing can be divided into two parts: a “transition phase” and a “stable tumbling phase”. In the transition phase, the wing is accelerated from the initial state, and the angular motion of the wing is not obvious while accelerating, mainly in translation motion. As the speed increases, the wing starts to tumble, and the vortex on the surface of the wing starts to shed, forming the vortex street. The frequency of vortex shedding is low in the transition phase. As the wing enters the stable tumbling stage, the frequency of vortex shedding accelerates and finally remains relatively stable.

The trend of the flow field during the transition phase of the example is shown in Figure 9. The formation of the vortex is similar to the results in Wan’s research [46]. Since the airfoil is symmetrical, the leading edge of the ellipse is defined as a side of the long-axis of the ellipse which is upwind and, similarly, the trailing edge can be defined as that downwind. In the later stage of the accelerated downward phase from the initial position, along with the pitching motion (Figure 9a), a vortex around the leading edge gradually forms and gradually strengthens; meanwhile, the force in the Z direction of the wing increases. When the pitch angle increases, the flow separates from the leading edge, but the vortex is still attached to the leading edge (Figure 9b), and the force in the Z direction peaks. The leading edge vortex then begins to shed (Figure 9c,d) at the same time as $F_z$ decreases, and an opposite separation vortex also appears at the trailing edge (Figure 9e). Due to the wing’s rotation, the original trailing edge became the current leading edge, and the circulation around the leading edge increases again (Figure 9e–h). As the wing maintains the rolling motion, the flow separates again at the leading edge with the separate vortex shedding out, similar to the process in Figure 9a–c. This tumbling motion is repeated for the subsequent transition phases of the wing.

In Figure 7, it can be observed that, during the early stages of the transition phase, the CFD simulation is closely aligned with the quasi-steady analytical model and the experiment method. During this period, the frequency of flow separation is relatively low, and there are no complex wake structures in the domain. However, the CFD model enters a steady tumbling motion after approximately 1 s (for a detailed explanation of the steady tumbling motion, refer to Section 3). This occurrence is relatively earlier compared to what is observed in the quasi-steady analytical model and experimental observations. It could be attributed to factors such as the approximation assumptions, numerical methods, and the boundary condition selection employed in the CFD model.

3. Effect of MOI

3.1. Quasi-Steady Analytical Model

Let the moment of inertia of the wing in the Y direction with units of kg·m${}^2$ be set to 0.00045, 0.0006, 0.0008, 0.001, 0.002, 0.004, 0.006, 0.008, and 0.009, respectively, while the mass and shape of the wing remain unchanged. By substituting these states into the quasi-steady analytical model and simulating at the same initial state, the trajectory of the freely falling wing can be observed as shown in Figure 10. In order to more clearly illustrate the simulation of the wing tumbling process, the change in orientation of the wing in the free-fall path is shown in Figure 11. This is similar to the observations made by Andersen et al. in a small, freely falling plate, as reported in [24].

As can be gleaned from Figure 10 and Figure 11, it is evident that, as the MOI about the $Y^{\prime }$-axis of the wing increases, the falling trajectory of the wing becomes steeper, resulting in a greater loss of height. Concurrently, the angular rate of tumbling $\dot{\theta }$ of the wing also decreases. In the phase diagrams of Figure 11b,e, the tumbling of the wing exhibits stable limit cycle oscillation. After entering the falling motion, the phase trajectory from the initial state is attracted to the limit cycle, eventually undergoing periodic motion within the limit cycle. As the moment of inertia (MOI) of the wing increases, the process of convergence from the initial state to the limit cycle becomes more protracted. In phase space, the phase locus takes longer to converge to the limit cycle in the $\dot{\theta }$ direction relative to the $v_{x^{\prime }}$, $v_{z^{\prime }}$ direction, and the average angular velocity decreases, which is a consequence of the increased MOI. From Equation (5), it can be deduced that the circulation generated by the angular motion decreases, resulting in a decrease in $F_z$, thereby augmenting the sink rate of the wing.

3.2. Experiment Results

The falling experiment was conducted using the experimental setup depicted in Figure 5 and Figure 6. The position of the clump weights were adjusted to alter the wing’s moment of inertia (MOI), and the wing was released at a height of 5 m with an initial attitude angle of ${\theta }_0=-{75}^{\circ }$. The acceleration, angular velocity, and attitude angle were measured using an inertial measurement unit (IMU) and they were used to derive the trajectory of the wing’s descent, as shown in Figure 12. The external forces acting on the wing were also estimated.

It can be seen that, as the wing’s MOI increases, the falling trajectory becomes steeper and the transition phase from the beginning of the fall to the stable tumbling phase becomes longer, which is consistent with the conclusions obtained from the analytical model in the previous sections. The experimental results also indicate that the analytical model can accurately predict the tumbling behavior of a freely falling wing. This is likely due to the fact that some of the coefficients in the analytical model are determined by experiments conducted in previous work, allowing the model to capture some of the three-dimensional effects during the falling motion. However, discrepancies between the simulated and experimental data may be attributed to sensor noise, the analytical model’s inability to account for complex nonlinear and unsteady aerodynamic effects, and external perturbations in the falling experiment. The wing tumbling motion in the experiment, shown in Figure 12, did not occur strictly in the two-dimensional plane, also due to these perturbations.

4. Effect of COM Position on Freely Falling Wing

The coordinate system $O{X}^{\prime }clump$ is established at the center of the stalk in Figure 3, with the $O{X}^{\prime }clump$ axis aligned with the $O^{\prime }{X}^{\prime }$ axis. The coordinates of the two clump weights on the $O{X}^{\prime }clump$ axis are denoted as $C_{xl}$ and $C_{xr}$, respectively.

By altering the relative position of the counterweight on the stalk in Figure 3, the position of the COM can be adjusted while maintaining the mass constant. The relationship between the clump weights’ position, the COM position $X_c$, and the moment of inertia $I_c$ is presented in

Table 1. The relationship between the coordinates of the clump weights, the COM position, and the MOI.

Case	${\mathit{C}}_{\mathit{xl}}$ (mm)	${\mathit{C}}_{\mathit{xr}}$ (mm)	${\mathit{X}}_{\mathit{c}}$ (mm)	${\mathit{I}}_{\mathit{c}}$ (kg·m${}^2$)
1	−400	500	4.547	0.008
2	−300	500	10.126	0.007
3	−125	500	19.89	0.006
4	100	500	32.44	0.006
5	250	500	40.813	0.006
6	400	500	49.138	0.008

:

4.1. Simulation of the Analytical Model

By incorporating the above model parameters into the analytical model, the simulation trajectory can be obtained, as depicted in Figure 13. The phase diagrams of the freely falling wing’s $v_{x^{\prime }}$, $v_{z^{\prime }}$, and $\dot{\theta }$ can also be generated.

According to the simulation results, it can be observed that, as the COM position moves forward, the transition phase from the initial state to the limit cycle of tumbling motion takes longer and the phase trajectory of the transition phase becomes more complex. As seen in Figure 13, when the COM position is shifted approximately 40 mm from the geometric center of the wing, corresponding to case 5, a new limit cycle appears in the phase space. The projection of the phase trajectory in the $v_{x^{\prime }}$, $v_{z^{\prime }}$ plane and the $v_{z^{\prime }}$, $\dot{\theta }$ plane in Figure 13f is shown in Figure 14, respectively. It can be observed that the curve starts from the initial state and quickly enters the new limit cycle, undergoing quasi-periodic motion in the ring. The stability of the $v_{x^{\prime }}$, $v_{z^{\prime }}$ direction within the limit cycle is evident, while the $\dot{\theta }$ direction is not stable. As a result, the phase trajectory eventually leaves the new limit cycle as the number of oscillations increases, as shown in Figure 14b. After leaving the new limit cycle, the curve is attracted to the limit cycle corresponding to the original tumbling motion and eventually converges to this limit cycle. The falling motion of the wing is characterized by initial fluttering, with the flutter angular velocity’s amplitude increasing before the wing enters the tumbling falling phase. The trajectory of the wing’s fall is depicted in Figure 15. Figure 15a shows the fluttering of the wing during the fall represented by the newly emerged limit cycle, and Figure 15b shows the transition from fluttering to tumbling motion. As the COM moves further forward, the limit cycle corresponding to the tumbling motion disappears and the wing falls while fluttering periodically according to the limit cycle shown in Figure 13g.

4.2. Experiment Results

In this part, the experimental method was employed to observe the actual falling trajectory of the wing after release. Using the wing shown in Figure 3, the position of the clump weights block is varied according to the scheme in Table 1, and the falling experiments are conducted for Case 1 to Case 6, respectively. The measured falling trajectory of the wing is shown in Figure 16.

In the experiment, due to the asymmetrical disturbance on the wing, the actual falling trajectory is difficult to maintain in the $XZ$ plane, and the wing typically falls along a helix. The tumbling motion is depicted in Figure 16 as the wing tumbles around the $Y^{\prime }$ axis of the body coordinate system in the helix. It can be seen that, for case 1 to case 5 (Figure 16a–e), the tumbling motion occurs before the wing touches the ground, with the transition phase from the initial state to the tumbling motion of the wing becoming longer as the position of the COM moves forward. However, for case 6 (Figure 16f), the wing does not tumble until it touches the ground and falls more smoothly along the helix. Projecting the spiral line falling trajectory in the transition phase of the wing into the two-dimensional plane is similar to the trajectory corresponding to case 5 in Figure 15, which demonstrates a heaving/pitching composite motion along the spiral line in the actual three-dimensional motion. The discrepancy between the experimental trajectory and the trajectory simulated by the two-dimensional analytical model is attributed to the asymmetry of the wing, the asymmetry of external disturbance, and the asymmetrical shedding of vortex on the wing surface during the actual falling process.

5. Results and Discussion

5.1. Quasi-Steady Analytical Model Results

1.

Effect of MOI: In the case of wings with different MOIs, there exists a limit cycle of tumbling motion that all wings will eventually converge to after releasing, reaching a stable state shown in Figure 11. The phase from the initial state to the final stable tumbling limit cycle is defined as the transition phase. As the MOI of the wing increases, the trajectory of the free-falling wing becomes steeper, and the transition phase becomes longer, meaning that the wing will take more time to converge to the limit cycle.

2.

Effect of COM: When the COM moves forward, the transition phase from the initial state to the tumbling motion limit cycle becomes longer, and the phase trajectory during the transition phase becomes more complex. A new limit cycle emerges when the COM is positioned 40 mm ahead of the wing’s geometric center, as can be seen in Figure 13f, corresponding to the quasi-periodic fluttering motion of the wing. At this point, the limit cycle of the fluttering motion is unstable, and the phase trajectory eventually diverges from this limit cycle, converging towards the limit cycle of the tumbling motion; Figure 15 shows this process well. As the COM continues to move forward, the limit cycle corresponding to the tumbling motion disappears, and the limit cycle of the fluttering motion becomes a stable limit cycle, resulting in the wing freely falling with periodic fluttering motion.

5.2. Experimental Results

1.

Effect of MOI: By symmetrically altering the position of the clump weights block, the MOI of the wing is changed, and tumbling motion is observed under different conditions of MOI, as shown in Figure 12. As the MOI of the wing increases, the trajectory of the freely falling wing becomes steeper, and the transition phase from the initial descent to stable tumbling becomes longer. This observation aligns with the conclusions drawn from the analytical model.

2.

Effect of COM: Experimental results show that maintaining the actual wing trajectory in the $XZ$ plane is challenging because of asymmetrical disturbances. Typically, the wing descends along a spiral path, as shown in Figure 16. As the COM moves forward, the transition phase from the initial state to tumbling motion increases. Further forward movement of the COM results in the fluttering motion of the wing. However, no tumbling motion is observed in the experiment after shifting the COM forward by 40 mm, as shown in Figure 16d.

5.3. Comparison of the Analytical Model with Experimental Results

The quasi-steady analytical model accurately reveals the influence of changes in MOI and COM on the freely falling motion of the wing. These observed motion patterns are consistent with the conclusions obtained from the quasi-steady analytical model. However, due to unavoidable asymmetrical disturbances on the wing during experiments, maintaining the wing’s trajectory within a two-dimensional plane becomes challenging. Additionally, obtaining an accurate trajectory of the wing is difficult due to sensor errors.

6. Conclusions

In this study, the dynamics of a freely falling wing in the context of air-drop launch technology for UAVs was investigated. The mass, MOI, and COM position of the wing have a significant impact on the flight characteristics of UAVs; thus, this study examined the effect of changes in the moment of inertia and center of mass position on the falling motion while holding the mass constant. The aerodynamic forces acting on a freely falling wing are complex, nonlinear, and unsteady, and exhibit strong coupling with the wing motion, making them challenging to accurately model numerically. Extensive research has been conducted on the dynamics of falling objects, resulting in the development of various numerical simulation methods and simplified descriptive models. However, these studies largely focus on the free-fall of small objects with simple geometric shapes such as flat plates and discs in a laboratory setting, rather than larger objects such as wings. In this study, the Andersen–Pesavento–Wang analytical model was adapted to the analysis of the freely falling wing, and a correction term was added to account for the center of mass position. The quasi-steady analytical model was validated by CFD and experimental results. This study investigated the effects of the MOI and COM position on the dynamics of a freely falling wing, and drew the following conclusions:

1.

A quasi-steady analytical model was developed based on the Andersen–Pesavento–Wang model. The analytical model was derived in the two-dimensional plane and can reflect the dynamic behavior well in the longitudinal plane during the three-dimensional falling of a real wing.

2.

After deployment, the wing undergoes a transition phase before eventually entering a stable motion, which can be characterized as either tumbling or fluttering, as revealed by the quasi-steady analysis model and experiment method employed in this study. Through CFD analysis, it is observed that, during the transition phase, the shedding frequency of vortices gradually increases and stabilizes, resulting in periodic aerodynamic oscillations due to the cyclic generation and separation of vortices. These periodic oscillations eventually lead to the emergence of a stable periodic motion. Both of these stable motions exhibit relatively low translational velocities. In the case of rolling motion, the pitch rate of the wing remains approximately 15 rad/s. In the case of fluttering motion, the pitch rate oscillates between positive and negative values. Therefore, both of these stable motion patterns should be avoided in the air-drop launched UAVs.

3.

For a wing with elliptical airfoil, the higher the MOI about the $Y^{\prime }$-axis, the steeper the trajectory, the lower the angular rate of tumbling, and the shorter time for it staying in the air. For most aircraft, the pitch tumbling is not beneficial. Without changing the shape of the wing and the position of the COM, and only changing the MOI of the wing about the $Y^{\prime }$-axis, the initial drop attitude cannot prevent the pitch tumbling of the wing during the falling process.

4.

The position of the COM has a crucial influence on the handling performance and stability of the UAVs. The traditional flight dynamics theory suggests that the forward shift of the COM position can make the vehicle obtain better static stability. For the freely falling wing, the forward shift of the COM will delay the appearance of tumbling, i.e., the transition phase becomes longer. When the forward shift of the COM exceeds a certain value, a new relatively stable motion of falling appears, which is expressed as fluttering in the analytical model, while the wing shows a heaving/ pitching composite motion along the spiral line in the real three-dimensional falling experiment. In the future, for air-drop launch UAVs with relaxed longitudinal static stability, the possibility of its tumbling needs to be considered.

In the future, the analytical model can be considered to include corrections for different airfoil types and wing shapes, and to include the effects of rudder surfaces, so that the model can describe the freely falling process of a more general aircraft. Since the model also uses two-dimensional quasi-steady assumptions, extending the model to three-dimensional space and considering the unsteady aerodynamic force for better prediction of falling motion is worth studying. In addition, the analysis of the uncertainties of this complex process is a research direction of interest.