


This chapter presents a comprehensive study on optimising the DI design for the AeroTail to maximise its aerodynamic torque for robotic balancing; it demonstrates a specialised Computational Fluid Dynamics (CFD) simulations on a selection of shapes to explore how specific geometrical features influence aerodynamic drag, with detailed explanations of the observed effects; and it introduces an optimisation framework to identify the optimal frame design that maximises aerodynamic drag, while adhering to predefined constraints such as weight and dimensional limits. Finally, the chapter includes the implementation of practical experiments using the AeroTail test platform, equipped with the optimised DI designs, to assess their real-world performance on a rotational tail mechanism with variable wind speeds along its length. This involves a comparative analysis of the experimental findings against the CFD simulation results and the establishment of an analytical approach for quantifying the performance variation among different DI models.

\section{Methodology for drag-inducer design}\label{sec_method}

The primary objective in designing the DI is to achieve a balance between minimising weight and maximising the generated drag force. The first half of the objective can be simplified to finding the minimum material use for the application and then optimising the distribution; such that a fixed available weight for the material is maintained during the design process. In order to find the optimal DI design, two distinctive steps need to be taken; firstly, the theoretical best shape to generate aerodynamic drag is to be found through CFD simulation since there is insufficient related conclusions in previous studies; secondly, using this general shape, the optimal DI design for maximising torque generated aerodynamically can be calculated using a optimisation problem. 

For the first step, there are initial assumptions needed to be made for the simulation setup. From \eqref{aero_drag}, it can be concluded that there are two strategies of maximising drag; one is to maximise the frontal area ($A$), and the other is doing that for the drag coefficient ($C_{d}$). However, these two methods cannot be examined separately because both $A$ and $C_{d}$ affect drag proportionally. An important finding from previous studies is that cup-like concave shapes yield better drag coefficients than their convex or flat counterparts \cite{hoerner1965practical,aziz2008online}. To make the geometry of the shapes in question easily describable, we define a concave shape to be that with an 2 dimensional ``opening" surrounded by a rim and the rest of the body extending into the third dimension. If such shape is to be used to generate drag thus bearing load, it is usually structurally sound to only have rigid material in the rim to form a ``frame" while the rest of the shape can be made from less strong or even soft film-like materials, forming the ``film". In this setup, the effective area ($A$) is determined by the size of the frame's opening, while the specific shape of the film can significantly impact the drag coefficient ($C_{d}$). Additionally, the mass is primarily influenced by the frame material, as the film can be made from lightweight materials, thereby reducing the overall device mass. This `frame-film' setup enables us to explore drag-maximising strategies without introducing too many uncontrolled variables.

Even though a formulated optimisation problem is not necessary, there are some constraints to act as the controlled variables in the simulation studies. As mentioned before, the available mass for the frame must be a given and fixed value according to the condition of the overall device. For the simplicity of the study, the frame is envisaged as an extruded, cage-like structure with a uniform cross-section. The available length for the frame $l$ can be calculated using the following equation assuming the material has uniform and constant density:
\begin{equation}\label{frame_length}
    l = \frac{m_{0}N}{a_{0}\rho},
\end{equation}
where $m_{0}$ is the total mass of the balancing mechanism, $N$ is a design factor which determines the percentage of the mass that can be used on frame, $a_{0}$ is the area of the cross-section and $\rho$ is the density of the material used. 

When conducting a study about maximising drag with $A$ and $C_{d}$ being the subjects, it is only logical to set $\rho$ and $v$ as constants according to \eqref{aero_drag}. In this case, because $\rho$ is simply the density of air, is it assumed to be constant; flow speed $v$, however, must be a value calculated from the torque required ($\tau_{r}$) for the device to complete a specific task. Assuming the drag acts on a single point, the aerodynamic torque can be simply calculated as $\tau_{r}=Dl_{0}$ where $l_{0}$ is the lever arm. Therefore, the flow speed can be expressed using this equation:
\begin{equation}\label{flow_speed}
    v = \sqrt{\frac{2|\tau_{r}|}{\rho C_{d}Al_{0}}} \ for \ v>0,
\end{equation}
in which $v$ is the flow speed needed to aerodynamically produce the required torque $\tau_{r}$.

In \secref{sec_sim_cfd}, the CFD simulation study will be conducted according to the aforementioned methodology for undergoing the first step. In the simulation study, specific values from AeroTail-based scenarios are used in \eqref{frame_length} and \eqref{flow_speed} to obtain the constant setup for the simulation environment. As for the second step, an optimisation problem shall be setup according to the observations from the simulation. The optimisation problem aims to directly or indirectly maximise the aerodynamic drag that could be generated by the DI in a given scenario. The results from the optimisation problem can be one or several DI shape with mathematically definable parameters and they will undergo further examination until the best option is chosen. The details of this step are presented in \secref{sec_opt_frame}.



\section{CFD Simulation Study} \label{sec_sim_cfd}

\subsection{CFD Simulation Setup} \label{cfd_set}

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{pics/All shapes_1.jpg}
    \caption{(a)-(r) depict all the shapes involved in the initial phase of the CFD study. (s) illustrates the frontal view of these shapes, with the effective area marked in green and the frame outlined in black. (t) is a flat rectangular shape that has accurate wind tunnel data which can be used for result validation. A specific comparison for the Semi-cylinders shown in (l) is conducted, focusing on the height-to-length ratio ($H/L$), which varies while maintaining the overall configuration, thereby altering the shape's slenderness or fullness.}
    \label{all shapes}
\end{figure}

Even though the main purpose of the simulation study is to find the general characteristics of a shape that enables it to produce larger aerodynamic drag, the initial conditions and assumptions do need to be based on the AeroTail design \cite{sun2023aerotail} in the previous work since it is easier to evaluate the outcome given a scenario with fitting application. The estimated mass of the AeroTail is approximately $5$ g, and it is targeted that each DI's opening frame should not exceed $10\%$ of this total mass. The frame is assumed to be made from Tough 2000 Resin, a widely used material for 3D printing (which is later on used for the DI prototype in the \secref{sec_ex}) with a density of $1.11$ g/cm$^3$. The frame is set to have a cross-sectional diameter of $2$ mm (and its area is $3.14$ mm$^2$) which is the smallest value possible for the structure to stay integral during manufacturing and under load judging from past experience with the material. Substituting $m_{0}=5g$, $N=0.1$, $a_{0}=3.14$ mm$^2$, and $\rho=1.11$ g/cm$^3$ into \eqref{frame_length}, the maximum permissible frame length is calculated to be $14.3$ cm. Consequently, a total usable frame length of $12$ cm is selected for its mathematical convenience.

The initial stage of the simulation study focuses on shapes subjected to an opening circumference constraint of $12$ cm, which is the available frame length, since it is postulated that the frame structure is continuous and confined to the rim of the opening. Regular shapes are categorised based on the number of edges, and the most logical and intuitive film shapes are conceptualised for each (as depicted in \figref{all shapes}). These shapes are modelled as solid shells with a consistent thickness of $1$ mm. The setup specifications can be found in \tabref{tab_sim_setup}. The mesh's automatic grain size is set to $10$ mm (with a maximum of $20$ mm). The mesh around the shape's walls was refined to capture complex geometric features accurately. The grain sizes are obtained through a convergence test as part of the mesh independence study. The boundary walls were placed $150$ mm away from the shape, with the inlet positioned $100$ mm in front and the outlet $200$ mm behind the shape forming a typical rectangular cuboidal mesh. These distances were established through a sensitivity analysis. The boundary walls have a no-slip condition while the inlet has a set flow velocity. In order to choose an appropriate flow speed for the simulation, we assume the drag-induced torque of $0.0111$ Nm, which is the torque sufficient to reorient a $23$ g device from a $30^{\circ}$ angle to a vertical position using a DI with a frontal area of $5$ cm$^2$. A Semi-sphere shaped DI is selected as the benchmark, as its drag coefficient data is supported by plenty of existing wind tunnel tests. Substituting $\tau_{r}=0.0111$ Nm, $A=5$ cm$^2$, $C_{d}=1.43$ (for semi-sphere) and $l_{0}=10$ cm (the lever arm for AeroTail) into \eqref{flow_speed}, the flow speed required for the task is calculated to be $16$ m/s. We utilised the  k-$\epsilon$ turbulence model for its robustness and suitability and the simulation stops after a set number of iterations to ensure steady-state solution.

\begin{table}
\centering
\caption{Simulation set-ups}
%\vspace{3mm}
\begin{tabular}{lccl}
\toprule
\vspace{1mm}
Setups & Specification\\
\hline 
Element type \rule{0pt}{2.6ex} & Unstructured \\
%
Element size  & $10$ mm \\
%
Maximum element size  & $20$ mm \\
%
Boundary walls & $150$ mm $\times$ $150$ mm \\ 
%
Inlet/Outlet & $100$ mm / $200$ mm \\
%
Boundary wall condition & No-slip \\
%
Time step & Auto-generated \\
%
Turbulence model & k-$\epsilon$ \\
%
Convergence criteria & Specified iteration \\
\bottomrule
\end{tabular}
\label{tab_sim_setup}
% \vspace{-2mm}
\end{table}


\subsection{CFD Study Results} \label{cfd_result}

\begin{table}
\centering
\resizebox{\columnwidth}{!}{%
\begin{threeparttable}
\caption{Overview and Results of Initial CFD Simulation for Various Shapes}
%\vspace{3mm}
\label{tab1}
\begin{tabular}{lcccccl}
\toprule
\vspace{1mm}
Shapes & Area (cm$^2$) & Force (N) & $C_{d}$ & $C_{d}*$ & Sources\\
\hline 
(a) Semi-sphere \rule{0pt}{2.6ex}& $11.46$ & $0.257$ & $1.43$ & 
\begin{tabular}{@{}c@{}} $1.42$ \\ $1.40$ \end{tabular} &
\begin{tabular}{@{}c@{}} \cite{aziz2008online} \\ \cite{white1966fluid}\end{tabular} \\ 
%\hline 
(b) Cone & $11.46$ & $0.244$ & $1.36$ & $2.10$ & \cite{aziz2008online},\cite{munson1995fundamentals}\\
%\hline
(c) Hollow cylinder & $11.46$ & $0.237$ & $1.32$ & -\\
%
(d) Tetrahedron & $6.93$ & $0.147$ & $1.35$ & -\\
%
(e) Triangle-prism & $6.93$ & $0.120$ & $1.10$ & -\\
%
(f) Pyramid & $9.00$ & $0.190$ & $1.34$ & (i) \\
%
(g) Angled hollow cube & $9.00$ & $0.191$ & $1.36$ & (i) \\
%
(h) The ``hat" & $9.00$ & $0.198$ & $1.40$ & -\\
%
(i) The ``book" & $9.00$ & $0.190$ & $1.34$ & $2.20$ & \cite{hoerner1965practical} \\
%
(j) Hollow cube & $9.00$ & $0.183$ & $1.30$ & $2.05$ &
\cite{aziz2008online},\cite{munson1995fundamentals}\\
% 
(k) Semi-cylinder shallow & $8.00$ & $0.175$ & $1.40$ & $2.30$ & 
\cite{hoerner1965practical},\cite{white1966fluid}\\
%
(l) Semi-cylinder & $9.00$ & $0.195$ & $1.39$ & (k) \\
%
(m) Semi-cylinder deep & $8.00$ & $0.171$ & $1.36$ & (k)\\
%
(n) Elongated cylinder & $9.00$ & $0.193$ & $1.34$ & (k)\\
%
(o) Open Semi-cylinder & $9.00$ & $0.188$ & $1.33$ & (k)\\
%
(p) Long triangle & $4.90$ & $0.095$ & $1.24$ & -\\
%
(q) Capsule & $10.2$ & $0.230$ & $1.44$ & (a) \\
%
(r) Rhombus & $7.79$ & $0.164$ & $1.34$ & -\\
%
(t) Plank & $4.00$ & $0.076$ & $1.23$ & 
\begin{tabular}{@{}c@{}} $1.17$ \\ $1.20$ \end{tabular} & 
\begin{tabular}{@{}c@{}} \cite{hoerner1965practical} \\ \cite{mohammed2016theory} \end{tabular}\\
\bottomrule
\end{tabular}
\begin{tablenotes}
    \item $^*$ This column presents experimental data from previous sources. If a single value with multiple sources is listed for a shape, the sources concur on this value. When the $C_{d}$ refers to a different shape, it implies that the shapes are sufficiently similar, differing only in their dimensions.
\end{tablenotes}
\end{threeparttable}
}

% \vspace{-2mm}
\end{table}

The outcomes of the CFD simulation study are collated in \tabref{tab1}, with the shapes named according to their corresponding letter markings in \figref{all shapes}. Given the opening circumference constraint, there is a significant variance in the effective areas of these shapes. Both the drag force and drag coefficient values are computed by the simulation software. Analysis of these results reveals several insights. Predominantly, shapes with a larger effective area tend to generate greater drag, as the drag force is directly proportional to this area. However, the drag coefficients of these shapes exhibit relative unpredictability, fluctuating between $1.1$ and $1.5$, while effective areas may differ markedly. Therefore, optimising the effective area is crucial in identifying the shape capable of producing the highest drag. When comparing shapes (a) and (q), (e) and (p), or (k) and (l), it is evident that shapes with a slimmer frontal profile generally yield higher $C_{d}$ values. Moreover, shapes with closed ends typically exhibit larger $C_{d}$ than their open-ended counterparts, as seen in comparisons between (l) and (o), or (g) and (i). Additionally, an analysis of (a) versus (c) and (h) versus (j) indicates that arch-shaped profiles usually generate the highest $C_{d}$ for a given frontal shape.

The final column of \tabref{tab1} presents wind tunnel test data from previous studies. The simulation conditions for the semi-sphere (a) and plank (t) closely match those of the wind tunnel tests, leading to comparable $C_{d}$ outcomes, thereby affirming the simulated results' reliability. However, the remaining shapes display considerably different $C_{d}$ values. This discrepancy arises because the simulation setup, which more accurately reflects practical scenarios with shapes not in contact with boundary walls, differs from the wind tunnel setup (as outlined in \secref{sec_intro}). This difference confirms that the shapes under consideration vary significantly from the 2D ``prism" shapes in previous studies, whose $C_{d}$ values are not directly applicable for practical calculations.

\begin{table}
\centering
\caption{CFD results for the Semi-cylinders}
%\vspace{3mm}
\begin{tabular}{lcccl}
\toprule
\vspace{1mm}
$H/L$ (cm) & Area (cm$^2$) & Force (N) & $C_{d}$ \\
\hline 
$2.0/4.0$ \rule{0pt}{2.6ex} & $8.00$ & $0.175$ & $1.400$  \\ 
%\hline 
$2.4/3.6$ & $8.36$ & $0.188$ & $1.390$ \\
%
$2.7/3.3$ & $8.91$ & $0.193$ & $1.384$ \\
%
$2.9/3.1$ & $8.99$ & $0.195$ & $1.382$ \\
%
$3.0/3.0$ & $9.00$ & $0.194$ & $1.370$ \\
\bottomrule
\end{tabular}
\label{tab_semis}
% \vspace{-2mm}
\end{table}

Examining the Semi-cylinders with varying width-to-length ratios allows for a more distinct understanding of the impact of both area and drag coefficient on drag force. \figref{all shapes} (l) and \tabref{tab_semis} provide a comparative analysis of five Semi-cylinders, each sharing the same opening circumference but differing in width-to-length ratios, and consequently, in effective areas. The data indicates that as the shapes become more slender, their frontal area diminishes, yet there is a consistent trend of increasing drag coefficients. A $H/L$ ratio of $2.9/3.1$ yields the highest drag force, as the benefits of a higher $C_{d}$ outweigh the drawbacks of a smaller area. However, as the shape becomes increasingly slender, the diminishing frontal area begins to adversely affect the drag force, leading to a reduction in overall drag.

To summarise the previous paragraphs, the findings can be organised into these main points:
\begin{itemize}
    \item Shapes with a larger effective area tend to generate greater drag while the drag coefficients of these shapes exhibit relative unpredictability.
    \item  Shapes with slimmer frontal profiles, closed ends, and arch-shaped side profiles generally yield higher $C_{d}$ values. 
    \item The $C_{d}$ values of 2D ``prism" shapes in previous studies are not directly applicable for practical calculations.
\end{itemize}
The variation in drag coefficients among different shapes can be attributed to the phenomena of flow separation and vortex formation. For instance, comparing shapes (a) and (e), their streamline images from the simulation, depicted in \figref{stream}, reveal notable differences. The extended blue (low velocity) area behind shape (a) compared to (e) suggests longer flow separation, which in turn leads to the formation of vortices, observable as circular vectors around and behind the shapes. Larger and more intense vortices are generally indicative of higher drag \cite{ponitz2014aerodynamics, khosravi2017aerostructural}, corroborating the higher drag coefficient observed for shape (a). Additionally, the shape's wind-capturing capability plays a vital role, as these vortices originate from air pressure disparities at the shape's edges \cite{amini2019suppressing}. This elucidates why arch-shaped profiles and shapes with closed ends exhibit higher $C_{d}$: they effectively create higher air pressure in the pocket regions, resulting in a larger pressure differential with the surrounding airflow.

\begin{figure}
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/streamlines.jpg}
    \caption{The CFD simulation's streamline (left) and vector (right) images. The colour represents flow velocity with red denoting high and blue low velocities, indicated by the colour bar on the right. At the top is the Triangle-prism (e), and at the bottom, the Semi-sphere (a).}
    \label{stream}
\end{figure}


\section{Optimal Frame Design} \label{sec_opt_frame}

\subsection{Methodology of Frame Design}

In previous designs discussed in \secref{sec_intro}, the DI frame was limited to `closed' structures. However, if the primary constraint is the total frame length and the goal is to maximise effective area, a closed design can be inefficient. Our approach considers these limitations to develop a more material-efficient frame design. An alternative approach involves using segmented frames on the rim, each segment connected to the frame's centre by a ``rib", as illustrated in \figref{frames}. This design is feasible as the film material does not require support along its entire edge; supporting it in sections is sufficient for structural stability. In this configuration, $n$ represents the number of ribs in a polygonal frame, $a$ is the length of a rib, and $b$ is the length of a rim segment to which the film is attached. As the ``T" shapes are evenly spaced, the effective area typically forms regular polygons.

\begin{figure}
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/frame definition.jpg}
    \caption{The structures of segmented frames with varying numbers of ribs. Each ``T" shape, comprising a rim segment and a rib, extends outward from the centre (indicated by a red dot) to form a carousel configuration. The unsupported film edges are represented by dotted lines, and the effective area is highlighted in green.}
    \label{frames}
\end{figure}

The area of such a polygon is given by the equation:
\begin{equation} \label{area}
    \begin{split}
        A = & n(A_{1}+2A_{2}-2A_{3}) \\
        = & n(\cos(\frac{\pi}{n})a\cdot\sin(\frac{\pi}{n})a+2\sin(\frac{\pi}{n})a\cdot\sin(\frac{\pi}{n})(\frac{b}{2}) \\
        & -\cos(\frac{\pi}{n})(\frac{b}{2})\cdot\sin(\frac{\pi}{n})(\frac{b}{2})) \\
        = &  n(\frac{1}{2}\sin(\frac{2\pi}{n})(a^2-\frac{b^2}{4})+2\sin^2(\frac{\pi}{n})(a\cdot\frac{b}{2})),
    \end{split}
\end{equation}
where $A$ is the effective area. The variables $n$, $a$ and $b$ are interconnected as follows:
\begin{equation} \label{abln}
    a=\frac{l}{n}-b,
\end{equation}
where $l$ is the total available length of the frame ($12$ cm in this study). Consequently, Equation \eqref{area} can be reformulated as a function of $n$ and $b$, with $l$ as a constant:
\begin{equation} \label{area2}
    \begin{split}
         A(n,b)= & n(\frac{1}{2}\sin(\frac{2\pi}{n})(\frac{3}{4}b^2-\frac{2l}{n}b+\frac{l^2}{n^2}) \\
        & +\sin^2(\frac{\pi}{n})(\frac{l}{n}b-b^2)).
    \end{split}
\end{equation}

In the top-left corner of \figref{frames}, three additional parameters, $w$, $h$, and $d$, are defined. These will serve as constraints in the optimisation problem. A coordinate system is placed at the frame's centre (red dot), with the x-axis aligned with one of the ``T" shapes (in yellow) and the y-axis perpendicular to it. $h$ is defined as the length of the frame in the x-direction from its centre to the opposite end of the yellow ``T". $w$ is the frame's width in the y-direction, and $d$ is the distance between the ends of two adjacent rim segments. These parameters can be calculated using the following equations:
\begin{equation}
\begin{split}
     w = & 2\sin(\pi-\frac{2\pi}{n})a+\cos(\pi-\frac{2\pi}{n})b \\
     h = & \cos(\pi-\frac{2\pi}{n})a+\frac{1}{2}\sin(\pi-\frac{2\pi}{n})b \\
     d = & 2\sin(\frac{\pi}{n})a-\cos(\frac{\pi}{n})b.
\end{split}
\end{equation}


\subsection{Optimisation}

For simplistic structures, the design optimisation can be approached by a step-based qualitative method \cite{lai2011research}; however, since the properties of the frame have being quantified, a classic constrained optimisation problem (inspired by \cite{guodong2018kinematics}) is setup for the frame design. The aim of the optimisation problem is to identify the optimal values for $n$ and $b$ that maximise $A$, the effective area:
\begin{equation} \label{opt_eq}
    \begin{aligned}
        \max_{n,b} \quad & A(n,b) \\
        \textrm{s.t.} \quad & 2 \leq n , n\subset\mathbb{Z} \\
        & 0.2 \leq b  \\
        &  w \leq 5, d > 0 \\
        & \frac{l}{n}-b+h \leq 5
    \end{aligned}
\end{equation}
where $A(n,b)$ is defined in Equation \eqref{area2}. Both $n$ and $b$ are subject to specific constraints, and additional constraints are imposed by the parameters $w$, $h$, and $d$. The minimum value for $n$ is set at $2$, as having fewer than two ribs does not provide a structurally stable support for the film. Furthermore, $n$ must always be an integer. The lower bound for $b$ is established at $0.2$ cm, to ensure that the length of a single rim segment exceeds the frame's cross-sectional thickness. The frame's width $w$ is constrained by the prototype platform, limited to a maximum of $5$ cm. Similarly, the sum of $a$ and $h$ cannot exceed $5$ cm due to the physical limitations of the prototype (as indicated in the last line of \eqref{opt_eq}). Lastly, to avoid the rim segments intersecting or overlapping, the distance $d$ between their tips must always be greater than zero.

The optimisation problem is solved using the ``fmincon" function in MATLAB. The results indicate that the optimal value for $n$ is $2$, with a corresponding $b$ of $3.50$ cm, yielding a maximum area of $17.5$ cm$^2$. Although this result is definitive, outcomes for $n$ values of $3$, $4$, and $5$ are also documented in \tabref{tab_opt} for additional analysis. Notably, when the frame has five ribs, $b$ reaches its lower bound, rendering any further increase in $n$ impractical. It is important to recognise that while the effective area is a primary factor influencing aerodynamic drag, the drag coefficient also plays a significant role. In the case of the two-rib ``I" shape frame, different configurations of the resultant film shape can lead to substantially different drag coefficients. Therefore, further CFD simulations are essential to validate the superiority of the optimised frame designs over those proposed in \secref{cfd_set} and to facilitate a comparative analysis of the CFD results from \secref{cfd_result} with those of the optimised shapes.

\begin{table}
\centering
\caption{The parameter optimisation results. The drag force and drag coefficient of each shape are obtained from further CFD simulations.}
\vspace{3mm}
\begin{tabular}{lcccccl}
\toprule
\vspace{1mm}
Shapes & $n$ & $b$ (cm) & $A$ (cm$^2$) & Drag (N) & $C_{d}$\\
\hline 
(a) \rule{0pt}{2.6ex} & $2$ & $3.50$ & $17.50$ & $0.403$ & $1.47$  \\ 
%\hline 
(b) & $2$ & $3.50$ & $17.50$ & $0.394$ & $1.44$ \\ 
%\hline 
(c) & $3$ & $1.55$ & $15.57$ & $0.323$ & $1.32$ \\
%\hline
(d) & $4$ & $0.50$ & $14.88$ & $0.327$ & $1.40$ \\
%
(e) & $5$ & $0.20$ & $12.24$ & $0.267$ & $1.39$ \\ 
\bottomrule
\end{tabular}
\label{tab_opt}
\vspace{-2mm}
\end{table}

\subsection{Further CFD Simulation Study} \label{sec_further_sim}

\begin{figure}
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/polygon_cylinder.jpg}
    \caption{The shapes utilised in the CFD studies, inclusive of both the frames and the film behind them. The sizes are not to scale but are presented to illustrate the structural configurations. (a) Configuration 1 for $n=2$, with arches on the shorter edges; (b) Configuration 2 for $n=2$, with arches on the longer edges; (c) for $n=3$; (d) for $n=4$; and (e) for $n=5$.}
    \label{swshapes}
\end{figure}

The CFD study encompasses five configurations, as shown in \figref{swshapes}, each corresponding to the optimisation results listed in \tabref{tab_opt}, with two scenarios for $n=2$. The film shapes for these frames are designed in line with the principle of maximising drag given the specific frame structure. Therefore, they are fashioned with arch-shaped profiles and closed ends, aligning with the findings from \secref{cfd_result}. The conditions applied to these CFD simulations are consistent with the first stage outlined in \secref{cfd_set}.

The results indicate that configuration (a) produces a drag force of $0.403$ N, marginally surpassing the $0.394$ N of configuration (b), despite both having identical effective areas. This difference stems from the higher drag coefficient observed when arches are positioned on the shorter edges. As $n$ increases, there is a noticeable trend of diminishing drag force, with variations in drag coefficients. An exception is observed with shape (d), which, despite a smaller frontal area than (c), generates a higher drag force due to its substantially greater drag coefficient. Consequently, the most effective design for maximising aerodynamic drag is configuration (a), the ``I" frame closed semi-cylinder with the rim segment forming the arch diameter. An additional advantage of (a) is its longer rim segment compared to other configurations, which facilitates easier and more secure attachment of the film.


\section{Experiments} \label{sec_ex}

The simulation studies in \secref{sec_further_sim} offered insights into the aerodynamic drag of various DI designs under uniform wind conditions. Building on these findings, we conducted experiments to validate the theoretical predictions in a practical setting.


\subsection{Experiment Setup}
The experiments are carried out using an existing AeroTail testing platform, designed to simulate a legged robot with an actuated inverted pendulum. The setup includes two adjustable restrictors on the sides, positioned on a T-slot frame segment that also serves as a rail for all components. The actuator is a 3-phase Brushless DC (BLDC) motor, managed by a standard Electronic Speed Controller (ESC). A magnetic encoder attached to the motor allows for accurate measurement of rotation speed.

The AeroTails tested are 3D-printed single units, featuring DI frames at their ends. The central mount of each tail is compatible with the motor on the testing platform. The films for the DIs are crafted from heat shrink polyvinyl chloride bags, shaped using 3D printed moulds, as illustrated in \figref{molds}. The various DI versions, as shown in \figref{di_pics}, are named according to their corresponding shapes in \figref{swshapes}. The numbering after the dot in each name indicates the connection style, with ``1" denoting the tail body situated between the ``T" shapes and ``2" indicating the tail body positioned directly atop them.

The following assumptions are adopted for the aerodynamic modelling and testing conditions used in this study:
\begin{itemize}
    \item The drag coefficient $C_d$ of each drag-inducing (DI) configuration is treated as constant during each test condition, based on CFD and experimental calibration.
    \item The target tail rotational speed during testing is approximately $150$ rad/s, representing the nominal operating condition used for performance comparison.
    \item The Reynolds number setting for CFD study in previous sections are automatically calculated by ANSYS using air density and viscosity (at sea level), characteristic length, and local relative airspeed. These parameters will stay consistent for any calculations done in this section. 
    \item Due to rotational motion of the tail, the local airspeed along the DI span is non-uniform, increasing approximately linearly with radial distance from the rotation axis.
\end{itemize}


\begin{figure}
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/manufacture.jpg}
    \caption{The manufacturing process involves several stages: enclosing the 3D printed mould within the film bag, applying heat using a heat gun; once cooled, the shaped film is cut from the mould; and finally, the film is secured onto the tail frame with adhesive tape.}
    \label{molds}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.99\columnwidth]{pics/DI overview.jpg}
    \caption{The top of the figure shows an AeroTail prototype, detailing all its components. Below are the various DI versions. Red tape is applied to the sides of each DI for enhanced visibility during experiments.}
    \label{di_pics}
\end{figure}

\subsubsection{Experiment 1}
In the first experiment, the objective is to measure the torques generated by different versions of the DI while the device body remains stationary. As you can see in \figref{ex1}, a weight is suspended from the end of the pendulum (aligned with the tail's centre) and placed on a scale. The spring is tensioned, and its linear behaviour implies that an upward force, or torque generated by the tail, will reduce the reading on the scale. The pendulum is fixed at a $30^\circ$ angle, achieved by adjusting the restrictors' height and spacing. For each DI, a specific throttle input is determined to achieve an angular velocity of $150$ rad/s for the tail. The torques generated by each DI version at this constant speed are then measured using the previously mentioned method.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/ex1setup.jpg}
    \caption{The figure on the left illustrates the overall setup for Experiment 1. At the top, a simplified diagram explains the experiment's internal logic, and the procedure is depicted below it.}
    \label{ex1}
\end{figure}

\subsubsection{Experiment 2}
The second experiment aims to assess the balancing capability of the pendulum equipped with various DIs, starting from a predefined angle to an upright position. As illustrated in \figref{ex2}, the pendulum begins at an angle $\theta$ and is halted in the upright position by the restrictor. Throughout this process, the tail's angular velocity is maintained at a constant $150$ rad/s, identical to Experiment 1. The angle $\theta$ is varied incrementally from lower to higher values to determine the maximum angle from which each DI can successfully balance the pendulum.

\begin{figure}[b]
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/Balance_ex.jpg}
    \caption{The sequence of the balancing process in Experiment 2, depicted from left to right. The pendulum is balanced from an angle $\theta$ to the upright position.}
    \label{ex2}
\end{figure}

\subsection{Experiment Results}
\subsubsection{Experiment 1}
The experimental results for torques align with the trends observed in the simulation study (see Video 1). There is a consistent decrease in drag torque from DI (a) to (e), with (d) being the only exception, as seen in the simulations. DIs (a) and (b) share the same frame structure but differ in the orientation of the arch on the film. The simulation concluded that (b) has a lower drag coefficient than (a), thus generating less drag with the same effective area. This is corroborated by the experimental results, with a more pronounced difference than in the simulation. Additionally, the way the DI is connected to the tail (connection style) significantly impacts torque performance. Across all DI versions, style 2 consistently produces better torque results than style 1, with DI a.2 outperforming a.1 by $20.8\%$, solely due to the difference in connection style. The drag coefficients for each DI version, derived from the experimental data, vary significantly from the simulated results. This discrepancy can be attributed to the non-uniformity of drag coefficients along the length of the DIs, and the assumption that drag force acts at the tail's ends when calculating $C_{dE}$. 

\begin{table}
\centering
\caption{Experiment Results}
%\vspace{3mm}
\begin{tabular}{lccccl}
\toprule
\vspace{1mm}
Version & $\tau_{e}$ (mNm) & $C_{dE}$ & $\theta$ & $\theta_{t}$\\
\hline 
a.1 \rule{0pt}{2.6ex} & $24.98$ & $4.14$ & $32^{\circ}$ & $31^{\circ}$ \\
%
a.2 & $29.82$ & $4.95$ & $38^{\circ}$ & $38^{\circ}$ \\
%
b.1 & $20.56$ & $3.41$ & $27^{\circ}$ & $25^{\circ}$ \\ 
%
b.2 & $23.02$ & $3.82$ & $30^{\circ}$ & $28^{\circ}$ \\
%
c.1 & $14.78$ & $2.76$ & $16^{\circ}$ & $18^{\circ}$ \\
%
c.2 & $15.46$ & $2.88$ & $17^{\circ}$ & $19^{\circ}$ \\
%
d.1 & $17.25$ & $3.36$ & $18^{\circ}$ & $21^{\circ}$ \\
%
d.2 & $17.93$ & $3.50$ & $20^{\circ}$ & $22^{\circ}$ \\
%
e.1 & $11.72$ & $2.78$ & $14^{\circ}$ & $14^{\circ}$ \\ 
%
e.2 & $12.49$ & $2.96$ & $15^{\circ}$ & $15^{\circ}$ \\
\bottomrule
\end{tabular}
\label{tab_ex}
% \vspace{-2mm}
\end{table}

\subsubsection{Experiment 2}
The findings from Experiment 2 largely corroborate those from Experiment 1. As can be seen in \tabref{tab_ex} (and Video 1), the experimental results ($\theta$) fit the theoretical value ($\theta_{t}$) calculated from $\tau_{e}$ with a error of $6\%$ on average. A higher torque generated in a stationary position corresponds to better performance in terms of the maximum angle from which the pendulum can be balanced. Notably, the quantitative differences in DI performance are more pronounced at higher angles and less so near the upright position. This variation is primarily due to the trigonometric nature of the gravitational torques resulting from the device's mass. The purpose of Experiment 2 is to validate the conclusions and insights regarding the torque-generating capabilities of the DIs from a stationary setup in a dynamic pendulum-balancing scenario.


\section{Comparative analysis} \label{discuss}

From the results of the experiments came new observations: the numerical difference between the simulation and experimental results cannot be ignored. Torques ($\tau_{s}$) calculated from simulation results, assuming drag force acts at the tail's ends, are shown in \tabref{tab_tau_cal}. They differ largely from the experimental results ($\tau_{e}$ in \tabref{tab_ex}) and the reason behind that was previously assumed in \secref{sec_ex} to be the rotational nature of the tail causing the wind speed to vary along its length. However, this difference was never quantified and thus mathematically proven. If a method could be devised to predict the  experimental results given setup conditions, it would serve as the proof that the differences between $\tau_{s}$ and $\tau_{e}$ are indeed caused by the rotational nature of the tail like previously speculated. In this section, such method will be demonstrated and new findings will be explained accordingly.

\subsection{The analytical method}
To reconcile the discrepancies between the CFD simulation results and experimental outcomes, an analytical method has been devised. This method calculates the theoretical torque each DI shape can produce, considering the distribution of aerodynamic drag along the DI's surface. This approach accounts for the fact that areas further from the tail's centre contribute more significantly to drag generation, due to both higher wind speeds (as drag force is proportional to the square of wind speed) and longer lever arms for torque production. 

\begin{figure}
    \centering
    \includegraphics[width=0.5\columnwidth]{pics/Discussion figure.jpg}
    \caption{Using DI c.2 as an example, an analytical method for calculating torque is derived.}
    \label{dis_fig}
\end{figure}

For DI c.2, the equation to calculate aerodynamic torque for the AeroTail is given by \cite{norby2021enabling}:
\begin{equation}\label{aero_drag_torque}
    \begin{split}
        |\tau_{D}|&=\int_{A} Dl\,dA \\
        & =\int_{L_{1}}^{L_{0}} \frac{1}{2}\rho C_{d}w(\omega l)^2l\,dl
    \end{split}
\end{equation}
where $D$ is the drag force, $l$ represents the distance from a point on the DI to the tail's centre, and $A$ is the frontal area. Variables include $\rho$ (air density), $C_{d}$ (drag coefficient), $w$ (DI width), and $\omega$ (tail's angular velocity). $L_{0}$ and $L_{1}$ are the distances from the tail's centre to the farther and nearer edges of the DI, respectively (as shown in \figref{dis_fig}). The width of the DI, $w$, is not constant but varies along the tail's length. Therefore, the equation can be restructured as:
\begin{equation}\label{aero_drag_torque2}
    |\tau_{D}| = \frac{1}{2}\rho C_{dE}\omega^2\int_{L_{1}}^{L_{0}} l^3w(l)\,dl.
\end{equation}
In this experiment, $\omega$ is a constant at 150 rad/s. The drag coefficient ($C_{dE}$) used in the calculations is taken from measured data in \tabref{tab_ex}. The width function $w(l)$ is expressed using the following equations:
\begin{equation} \label{w_l}
   w(l) = 
   \begin{cases}
   \begin{aligned}
         2(l-B_{1})\sin(\frac{\pi}{2n})+b, \ for \ B_{1} \leq l \leq B_{2} 
   \end{aligned} \\
   \begin{aligned}
        2a\sin(\frac{\pi}{n})-b\sin(\frac{\pi}{2n})-2(&l-B_{2})\sin(\frac{\pi}{2n}), \\
       & for \ B_{2} < l \leq B_{3} 
   \end{aligned}
   \end{cases}
\end{equation}
where $n$ is the number of "T" shapes, with $n=3$ in this case. $B1$, $B2$, and $B3$ define the ranges of $l$ where the width expressions vary. They are determined by:
\begin{equation}\label{b123}
    \begin{split}
        B_{1} & = 5-a \\
        B_{2} & = 5+a\cos(\frac{\pi}{n})-\frac{b}{2}\cos(\frac{\pi}{2n})\\
        B_{3} & = 5+a\cos(\frac{\pi}{n})+\frac{b}{2}\cos(\frac{\pi}{2n}).
    \end{split}
\end{equation}
By integrating \eqref{w_l} and \eqref{b123}, the width equation for the DI is simplified as:
\begin{equation} \label{w_l_2}
    w(l) = 
    \begin{cases}
     1.155l-1.395, \qquad 2.550 \leq l \leq 5.554 \\
     11.43-1.155l, \qquad 5.554 < l \leq 6.896
    \end{cases},
\end{equation}
in which the units are cm. Substituting \eqref{w_l_2} into \eqref{aero_drag_torque2}, the theoretical torque produced by DI c.2 is calculated as $17.78$ mNm. Each DI version has a unique expression for $w(l)$, but for brevity, only the full equations for c.2 are shown here. Using the respective width expressions in \eqref{aero_drag_torque2}, the theoretical torques ($\tau_{a}$) for each DI version can be computed and are summarised in \tabref{tab_tau_cal}. 

\subsection{Discussion for analytical results}

\begin{table}
\centering
\caption{Theoretical Torque Results}
%\vspace{3mm}
\begin{tabular}{lccl}
\toprule
\vspace{1mm}
Version & $\tau_{a}$ (mNm) & $\tau_{s}$ \\
\hline 
a.1  \rule{0pt}{2.6ex} & $28.02$ & $17.71$ \\
%
a.2  & $37.30$ & $17.71$ \\
%
b.1  & $23.08$ & $17.31$ \\ 
%
b.2  & $28.79$ & $17.31$ \\
%
c.1  & $17.14$ & $14.19$\\
%
c.2 & $17.78$ & $14.19$ \\
%
d.1  & $19.75$ & $14.37$\\
%
d.2 & $20.57$ & $14.37$ \\
%
e.1  & $12.60$ & $11.73$\\ 
%
e.2 & $12.27$ & $11.73$ \\
\bottomrule
\end{tabular}
\label{tab_tau_cal}
% \vspace{-2mm}
\end{table}

A visualised comparison of experimental results, analytical answers, and values derived from simulated results is shown in \figref{torque_comp}, with the numerical results shown in \tabref{tab_tau_cal}. It is observed that $\tau_{a}$ and $\tau_{s}$ differ significantly in their magnitudes and $\tau_{e}$ is somewhere between them; compared to results in $\tau_{a}$, those in $\tau_{s}$ do not differentiate between connection styles since there are only five sets of simulated results (one for each DI shape). Since the experimental results ($\tau_{e}$) represent the actual performance, the differences observed in the simulated ($\tau_{s}$) and analytical ($\tau_{a}$) results stem from oversimplified assumptions:
\begin{itemize} 
\item Simulated torques ($\tau_{s}$) are smaller because the simulation assumes a uniform velocity across the DI, neglecting the linear velocity distribution along the tail’s length and the impact of DI orientation. This simplification leads to an underestimation of the aerodynamic torque. 
\item Analytical torques ($\tau_{a}$) are larger because the model uses a constant drag coefficient across the DI, not accounting for the variation in $C_{dE}$ along the length of the DI. This assumption causes an overestimation of the torque compared to the experimental results. 
\end{itemize} 
The analytical results, while not matching the experimental magnitudes precisely, accurately predict the relative differences between the DI versions, a feature the simulated results do not possess. This confirms that the differences between experimental and simulated results, as well as those between the two connection styles, are indeed due to the rotational nature of the AeroTail and the resulting variable wind speeds. A discrepancy is noted in e.1 and e.2, likely due to manufacturing imperfections, as the frontal shape of the DI becomes more circular with increasing $n$, reducing the differences between connection styles.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/torque_compare.jpg}
    \caption{Comparison of experimental results ($\tau_{e}$, orange), analytical answers ($\tau_{a}$, yellow), and values from simulated results ($\tau_{s}$, blue).}
    \label{torque_comp}
\end{figure}

\subsection{An improved analytical method}

A major problem with the results calculated from simulation in using the previous analytical method is that the "average" wind speed is simply taken as the speed at the centre of the DI (end of tail body). This does not reflect reality well since the flow speed is overwhelmingly concentrated on the outside of the DI due to the rotational nature of the tail. Moreover, the drag force is assumed to be acting on the DI centre which produced unrealistic torque because the lever arm is of unrealistic length. Furthermore, since the drag coefficient used for the analytical methods is calculated from simulation results as well, this inaccurate representation of average wind speed undermines the accuracy of the analytic results as well. Therefore, we need a more realistic representation of the average flow velocity. 

\begin{figure}
    \centering
    \includegraphics[width=0.75\linewidth]{pics/Mass_average_velocity.png}
    \caption{Mass-averaged velocity explained. The realistic average wind speed and the lever arm are not at the DI centre like previously assumed.}
    \label{mass_average}
\end{figure}

The average flow velocity in a wind tunnel with varying values can be obtained as the "mass-averaged velocity" (AMV), which takes flow mass into consideration:
\begin{equation}
   \begin{split}
        U_{avg} & = \frac{\int_{A}\rho U d \dot m}{\int_{A} d \dot m} \\
        & = \frac{\int_{A} U^{2} d A}{\int_{A} U d A}
   \end{split}
\end{equation}
where \(U_{avg}\) is the mass-averaged velocity, $U$ is the steady-state flow velocity in the wind tunnel, A is the area of the inlet/plane the flow is coming through, and \(\dot m\) is mass flow rate. Thus, the velocity is weighted in this case. Usually, the \(U_{avg}\) is measured by built-in functions of the CFD; however, if the flow velocity varies linearly in only one axis, it is easy to derive an analytical expression for the mass-averaged velocity. If the gradient $k$ of the varying velocity is known, then the AMV can expressed as:
\begin{equation}
    U_{avg} = \frac{U_{0}^{2}H + U_{0}kH^{2} + \frac{1}{3}k^{2}H^{3}}{U_{0}H + \frac{1}{2}kH^{2}}
\end{equation}
where $U_{0}$ is the velocity at the bottom of the DI, $H$ is the height of the DI and $K$ here corresponds to the angular velocity of the tail. The variables are visualised in \figref{mass_average}. Now the accurate lever arm length $R_{la}$ can be found at:
\begin{equation}
    R_{la} = \frac{U_{avg} - U_{0}}{k}
\end{equation}
from the tail centre. 

Using the newly obtained average velocity and lever arm length, the results are shown in \figref{compare_new} while the original results are in \figref{torque_comp}. We can see the direct calculation results are much closer to the true value for most DI versions. Unlike the old method, the results now differentiate the orientations 1 and 2. However, there are some outliers: DI c.1 and e.1. This is because the direct calculation still ignore the width thus the effective area of the DI, therefore, the analytical methods is still important. Since the drag coefficient $C_{d}$ used in the analytic method is also calculated from the simulation results, more accurate average velocity yields more realistic estimated $C_{d}$ as well. With AMV applied, the average error of the analytical method drops from $14.4\%$ to $6.6\%$ and the large gaps in previous studies cease to exist. 

\begin{figure}
    \centering
    \includegraphics[width=0.75\columnwidth]{pics/DI_analytic_improved.jpg}
    \caption{The new estimated results compared to experimental results. Here, the mass-averaged velocity is applied to both the direct calculation and analytical method. The results are way closer to the true (experimental) results.}
    \label{compare_new}
\end{figure}



\section{Chapter summary} \label{sec_con}

This chapter presents a comprehensive study on the optimisation of drag-inducing device designs for the AeroTail, incorporating CFD simulations and experimental validation. The primary goal is to maximise aerodynamically induced torque for robotic balancing, adhering to specific dimensional and weight constraints. Given the scarcity of analytic and experimental research on optimal shapes for DIs, CFD simulations emerge as a foundational step for this investigation. The simulations reveal that among the evaluated shapes, the semi-sphere demonstrates superior performance. It is also noted that DIs with slender frontal surfaces, arch-like profiles, and closed ends generally yield higher drag coefficients. However, it's established that a complete frame design results in material wastage, leading to the introduction of the more efficient ``T" shaped frame. Subsequently, the DI frame design is optimised to maximise frontal area, which is a more reliable attribute to enhance within the current frame structure. While a strong theoretical foundation is provided by the simulation studies, these insights cannot be validated without experimental trials. In the following experiments, the results confirmed that the most torque is generated by DI with a semi-cylindrical film of configuration (a) and connection style 2, outperforming style 1 by $20.8\%$. This indicates that connection style plays a huge role in the difference of torque-generating performances. Subsequently, an analytical method developed for torque calculation is used to elucidate the disparities in quantitative experimental outcomes.

This chapter sees the end of the balancer side of the research. In the next chapter, the C-SEA jumping mechanism are explained. Both theoretical simulation and experimental validations on the findings and performance will be included. 

%Looking ahead, future improvements to the optimisation problem could consider the varying wind speeds along the tail's length. The method for calculating theoretical torque, as discussed in \secref{discuss}, demonstrates that the wind speed variance across different parts of the DI, due to the rotational nature of the tail, significantly affects torque output. Incorporating this factor into the optimisation problem--especially prioritising the concentration of effective area towards the tail's outer regions--could yield even more definitive solutions.
