A Neural Network Approach to Estimate Transient Aerodynamic Properties of a Flapping Wing System

Abstract

Understanding the causal impacts among various parameters is essential for designing micro aerial vehicles (MAVs). The simulation of computational fluid dynamics (CFD) provides us with a technique to calculate aerodynamic forces precisely. However, even a single result regularly takes considerable computational time. Machine learning, due to the advance in computer hardware, shows another approach that can speed up the analysis process. In this study, we introduce an artificial neural network (ANN) framework to predict the transient aerodynamic forces and the corresponding energy consumption. Instead of considering the whole transient changes of each parameter as inputs, we utilised the technique of Fourier transform to simplify the ANN structure for minimising the computation cost. Furthermore, two typical activation functions, rectified linear unit (ReLU) and sigmoid, were attempted to build the network. The validity of the method was further examined by comparing it with CFD simulation. The result shows that both functions are able to provide highly accurate estimations that can be implemented for model construction under this framework. Consequently, this novel approach makes it possible to reduce the complexity of analysis, study the flapping wing aerodynamics and enable a more efficient way to optimise parameters.

1. Introduction

While a human can fly into the sky with a machine, the mechanism of insect flight remains yet a mystery of sorts. Unlike conventional artificial aircraft, an insect exhibits its fascinating aerial manoeuvrability by repeatedly flapping its wings. This particular mechanism has recently been extensively investigated to develop an improved micro aerial vehicle (MAV). Furthermore, aerodynamics at a small Reynolds number provides a more efficient flight [1], which allows a MAV to cruise at a low speed to execute examination tasks [2]. As MAVs can overcome terrain constraints, they are expected to search for victims in narrow buildings or explore dangerous environments by employing various sensors [3,4]. However, as the flapping wing system is a relatively novel concept compared with other aircraft, its mechanism has not been fully revealed yet. Considerable time is therefore required to examine the impact of various variables.

Among various methods, some studies have reported their findings through biological observations. Ellington [5] claimed that wing paths had no consistent patterns among numerous insects. Wakeling and Ellington [6] displayed exceptional steady-state aerodynamic property of dragonfly wings and utilised it to predict the parasite drag. Josephson and Stevenson [7] measured the oxygen consumption from insects to evaluate the energy efficiency of various flight patterns; Dial et al. [8] also presented the measured power consumption of birds that flew at different speeds by examining the electromyograms (EMGs).

As it is tough to reproduce the same experiment due to the individual differences and the uncontrolled environmental variables, some studies consequently built flapping wing mechanisms to clarify the relationship between various parameters. By utilising a rack-pinion mechanism, T.A. Nguten et al. [9] investigated the effect of parameters such as wing aspect ratio and flapping frequency; Sato et al. [10] built a flapping wing robot and discussed the control strategy of altering the flapping amplitude of wings. Miyoshi et al. [11] found that the asymmetric wingbeat amplitude only affects the pitch and roll moments. These mechanism-based approaches enable us to regenerate the same flying condition readily. Furthermore, by controlling the variables, we can figure out the impact of each parameter on flight performance.

On the other hand, without utilising the technique of particle image velocimetry (PIV), aerodynamic information such as the leading edge vortex (LEV), which produces high transient lift [12], can be directly visualised by computational approaches. From the simulation, we can further obtain precise details such as vorticity to explain the mechanism behind it [13]. Zou et al. [14] and Lai et al. [15] unveiled impacts caused by phase lag and the wing–wing interaction through computational fluid dynamics (CFD) simulation; Johansson and Henningsson [16] compared the effect of the clap mechanism between rigid and flexible wings. This numerical approach allows researchers to simulate how the airflow interacts with various objects without spending extensive time and effort creating experimental environments [17,18,19]. To obtain precise aerodynamic outputs, building a fine mesh of calculation field is necessary. Nevertheless, it takes several days or weeks for a machine to complete the calculation.

With the advances in computer hardware, machine learning shows another approach for modelling. Unlike conventional computation, related mechanisms are not required (e.g., governing equations). Recently, some studies have shown that machine learning methods can be implemented in aerospace science, including flight pattern recognition [20] and aeroelastic estimation [21]. In the study of flapping mechanics, some researchers have tried to adopt this new technique into control systems [22,23,24]. Nevertheless, the application for flapping flight aerodynamic analysis has not yet been fully developed. Two studies [25,26] utilised this method to predict the net forces produced in a single flapping cycle, but the transient analysis has not been completed yet. Therefore, researchers still relied on conventional approaches to obtain details, such as the changes during upstroke and downstroke, to explain the impacts of different flight modes.

In this study, we introduce a neural network approach accompanied by CFD simulation to shorten the considerable computational time of transient analysis. We first collected flapping kinematics by biological observations and utilised them to build a model for CFD simulations. Afterwards, the data were split into training, validation and testing groups. We utilised the first two groups of data to create neural network models and evaluated them with the testing data. As reports [27,28] have shown that wing and body kinematics are the main parameters affecting butterfly flight performance, we considered wing rotation and body oscillation as inputs and utilised them to predict the corresponding changes in aerodynamic forces and power consumption. The effectiveness of the method was examined by comparing it with the results obtained by CFD simulation afterwards.

2. Data Collection

2.1. Biological Experiment

To record the flapping motion, we utilised the blue tiger butterfly (Tirumala septentrionis) as a reference (Figure 1a). Before the measurement, objects were frozen at −7 °C for 24 h. The morphological parameters were then measured ($N=4$), including mass, body length and wing area, span and mean chord (see

Table 1. Measured ($N=4$) and simulation model’s morphological parameters of the butterfly.

Parameters	Measurement	Simulation
Mass (mg)	$379.93\pm 72.99$	350
Wing Area (mm${}^2$)	$923.67\pm 42.55$	925
Wingspan (mm)	$46.00\pm 1.23$	46.04
Mean Chord (mm)	$20.07\pm 0.39$	20.11
Aspect Ratio	$2.29\pm 0.02$	2.29

). Afterwards, we calculated the aspect ratio $AR$ of a single wing with the equation:

$$AR=\frac{S}{\overline{c}},$$

(1)

in which S and $\overline{c}$ represent the wingspan and mean chord, respectively. These parameters were considered a reference to build a simulation model later.

To capture the details of a butterfly’s body action, we utilised two high-speed cameras (Phantom v7.3 and v310, Vision Research, Wayne, NJ, USA). Both cameras were operated at 1000 frames per second with a resolution of 1024 × 1024 and placed outside a transparent chamber (size: 1 m × 0.35 m × 0.35 m) in which a butterfly could fly inside freely. As these cameras were placed orthogonally, we defined the direction from the chamber to the two cameras as y- and z-axis (see Figure 1b); the synchronised photographed films were examined to determine the angles between various body parts. To attract a butterfly to fly forward, we placed a light source on one side of the chamber. As the study focused on the flapping motion of forwarding flight, we abandoned the data for which the butterfly flew with turning. On counting the number of frames, we deduced the wingbeat frequency to be around 11.020 ± 1.076 Hz ($N=15$).

To obtain the flapping kinetic equations, we utilised the coordinates measured from the experiments to rebuild them. Figure 2a shows the heading direction $\overrightarrow{b}$ was determined by the body from $B_0$ to $B_1$. The angle between the x-axis and $\overrightarrow{b}$ was defined as the body pitching angle $\theta$. Still, the vector $\overrightarrow{w_1}$ was determined by the root $W_0$ and the wingtip $W_1$; the sweeping angle $\eta$ was obtained by calculating the complementary angle of the angle between $\overrightarrow{w_1}$ and $\overrightarrow{b}$. Additionally, the wing plane was constructed with $\overrightarrow{w_1}$ and $\overrightarrow{w_2}$ (the direction pointing from $W_0$ to $W_2$, with $W_2$ located on the trailing edge of the hindwing). Lastly, the wing rotation angle $\alpha$ was computed from the normal vector of the wing $\overrightarrow{w_N}$ (i.e., the direction of $\overrightarrow{w_2}\times \overrightarrow{w_1}$) and $\overrightarrow{r}\times \overrightarrow{b}$; the angle between the latter and the negative z-axis was defined as its flapping angle, where $\overrightarrow{r}$ represents the axis of rotation. These parameters were utilised to describe the flapping flight behaviour of a butterfly.

Figure 2b shows the variation of the four angles in a flapping period measured from the experiments ($N=15$) with 95% confidence intervals (shaded areas). These curves were recognised as input parameters for subsequent simulation afterwards.

2.2. Computational Fluid Dynamics Simulation

To deduce the aerodynamic interactions that were hard to observe directly, we chose CFD simulation to generate highly accurate datasets. The obtained morphological parameters and kinematic equations were further utilised to build a butterfly model and regenerate the flight behaviour through a commercial solver (Fluent, Ansys, Canonsburg, PA, USA). The SIMPLEC algorithm was applied to resolve pressure and velocity fields [29]. To reduce the computational cost, we simulated the flow field under a relative condition frame. Therefore, the butterfly was flying at the centre of a sphere (see Figure 3a). The butterfly would thus encounter an incoming airflow with a virtual acceleration a created by its flight. To avoid inaccurate outcomes affected by the wall effect, the diameter of the spherical flow domain was set to $20S$ (around 900 mm). As the Reynolds number of a natural flying butterfly is around ${10}^3$–${10}^4$[30], the medium was considered an incompressible Newtonian laminar airflow, with a density of $\rho =1.23$ kg/m${}^3$ and viscosity of $\mu =1.79\times {10}^{-5}$ Pa·s. Furthermore, we utilised the following governing equations for computation:

$$\begin{array}{c}\nabla ·\mathbf{u}=0\end{array}$$

(2)

$$\begin{array}{c}\rho {\displaystyle \frac{D\mathbf{u}}{Dt}}=-\nabla p+\mu {\nabla }^2\mathbf{u}+\rho (g+a)\end{array}$$

(3)

in which u, t, p and g represent the velocity, time, pressure and gravity, respectively. The entire domain was split into two for giving different boundary conditions [29,31,32]. The front of the butterfly was considered to be the velocity inlet with an incoming flux accompanied by a. As a result, the stream fluctuated with the flight at each time step. The rear sphere was the pressure outlet. These conditions were defined as:

$$\begin{array}{c}\mathbf{u}=-\widehat{a_x}\dot{x}-\widehat{a_y}\dot{y}\end{array}$$

(4)

$$\begin{array}{c}p=-\rho gy,\end{array}$$

(5)

in which $\widehat{a_x}$ and $\widehat{a_y}$ indicate the unit vector along the x- and y-axis, respectively.

The grid size of the butterfly’s surface was 0.5 mm with a no-split condition (the air had zero velocity relative to the boundary), whereas its surroundings were set to 1 mm. This setting can enhance the precision of the simulation results [31,33]. We also utilised the lift force to do the grid convergence test. The result converged as the number of grids increased. Figure 3b illustrates the comparison among three different settings: the coarse grid (blue solid line, 7 million), the medium grid (red dashed line, 8 million) and the fine grid (amber dotted line, 12 million). From the result, we found that the maximum value difference between the coarse and fine settings was about 0.5% merely. Considering the balance between accuracy and computational cost, we hence selected the medium setting for the following computation.

We also chose the method of dynamic mesh (smoothing and remeshing) to prevent a negative mesh volume [34]. A flapping cycle was divided into 250 calculation time steps. The simulation outputs were the horizontal force $F_h$, vertical force $F_v$, normal force acting on a single wing $F_w$ and power consumption P at the 10^th flapping cycle (stable flight). All these values were nondimensionalised by the following equation:

$$\left[\begin{array}{c}{F}_h\\ F_v\\ F_w\\ P\end{array}\right]={\displaystyle \frac{1}{{\displaystyle \frac{1}{2}}\rho V^2\left(2S\overline{c}\right)}}\left[\begin{array}{c}{f}_h^{*}\\ f_v^{*}\\ f_w^{*}\\ p^{*}/V\end{array}\right],$$

(6)

in which the wingtip velocity $V=4\varphi fS$; $f_h^{*},f_v^{*},f_w^{*}$ and $p^{*}$ were the corresponding untransformed forces and power. The Reynolds number of the simulation result was 6130, which was close to the experimental result of 6050.

As the study focused on the influence caused by amplitudes of body oscillation and wing rotation, we altered the simulated conditions of body pitching angle and wing rotation angle by multiplying them with scalers b and w. In total, 25 flapping data were collected (see Figure 4).

3. Artificial Neural Network

An artificial neural network (ANN) is constructed by several connected computing cells (artificial neurons) inspired by the human brain as shown in Figure 5. The output y of the cell n can be calculated by the following equation:

$$y_n=\sigma (\sum _{i=1}^m{w}_{ni}{x}_{ni}+b_n),$$

(7)

in which $x_{n1},x_{n2},...,x_{nm}$ are the input signals; $w_{n1},w_{n2},$...$,w_{nm}$ are the respective weights; $b_n$ is the bias; and $\sigma$ is the activation function (transfer function) [35]. The inputs were the scalers b and w. To increased the training performance, we normalised the values in the interval of $[0,1]$ [36]. On the other hand, instead of predicting 250 time steps for each aerodynamic property, we simplified it by fitting each parameter by a Fourier series:

$$F\left(t\right)=A_0+\sum _{n=1}^m{A}_{n}cos2\pi nft+B_{n}sin2\pi nft,$$

(8)

where m is the order of the Fourier series. When choosing $m=5$ to fit the data, each aerodynamic property could be described by a vector formed by 11 parameters. We hence considered these 44 parameters (four aerodynamic properties, $F_h,F_v,F_w$ and P) as the output of an ANN model. As the number of outputs was shrunk from 1000 to 44, we could consider utilising a smaller number of hidden cells for the following calculation. The training process hence can be curtailed at the same time.

Figure 6 shows two common activation functions utilised for neural networks. Recently, Rectified Linear Unit (ReLU) has been chosen as the activation function in various applications [37]. The function is defined by the formula:

$$ReLU\left(x\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}x<0,\hfill \\ x\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(9)

Although the vanishing gradients problem is a critical issue when training an ANN, ReLU has a constant gradient of 1 when the input is greater than 0. Therefore, it is widely applied to ANNs. However, previous flapping wing studies chose the sigmoid function to estimate the mean values of aerodynamic coefficients [25,38]. The function can be described by:

$$Sigmoid\left(x\right)=\frac{1}{1+e^{-x}}.$$

(10)

As the s-shaped output of the sigmoid varies continuously in the interval of $[0,1]$, the activation values hence do not disappear. Because we aimed to predict transient aerodynamic statuses rather than their mean values, we need to evaluate which function can provide precise estimations. We consequently constructed two models based on these functions and compared the results afterwards.

To train the model, we utilised Adam optimiser [39], an advanced backpropagation method [40], with the loss function of mean square error (MSE) to minimise the error. Although many studies [41,42,43,44,45] had proposed their rules of thumb to select the number of hidden units, the main idea is to improve the prediction accuracy and minimise errors. Therefore, we randomly picked up 80% of the data to train two types of models through a various number of hidden units iteratively. Among these datasets, 80% (i.e., 64% of the entire data) were utilised for tuning weights and bias; 20% (i.e., 16% of the entire data) were for validation. Considering the balance between accuracy and computational time, we have examined the loss when the number of hidden neurons was between 20 to 40. To avoid the setting merely benefiting specific cases, each training was repeated 30 iterations with the same settings but different training data selections. Figure 4 shows one of the random states that we have utilised. The learning rate and the epoch size were set to 0.001 and 60,000.

Figure 7 depicts the mean convergent curves of each case during the training process. From the testing result, we found that though the training loss of ReLU model kept decreasing when epoch size was greater than 10,000, the validation loss of it remained around ${10}^{-3}$. We further found that as the number of hidden units was 40, which had the smallest loss, its validation loss was $8.28\times {10}^{-4}$ and $8.05\times {10}^{-4}$ when epoch sizes reached 30,000 and 60,000. As the difference between these two values was less than 3%, we considered the 40 hidden neurons and 30,000 epochs to be an adequate setup. On the other hand, the training and validation loss of sigmoid model plots reveal that a larger number of hidden neurons did not promise a smaller loss. When the number of hidden units was 30, the validation loss difference between epochs 50,000 and 60,000 was less than 3% (epoch 50,000: $8.12\times {10}^{-6}$, epoch 60,000: $7.93\times {10}^{-6}$). Consequently, we considered the number of hidden units and epochs to be 30 and 50,000, respectively.

It was unsure if the accuracy of the ReLU activation function was sufficient to predict the four aerodynamic properties, though the validation loss of the sigmoid-based model was even lower. We hence utilised both models for the following analysis. To avoid overfitting, the unused 20% of the data were assigned as testing datasets. As all the data generated by the CFD simulation were from the same system, the model that fits one group should fit the other as well. Because these cases were not utilised for previous training, they were considered to be unseen data for the following model evaluation. If the number of datasets is sufficient, we can utilise a large number of cases to evaluate our models. However, it is not easy to collect hundreds of results as it generally takes about 4 to 7 days to complete a single simulation. To deal with a limited amount of data, k-fold cross validation was implemented for obtaining a reliable evaluation [46]. To implement this method, the database needs to be randomly split into k groups first so that we can evaluate a model k times. In each iteration, we only utilise $k-1$ groups to train our network. As there is one group that does not join the training process, it can be utilised for testing. Therefore, this can be utilised to protect our models from overfitting as the model will be examined by k different combinations. To keep 20% of the entire data for testing, we chose k to be 5 in this study (see Figure 8).

4. Results and Discussion

4.1. Model Comparison

As the outputs of our models were coefficients of the four aerodynamic parameters, we converted the signal back by utilising Equation (8). The prediction can hence be compared with the original curves by calculating the coefficients of determination $R^2$. The mean and standard deviation (SD) of the two models are given in

Table 2. Statistics of the ReLU- and sigmoid-based models (mean ± SD).

Parameters	ReLU		Sigmoid
	${\mathit{R}}_{\mathit{Train}}^{\mathbf{2}}$	${\mathit{R}}_{\mathit{Test}}^{\mathbf{2}}$	${\mathit{R}}_{\mathit{Train}}^{\mathbf{2}}$	${\mathit{R}}_{\mathit{Test}}^{\mathbf{2}}$
${\mathit{F}}_{\mathit{h}}$	$0.995\pm 0.018$	$0.989\pm 0.025$	$0.999\pm 0.003$	$0.998\pm 0.005$
${\mathit{F}}_{\mathit{v}}$	$0.998\pm 0.007$	$0.995\pm 0.012$	$0.999\pm 0.003$	$0.999\pm 0.003$
${\mathit{F}}_{\mathit{w}}$	$0.999\pm 0.004$	$0.997\pm 0.006$	$0.999\pm 0.002$	$0.999\pm 0.002$
$\mathit{P}$	$0.996\pm 0.016$	$0.989\pm 0.034$	$0.999\pm 0.003$	$0.997\pm 0.009$

. We first checked if the overtraining occurred on our model. Since the training and testing $R^2$ values are all close to 1, the models encountered neither underfitting nor overfitting problems. Although the MSE loss of the sigmoid-based model was smaller than the other one, as they both have mean values greater than 0.99, the differences were not obvious. However, the SD values of the ReLU-based model were slightly higher. The result implies that the sigmoid activation function generally provided higher precise estimations.

As it is difficult to understand the difference between the two models by just viewing the $R^2$ values, we picked up a case from the testing group, which had a maximum $R_{test}^2$ value difference of 0.005, for comparison. Figure 9 presents the output aerodynamic properties as the b and w were assigned to 1 and 0, respectively. By comparing these curves, we can find that the red dashed line (sigmoid) almost overlaps with the black solid line (CFD simulation). On the other hand, the blue dashed line (ReLU) is slightly off in the peaks. By way of illustration, the blue dashed line shifted marginally to the right at $T=0.3$ in the $F_h$-T plot. Nevertheless, the trends of four aerodynamic properties can still be identified.

4.2. Aerodynamic Performance

While the traditional CFD simulation takes several days to complete a single computation, these two ANN models take less than a second merely. Due to the high $R^2$, both methods achieve elevated estimation accuracy and can be implemented to accelerate the analysis process. The networks hence can be utilised to investigate the interactions between body oscillation and wing rotation behaviours. Owing to the limited payload, a MAV cannot carry a battery of substantial volume. Consequently, how to efficiently generate lift is a critical issue. We can hence utilise this neural network approach to optimise the flapping motion, such as searching the efficient flying kinematics.

In our study, the mean lift efficiency ${\overline{E}}_v$ can be defined by:

$${\overline{E}}_v=\frac{{\overline{F}}_v}{\overline{P}},$$

(11)

where ${\overline{F}}_v$ and $\overline{P}$ are the mean values of $F_v$ and P in a single flapping period, respectively. We could hence obtain a quick result by utilising the ANN method. Figure 10 depicts the corresponding ${\overline{F}}_v$ and ${\overline{E}}_v$ when b and w are in the interval of [0, 2]. Because the two models presented similar outcomes, we only display the sigmoid-based calculation in Figure 10. The result illustrates that the maximum ${\overline{F}}_v$ appears at $\left[\begin{array}{cc}b,& w\end{array}\right]=\left[\begin{array}{cc}2,& 2\end{array}\right]$. Nevertheless, the lift efficiency can still be improved by reducing w to 1.2. To verify if it fits the actual circumstance, we ran the simulation and compared the results (see Figure 11). The MSEs of $F_v$ and P were $5.19\times {10}^{-4}$ and $8.35\times {10}^{-4}$, respectively. Considering the high $R^2$ and low MSE values, we believe this neural network approach can provide precise estimations which can be utilised for studying the flapping wing system.

Although it is still challenging to provide details such as flow fields with the current framework, we can utilise it to make a quick analysis and perform full CFD simulations of specific cases to clarify abnormal phenomena. This technique hence will be a great advantage when dealing with a complex system. Moreover, as we verified that the network could successfully estimate transient aerodynamic properties, we can extend the framework to include more complicated flapping motions as inputs in the future.

5. Conclusions

In this study, we introduce a novel neural network approach to speed up the transient analysis of flight mechanics. For evaluation, we analysed the butterfly’s flapping motions by CFD simulation and trained the model with these datasets. To simplify the model, we further utilised Fourier transform to reduce the number of neural cells. Through a series of tests, we found that both ReLU- and sigmoid-based models can accurately predict these coefficients, which can be utilised to obtain the original transient results. This enables us to rapidly estimate the corresponding aerodynamic properties with the given inputs.

The series of work conducted in this study aims to reduce the computational time cost. As the complex structure of flapping motion aerodynamics is non-linear, conventional approaches take excessive effort to analyse the unsteady system. With the aid of neural networks, we do not need to stick to CFD simulation for each case but can still obtain precise results without spending plenty of time. The technique reported here sheds new light on the development of flapping wing systems. This is a great advantage when a large amount of computing is required and provides us with a more efficient way to discuss the interactions among various parameters. Prior to this study, it was difficult to verify that the selected parameters were optimal. The technique introduced in the study makes it possible. We believe the framework will provide an efficient way to delve deeper into the flight mechanism and design a more efficient MAV.