In-Flight Solar Radiation Intensity Measurement Using A Small Unmanned Aerial Vehicle

Up to this time, particularly in recentyears, numerous trials have been conducted to develop small solar-powered UAVs. However, only a few prototypes realized the perpetual flight. A research team from Switzerland performed one of the most advanced research studies. Noth [1–3] successfully developed the Sky-Sailor aircraft. Oettershagen made essential contributions to the field of small solar-powered UAVs, which led to the development of AtlantikSolar [4,5]. In [6,7], he presented a set of flight data logs from perpetual flights. In [8], he validated the proposed solar power model. Dwivedi et al. published a series of papers [9–12] describing MARAAL aircraft. Dantsker etal. [13,14] developed Solar Flyer UAV and presented test results from long-endurance flights. Weider, with his team [15] developed SunSailor UAV. In Poland, there are only a few examples of research concentrated on this technology. Among them, the AGH Solar Plane [16] is notable.

To design the solar-powered aircraft correctly, the amount of solar radiation that can be collected must be precisely estimated. The UAV can perform complex maneuvers, and the solar cells can change the orientation with respect to the direction of the incoming radiation. Up to this time, several researchers have addressed the need to develop an accurate energy harvesting model for solar-powered UAVs [17–20]. Li et al. [21] presented a simulation that could predict the amount of collected solar energy and was partially validated by flight trials. Hrovatin and Žemva [22] presented and validated a model in a flight test using a small UAV named Bramor ppX. Rajendran and Smith [23] also proposed a model that can be used for predicting solar radiation and sizing UAVs. Morton and Papanikolopoulos [24], as well as Sineglazov and Kara-betsky [25], proposed a simple approach, suggesting that the captured solar power could be modeled using a sinusoid. However, the presented results were not validated against the experimental data in sufficient detail.

It is essential to note that creating reliable energy models for flying platforms requires considering several key phenomena. Such a model should be integrated, at the very least, with the kinematic description of motion. Often, the simple point mass models of the aircraft are adopted [26,27]. The geometry of the system must be taken into account. Zhang et al. [28] included each wing segment’s dihedral in the model’s bank angle. Ji et al. [29] investigated the mismatch loss problem in detail. Wu et al. [30] considered the equipping of a Z-shaped wing in solar cells. Most of the models are developed under the assumption of a clear sky [31]. For example, Huang et al. [32] reported the energy harvesting model and simulated the results for several cities in China, but again without experimental verification. However, neglecting clouds could significantly overestimate the total collected energy. For HALE aircraft, the reflected component could be neglected [33], but for small UAVs flying at low altitudes, this issue is not apparent. Gao et al. [34] suggested that for UAVs operating at very high altitudes, the influence of the roll and pitch angles on the amount of collected energy can be neglected. This assumption seems to be very rough. Martinez et al. [35] concluded that the roll angle can be ignored for predicting the amount of solar energy that UAVs can collect when the Sun is at high elevation, but not in the morning and evening (when the Sun is low above the horizon). Lee and Yu [36] considered the atmospheric transmittance. A detailed energy model for solar-powered UAVs was described by Wu et al. [37,38]. Brizon [39] presented a solar energy model for HALE platforms. He investigated the influence of several factors (latitude, longitude, and trajectory type) on the amount of collected energy. Meyer et al. [40] presented a simplified approach for predicting the amount of solar energy that UAVs can collect, but the method was not experimentally validated. Park et al. [41] developed a ground-based virtual flight system that mimics the real UAV’s dynamics and compared results with those from simulation.

The technology of LALE UAVs is not fully mature. Many of the aforementioned models were not sufficiently validated. To make the model realistic and practically usable, its parameters must be carefully tuned.

The primary contribution of this research is the description of solar radiation measurements evaluated in flight using a small UAV. The developed equipment is small, lightweight, and low-cost, making it suitable for adoption on other platforms. The aircraft and solar radiation models were developed and implemented in MATLAB/SIMULINK. These models were validated with the obtained data. The developed simulation might be used to design photovoltaic installations on solar-powered aircraft.

The structure of the paper is as follows. Chapter 2.1 describes the test platform. The solar radiation measurement system was presented in part 2.2. Section 2.3 shows mathematical models of the aircraft and solar radiation models. In part 2.4, the ground measurements of the power consumption were reported. Chapter 3 presents the results of the flight tests and model validation. The paper concludes with a summary of the main findings and their implications. Finally, further research directions are suggested.

The Multiplex Funcub NG RC airplane was used as a test platform (Fig. 3). The aircraft features a traditional design, characterized by a high-mounted wing, single engine, fixed landing gear, and a tailwheel configuration. It is characterized by ease of use and is frequently used as a training platform. The airplane has been modified in comparison to its original equipment.

Figure 3.

Test platform used in the experiments

Takeoff mass was m = 1,96 kg, wing area was S_w = 0,399 m², span b_w = 1,41 m, chord c_w = 0,23 m, and overall length 1,050 m. Moments of inertia were obtained experimentally using a bifilar pendulum and are equal: I_xx = 0,07598 kg·m², I_yy = 0,09504 kg·m², I_zz = 0,12134 kg·m². The product of inertia I_xz = 0,00259 kg·m² was calculated using the CAD model. The airplane’s standard equipment consists of a ROXXY BL C35-42-930 KV electric motor, a ROXXY BL Control 740S-BEC electronic speed controller, an APC 12x6E propeller, and a set of Hitec HS 55+ andHS 65HB servos. Power is provided by a single 3S package of 3200mAh capacity produced by Dualsky. The airplane is controlled using a Futaba T14SG radio controller and a Futaba R7008SB receiver module. Originally, the model weighed 1,380 kg with a wing loading of 3,46 kg/m². Modifications have been performed to adapt the aircraft for test platform functionality. Additional equipment (equivalent to 42% of its initial weight) caused an increase in platform weight to 1,960 kg, resulting in a wing loading increase up to 5,66 kg/m². Weight change is naturally connected with the reduction of maximum flight duration.

Among additional equipment, it is possible to distinguish the following hardware: the Cube Orange, Here2 GPS receiver, airspeed sensor, ultra-long range radio modem RFD 868x, FPV camera with vision transmitter AV TS832, 2 x Adafruit BH1750 light sensors, Arduino microcontroller with Datalogger Shield V1 equipped with RTC DS1307 and SD memory card slot, 9V backup battery for powering the Arduino and additional ballast needed to balance the test platform properly To reduce unpredictable platform behavior, it has been decided to disable the flaps operation.

The solar radiation measurement system was developed in cooperation with the School of Engineering and Built Environment (and Aviation) at Griffith University (in Brisbane, Australia).

Typically, the solar irradiance is measured using pyranometers or photovoltaic reference cells. However, such sensors are expensive and large, making it problematic to install them on small UAVs, such as Funcub. Due to this reason, it was necessary to develop a customized, relatively low-cost solution that could serve as an alternative to high-grade sensors.

The system consists of the following components: two identical light sensors BH1750 (Fig. 4) located at the end of the wings (at the top surface), 3D printed sensor holders mounted with Velcro (Fig. 5), Arduino microcontroller, DataLogger shield V1.0 RTC DS1307 with SD card reader dedicated for Arduino, 32GB microSD card, 9V rechargeable lithium-ion battery and connection cables.

Figure 4.

BH1750 light sensor

Figure 5.

Sensor mounted on the surface of the left wing

The data recorder is based on the Arduino board (Fig. 6).

Figure 6.

Completed set of solar radiation measurement system

The measurement system has been designed to be completely independent of autopilot hardware. To synchronize data registered by the autopilot and the developed system, a signal responsible for aileron deflection has been recorded to a microSD card.

Measurement system functionality has been obtained by programming the Arduino microcontroller using dedicated libraries. The data acquisition system communicates with light sensors and saves recorded information to the microSD card.

The first stage of the communication protocol checks was completed using two microcontrollers (Fig. 7).

Figure 7.

Microcontroller protocol communication test

Integration of the light sensor required setting individual IP addresses for each device. In further steps, a complete system was assembled on the stationary test bench (Fig. 8a).

Figure 8.

(a) Stationary test stand during experiments in Brisbane (b) Solar Survey 100 irradiance meter

The developed system estimates solar irradiance (a radiometric quantity) using information about illuminance (a photometric quantity).

The BH1750 light sensor outputs the illuminance (expressed in lux) directly. A set of preliminary measurements has been completed, providing that the sensor’s reading has been saturated at around 32000 lux. However, in outdoor conditions at noon, the illuminance can achieve more than 110000 lux. For this reason, it was necessary to modify the default setup and increase the sensor’s reading range to the required value. Improvements have been made with the change of device resolution, focusing on software code modification. Next, the procedure included a light sensor calibration from lux to power expressed as W/m² delivered by the Sun to Earth’s surface, which was accomplished using the Solar Survey 100 irradiance meter (Fig. 8b). Finally, data were logged on the microSD at a frequency of 10 Hz. Additionally, the data stream was synchronized in time with the selected single reference signal from the autopilot (right aileron deflection angle). This allowed for the simultaneous analysis of the flight logs (from the autopilot) and the solar radiation measurement system.

The system successfully completed the stationary tests. It has been integrated with the RC airplane research platform in the next step. Sensors were mounted on the upper wing surface. The distance between devices is equal to 1,16 m. Additionally, the surfaces of both sensors were aligned with respect to each other. Minor mounting inaccuracies might influence the quality of the results. Nevertheless, this method was deemed sufficient for testing purposes. Because the airplane wing has a small inclination, it is expected to obtain slightly different readings from the light intensity sensor, as its vectors are not parallel.

Due to the fact that wires and sensor holders are located on the wing’s upper surface and are not integrated into the structure, the airplane’s aerodynamic characteristics may be slightly influenced, causing increased drag force.

A detailed mathematical model of the aircraft was developed to predict the amount of solar energy that can be collected during the flight.

2.3.1.

Equations of Motion

It was assumed that the aircraft is a rigid body with 6 degrees of freedom. The coordinate systems used in the model are presented in Fig. 9. The navigation frame O_nx_ny_nz_n is used to describe the position of the aircraft. Origin O_b of the body-fixed frame O_bx_by_bz_b is attached to the center of mass of the plane. Additionally, O_gx_gy_gz_g was used as the intermediate frame (axes of the O_gx_gy_gz_g are parrarel to O_nx_ny_nz_n). Moreover, to define the aerodynamic properties wind frame O_wx_wy_wz_w and wind-laboratory frame O_wlx_wly_wlz_wl were used.

Figure 9.

Coordinate systems used in the model

The equations of translational motion are as follows: 1[ U˙V˙W˙ ]=1m[ XbYbZb ]−[ 0−RQR0−P−QP0 ][ UVW ]\left[ {\matrix{ {\dot U} \cr {\dot V} \cr {\dot W} \cr } } \right] = {1 \over m}\left[ {\matrix{ {{X^b}} \hfill \cr {{Y^b}} \hfill \cr {{Z^b}} \hfill \cr } } \right] - \left[ {\matrix{ 0 & { - R} & Q \cr R & 0 & { - P} \cr { - Q} & P & 0 \cr } } \right]\left[ {\matrix{ U \cr V \cr W \cr } } \right] where U, V, W – components of the velocity vector (in body axes), P, Q, R – roll, pitch, and yaw angular rates, X^b, Y^b, Z^b – axial, side, and normal forces (the upper-right superscript defines the coordinate system in which the quantity is expressed). The equations for rotational motion are: 2[ P˙Q˙R˙ ]=I−1[ LbMbNb ]−I−1[ 0−RQR0−P−QP0 ]I[ PQR ]\left[ {\matrix{ {\dot P} \cr {\dot Q} \cr {\dot R} \cr } } \right] = {{\bf{I}}^{ - 1}}\left[ {\matrix{ {{L^b}} \cr {{M^b}} \cr {{N^b}} \cr } } \right] - {{\bf{I}}^{ - 1}}\left[ {\matrix{ 0 & { - R} & Q \cr R & 0 & { - P} \cr { - Q} & P & 0 \cr } } \right]{\bf{I}}\left[ {\matrix{ P \hfill \cr Q \hfill \cr R \hfill \cr } } \right] where L^b, M^b, N^b – total external moments with respect to the aircraft center of mass, I - inertia matrix of the form: 3I=[ Ixx−Ixy−Ixz−IyxIyy−Iyx−Izx−IzyIzz ]{\bf{I}} = \left[ {\matrix{ {{I_{xx}}} & { - {I_{xy}}} & { - {I_{xz}}} \cr { - {I_{yx}}} & {{I_{yy}}} & { - {I_{yx}}} \cr { - {I_{zx}}} & { - {I_{zy}}} & {{I_{zz}}} \cr } } \right]

It was assumed that the aircraft possesses the vertical plane of symmetry, so I_xy = I_zy = 0.

Quaternions were used to describe the altitude of the aircraft: 4[ e˙0e˙1e˙2e˙3 ]=−12[ 0PQR−P0−RQ−QR0−P−R−QP0 ][ e0e1e2e3 ]−kE[ e0e1e2e3 ]\left[ {\matrix{ {{{\dot e}_0}} \hfill \cr {{{\dot e}_1}} \hfill \cr {{{\dot e}_2}} \hfill \cr {{{\dot e}_3}} \hfill \cr } } \right] = - {1 \over 2}\left[ {\matrix{ 0 & P & Q & R \cr { - P} & 0 & { - R} & Q \cr { - Q} & R & 0 & { - P} \cr { - R} & { - Q} & P & 0 \cr } } \right]\left[ {\matrix{ {{e_0}} \hfill \cr {{e_1}} \hfill \cr {{e_2}} \hfill \cr {{e_3}} \hfill \cr } } \right] - kE\left[ {\matrix{ {{e_0}} \hfill \cr {{e_1}} \hfill \cr {{e_2}} \hfill \cr {{e_3}} \hfill \cr } } \right] where k – constant such as kh_int < 1 (h – the size of the numerical integration step). It was assumed that k = 1. In that way E=1−(e02+e12+e22+e32)E = 1 - \left( {e_0^2 + e_1^2 + e_2^2 + e_3^2} \right) is adjusted in each iteration, and it normalizes quaternion norm to 1.

Velocity components in O_nx_ny_nz_n are computed from the relation: 5[ x˙ny˙nzn ]=[ e02+e12−e22−e322(e1e2−e0e3)2(e0e2−e1e3)2(e0e3−e1e2)e02−e12+e22−e322(e2e3−e0e1)2(e1e3−e0e2)2(e0e1−e2e3)e02−e12−e22+e32 ][ UVW ]\left[ {\matrix{ {{{\dot x}_n}} \hfill \cr {{{\dot y}_n}} \hfill \cr {{z_n}} \hfill \cr } } \right] = \left[ {\matrix{ {e_0^2 + e_1^2 - e_2^2 - e_3^2} & {2\left( {{e_1}{e_2} - {e_0}{e_3}} \right)} & {2\left( {{e_0}{e_2} - {e_1}{e_3}} \right)} \cr {2\left( {{e_0}{e_3} - {e_1}{e_2}} \right)} & {e_0^2 - e_1^2 + e_2^2 - e_3^2} & {2\left( {{e_2}{e_3} - {e_0}{e_1}} \right)} \cr {2\left( {{e_1}{e_3} - {e_0}{e_2}} \right)} & {2\left( {{e_0}{e_1} - {e_2}{e_3}} \right)} & {e_0^2 - e_1^2 - e_2^2 + e_3^2} \cr } } \right]\;\left[ {\matrix{ U \cr V \cr W \cr } } \right]

2.3.3.

Aerodynamic Loads

The aerodynamic forces in the body-fixed frame are: 7F→ab=12ρVtot2STw→b[ −CXwlCYwl−CZwl ]\vec F_a^b = {1 \over 2}\rho V_{tot}^2S{{\bf{T}}_{{\bf{w}} \to {\bf{b}}}}\left[ {\matrix{ { - C_X^{wl}} \cr {C_Y^{wl}} \cr { - C_Z^{wl}} \cr } } \right] where ρ – air density, V_tot – airspeed, S_w – wing area, and CXwl,CYwl,CZwlC_X^{wl},C_Y^{wl},C_Z^{wl} – coefficients of the drag, side, and lift forces, respectively (in wind-laboratory frame O_wlx_wly_wlz_wl). The T_w→b is transformation matrix from O_wx_wy_wz_w frame to O_bx_by_bz_b frame: 8Tw→b=[ cosαcosβ−cosαsinβ−sinαsinβcosβ0sinαcosβ−sinαsinβcosα ]{{\bf{T}}_{{\bf{w}} \to {\bf{b}}}} = \left[ {\matrix{ {\cos \alpha \cos \beta } & { - \cos \alpha \sin \beta } & { - \sin \alpha } \cr {\sin \beta } & {\cos \beta } & 0 \cr {\sin \alpha \cos \beta } & { - \sin \alpha \sin \beta } & {\cos \alpha } \cr } } \right] where α, β – angles of attack and sideslip, respectively. Aerodynamic moments with respect to the center of mass are: 9M→a,Ob=12ρVtot2STw→b[ bwCLwcwCMwbwCNw ]\vec M_{a,O}^{\rm{b}} = {1 \over 2}\rho V_{tot}^2S{{\bf{T}}_{{\bf{w}} \to {\bf{b}}}}\left[ {\matrix{ {{b_w}C_L^w} \cr {{c_w}C_M^w} \cr {{b_w}C_N^w} \cr } } \right] where CLw,CMw,CNwC_L^w,C_M^w,C_N^w – rolling, pitching, and yawing moment coefficients, respectively (in wind frame). The aerodynamic coefficients of forces were calculated as [42,43]: 10CXwl=CX0wl+CXαwlα+CXα2wlα2+CXQwlQcw2Vtot+CXδEwlδEC_X^{wl} = C_{{X_0}}^{wl} + C_{{X_\alpha }}^{wl}\alpha + C_{{X_{\alpha 2}}}^{wl}{\alpha ^2} + C_{{X_Q}}^{wl}{{Q{c_w}} \over {2{V_{tot}}}} + C_{{X_{{\delta _E}}}}^{wl}{\delta _E} 11CYwl=CY0wl+CYβwlβ+CYPwlPbw2Vtot +CYRwlRbw2Vtot +CYδAwlδA+CYδRwlδRC_Y^{wl} = C_{{Y_0}}^{wl} + C_{{Y_\beta }}^{wl}\beta + C_{{Y_P}}^{wl}{{P{b_w}} \over {2{V_{{\rm{tot }}}}}} + C_{{Y_R}}^{wl}{{R{b_w}} \over {2{V_{{\rm{tot }}}}}} + C_{{Y_{{\delta _A}}}}^{wl}{\delta _A} + C_{{Y_{{\delta _R}}}}^{wl}{\delta _R} 12CZwl=CZ0wl+CZαwlα+CZQwlQcw2Vtot+CZδEwlδEC_Z^{wl} = C_{{Z_0}}^{wl} + C_{{Z_\alpha }}^{wl}\alpha + C_{{Z_Q}}^{wl}{{Q{c_w}} \over {2{V_{tot}}}} + C_{{Z_{{\delta _E}}}}^{wl}{\delta _E} where δ_E – elevator deflection (positive trailing edge down), δ_A – aileron deflection (positive trailing edge down) calculated as mean of deflections for left and right ailerons: 13δA=δAL−δAR2{\delta _A} = {{{\delta _{AL}} - {\delta _{AR}}} \over 2}

Additionally, δ_R – rudder deflection (positive trailing edge left).

The moment’s coefficients are: 14Clw=Cl0w+Clβwβ+ClPwPbw2Vtot +ClRwRbw2Vtot +ClδAwδA+ClδRwδRC_l^w = C_{{l_0}}^w + C_{{l_\beta }}^w\beta + C_{{l_P}}^w{{P{b_w}} \over {2{V_{{\rm{tot }}}}}} + C_{{l_R}}^w{{R{b_w}} \over {2{V_{{\rm{tot }}}}}} + C_{{l_{{\delta _A}}}}^w{\delta _A} + C_{{l_{{\delta _R}}}}^w{\delta _R} 15Cmw=Cm0w+Cmαwα+CmQwQcw2Vtot+CmδEwδEC_m^w = C_{{m_0}}^w + C_{{m_\alpha }}^w\alpha + C_{{m_Q}}^w{{Q{c_w}} \over {2{V_{tot}}}} + C_{{m_{{\delta _E}}}}^w{\delta _E} 16Cnw=Cn0w+Cnβwβ+CnPwPbw2Vtot +CnRwRbw2Vtot +CnδAwδA+CnδRwδRC_n^w = C_{{n_0}}^w + C_{{n_\beta }}^w\beta + C_{{n_P}}^w{{P{b_w}} \over {2{V_{{\rm{tot }}}}}} + C_{{n_R}}^w{{R{b_w}} \over {2{V_{{\rm{tot }}}}}} + C_{{n_{{\delta _A}}}}^w{\delta _A} + C_{{n_{{\delta _R}}}}^w{\delta _R}

The total speed is: 17Vtot=(U−UW)2+(V−VW)2+(W−WW)2{V_{tot}} = \sqrt {{{\left( {U - {U_W}} \right)}^2} + {{\left( {V - {V_W}} \right)}^2} + {{\left( {W - {W_W}} \right)}^2}} where U_w, V_w, W_w – wind velocity vector components in the body-fixed frame O_bx_by_bz_b. A detailed description of the wind model can be found in [44,45]. The angle of attack was calculated as: 18α=atanW−WWU−UW\alpha = {\mathop{\rm atan}\nolimits} {{W - {W_W}} \over {U - {U_W}}} and the angle of the sideslip is: 19β=asinV−VWVtot\beta = {\mathop{\rm asin}\nolimits} {{V - {V_W}} \over {{V_{tot}}}}

The aerodynamic characteristics of the plane were calculated using XFLR5 software. The coefficients (for undeflected control surfaces) are presented in Fig. 10.

Figure 10.

Aerodynamic coefficients of the airplane (a) drag (b) side force (c) lift force (d) rolling moment (e) pitching moment (f) yawing moment

Force coefficients are related in O_wlx_wly_wlz_wl frame, and moments (with respect to the center of mass) in O_wlx_wly_wlz_wl frame.

The stability and control derivatives of the Funcub plane are presented in Table 1.

Table 1.

Longitudinal	Value	Unit	Lateral–directional	Value	Unit
CX0WlC_{{X_0}}^{Wl}	0,01920	–	CY0wlC_{{Y_0}}^{wl}	0	–
CXαwlC_{{X_\alpha }}^{wl}	0,01653	1/rad	CYBwlC_{{Y_B}}^{wl}	–0,26240	1/rad
CXα2wlC_{{X_{{\alpha _2}}}}^{wl}	1,53900	1/rad2	CYPwlC_{{Y_P}}^{wl}	–0,07607	1/rad
CXQwlC_{{X_o}}^{wl}	0,29927	1/rad	CYRwlC_{{Y_R}}^{wl}	0,22320	1/rad
CXδEwlC_{{X_{{\delta _E}}}}^{wl}	0,00989	1/rad	CYδAwlC_{{Y_{{\delta _A}}}}^{wl}	0	1/rad
CZ0wlC_{{Z_0}}^{wl}	0,14840	–	CYδRwlC_{{Y_{{\delta _R}}}}^{wl}	0,11450	1/rad
CZαWLC_{{Z_\alpha }}^{WL}	4,22000	1/rad	Cl0WC_{{l_0}}^W	0	–
CZQwlC_{{Z_Q}}^{wl}	7,08930	1/rad	ClβwC_{{l_\beta }}^w	–0,00789	1/rad
CZδEwlC_{{Z_{{\delta _F}}}}^{wl}	0,49847	1/rad	ClPwC_{{l_P}}^w	–0,46079	1/rad
Cm0wC_{{m_0}}^w	0,03794	–	ClRwC_{{l_R}}^w	0,10362	1/rad
CmαwC_{{m_\alpha }}^w	–1,68900	1/rad	ClδAwC_{{l_{{\delta _A}}}}^w	–0,14324	1/rad
CmQwC_{{m_Q}}^w	–9,69450	1/rad	ClδRwC_{{l_{{\delta _R}}}}^w	0,01302	1/rad
CmδEwC_{{m_{{\delta _F}}}}^w	–1,20320	1/rad	Cn0wC_{{n_0}}^w	0	–
			CnβwC_{{n_B}}^w	0,14330	1/rad
			CnPwC_{{n_P}}^w	–0,02406	1/rad
			CnRwC_{{n_R}}^w	–0,10615	1/rad
			CNδAwC_{{N_{{\delta _A}}}}^w	0	1/rad
			CnδRwC_{{n_{{\delta _R}}}}^w	–0,05460	1/rad

2.3.4.

Propulsion Model

Data of the propeller were obtained from the dataset provided by the manufacturer (APC Propellers). The propulsive forces were calculated as: 20F→pb=[ Tp,0,0 ]\vec F_p^b = \left[ {\matrix{ {{T_p},} \hfill & {0,} \hfill & 0 \hfill \cr } } \right]

The propeller thrust was calculated using the formula: 21Tp=ρn2Dp4CT{T_p} = \rho {n^2}D_p^4{C_T} where ρ – air density, n – angular rate of the propeller, D_p – propeller diameter, and C_T – thrust coefficient (Fig. 11).

Figure 11.

Thrust coefficient

The moments with respect to the aircraft center of mass were also taken into account: 22M→p,Ob=r→pb×F→pb+Jpω→p×Ω→\vec M_{p,O}^b = \vec r_p^b \times \vec F_p^b + {J_p}{{\vec \omega }_p} \times \vec \Omega where r→pb=[ 0,370,0 ]T\vec r_p^b = {\left[ {\matrix{ {0,37} \hfill & {0,} \hfill & 0 \hfill \cr } } \right]^{\rm{T}}} – position vector of the propeller hub in body-fixed frame, J_p – moment of inertia of the rotating masses, ω→p=[ ω,0,0 ]{{\vec \omega }_p} = \left[ {\matrix{ {\omega ,} \hfill & {0,} \hfill & 0 \hfill \cr } } \right] – the angular rate of the propeller and Ω→=[P,Q,R]\vec \Omega = [P,\quad Q,\quad R] – vector of angular rates of the plane.

The presented mathematical model of the aircraft was implemented in MATLAB/SIMULINK R2023a with Update 5. The equations of motion were integrated numerically using a 4th-order fixed-step solver (Runge-Kutta) with a time step size of 0,001 s.

2.3.5.

Solar Radiation Model

The solar radiation model was developed to predict the amount of received solar radiation. The position of the Sun was parametrized using azimuth and elevation angles (Fig. 12) [23].

Figure 12.

Solar azimuth and elevation angles

Solar declination was calculated as follows [46,47]: 23δ=23,45sin(360365(284+n))\delta = 23,45\sin \left( {{{360} \over {365}}(284 + {\rm{n}})} \right) where n = 1… 365 is day of the year (1^st January = 1). Next, the Sun hour angle t_ST was estimated using the relation: 24tST=tLST+tET+(lST−lLocal )4min/°{t_{ST}} = {t_{LST}} + {t_{ET}} + \left( {{l_{ST}} - {l_{Local{\rm{ }}}}} \right)4min{/^^\circ } where: t_LST – Local Standard Time, t_ET – equation of time, l_ST – Standard Longitude and l_Local – Local Longitude. The equation of time is given in the form [47]: 25tET=229,2(0,000075+0,001868cosB−0,032077sinB−0,014615cos2B−0,04089sin2B)\matrix{ {{t_{ET}}} \hfill & { = 229,2(0,000075 + 0,001868\cos B} \hfill \cr {} \hfill & { - 0,032077\sin B - 0,014615\cos 2B} \hfill \cr {} \hfill & { - 0,04089\sin 2B)} \hfill \cr } where B=(n−1)360365B = (n - 1){{360} \over {365}}. Solar elevation angle α_el is [46, 48]: 26αel=asin(sinδsinφ+cosδcosφcosh){\alpha _{el}} = {\mathop{\rm asin}\nolimits} (\sin \delta \sin \varphi + \cos \delta \cos \varphi \cos h)

Azimuth angle a_az was obtained from the relation [48]: 27aaz=acos(sinδcosφ−cosδsinφcosh)cosαel{a_{az}} = {\mathop{\rm acos}\nolimits} {{(\sin \delta \cos \varphi - \cos \delta \sin \varphi \cos h)} \over {\cos {\alpha _{el}}}}

The inverse cosine trigonometric function returns values in the range from 0° to 180°. If the hour angle is negative, then a_az = 360°–a_az. In that way a_az falls in range from 0° to 360°.

The ASHRAE clear sky model was used to predict the amount of solar radiation received by the sensors. The changes in extraterrestrial radiation during the year are modeled as [46]: 28I=I0(1+0,034cos(360n365,25))I = {I_0}\left( {1 + 0,034\cos \left( {{{360n} \over {365,25}}} \right)} \right) where the Solar constant is I₀ = 1353 W/m². Optical air mass m was calculated in the following way [46]: 29m=1sinαel+0,50572(6,07995+αel)−1,6364m = {1 \over {\sin {\alpha _{el}} + 0,50572{{\left( {6,07995 + {\alpha _{el}}} \right)}^{ - 1,6364}}}} where α_el – solar elevation angle (this value should be expressed in degrees).

The total irradiance is the sum of the beam, diffused, and reflected components [49,50]. The beam’s normal irradiance per unit area, (which is perpendicular to the sun rays) [46,51] was obtained as: 30Ib,N=Ie−τbmb{I_{b,N}} = I{e^{ - {\tau _b}{m^b}}}

The diffuse horizontal irradiance per unit area on a horizontal surface [46,51]: 31Id,h=Ie−τdmd{I_{d,h}} = I{e^{ - {\tau _d}{m^d}}}

The coefficients τ_b and τ_d depend on the geographical location and vary during the year. Some tabulated values can be found, for example, in [46]. The coefficients b and d were obtained from the following relations [37,46,51]: 32b=1,219−0,043τb−0,151τd−0,204τbτdb = 1,219 - 0,043{\tau _b} - 0,151{\tau _d} - 0,204{\tau _b}{\tau _d} 33d=0,202+0,852τb−0,007τd−0,357τbτdd = 0,202 + 0,852{\tau _b} - 0,007{\tau _d} - 0,357{\tau _b}{\tau _d}

In O_bx_by_bz_b frame the unit vector along the z_b axis (in the negative direction) have the following components n^=[ 0,0,−1 ]T\hat n = {\left[ {\matrix{ {0,} \hfill & {0,} \hfill & { - 1} \hfill \cr } } \right]^{\rm{T}}} [28]. Next, the value of the cosine of the angle between the direction of the incoming solar rays and the line that is perpendicular to the surface of the wings was calculated as [29,37,39,51,52]: 34cosλ=n^·[ 1000cosΦ−sinΦ0sinΦcosΦ ][ cosΘ0sinΘ010−sinΘ0cosΘ ][ cosΨ−sinΨ0sinΨcosΨ0001 ][ cosαazcosαelsinαazcosαelsinαel ]\matrix{ {\cos \lambda = } \hfill & {\hat n\cdot\left[ {\matrix{ 1 & 0 & 0 \cr 0 & {\cos \Phi } & { - \sin \Phi } \cr 0 & {\sin \Phi } & {\cos \Phi } \cr } } \right]} \hfill \cr {} \hfill & {\left[ {\matrix{ {\cos \Theta } & 0 & {\sin \Theta } \cr 0 & 1 & 0 \cr { - \sin \Theta } & 0 & {\cos \Theta } \cr } } \right]} \hfill \cr {} \hfill & {\left[ {\matrix{ {\cos \Psi } & { - \sin \Psi } & 0 \cr {\sin \Psi } & {\cos \Psi } & 0 \cr 0 & 0 & 1 \cr } } \right]} \hfill \cr {} \hfill & {\left[ {\matrix{ {\cos {\alpha _{az}}\cos {\alpha _{el}}} \cr {\sin {\alpha _{az}}\cos {\alpha _{el}}} \cr {\sin {\alpha _{el}}} \cr } } \right]} \hfill \cr }

After several mathematical operations, the result is [51,52]: 35cosλ=cosαesinαs(cosΨsinΦ−cosΦsinΨsinΘ)−cosαecosαs(sinΦsinΨ+cosΦcosΨsinΘ)+cosΦsinαecosΘ\matrix{ {\cos \lambda = } \hfill & {\cos {\alpha _e}\sin {\alpha _s}(\cos \Psi \sin \Phi } \hfill \cr {} \hfill & { - \cos \Phi \sin \Psi \sin \Theta )} \hfill \cr {} \hfill & { - \cos {\alpha _e}\cos {\alpha _s}(\sin \Phi \sin \Psi } \hfill \cr {} \hfill & { + \cos \Phi \cos \Psi \sin \Theta )} \hfill \cr {} \hfill & { + \cos \Phi \sin {\alpha _e}\cos \Theta } \hfill \cr }

When the cos λ is known, it can be used to calculate the amount of received power from the unit area [32, 51,53].

The sensors are installed only on the upper surface of the wing. If the surface of the sensor is oriented in the direction opposite to the Sun, then the direct component is [26,37,54]: 36Pb={ Ib,Ncosλcosλ≥00cosλ<0{P_b} = \left\{ {\matrix{ {{I_{b,N}}\cos \lambda } \hfill & {\cos \lambda \ge 0} \hfill \cr 0 \hfill & {\cos \lambda < 0} \hfill \cr } } \right.

The diffuse irradiance is: 37Pd=Id,hcos2λ{P_d} = {I_{d,h}}{\cos ^2}\lambda

The solar radiation on the horizontal surface is: 38Ih=Ib,Nsinαe+Id,h{I_h} = {I_{b,N}}\sin {\alpha _e} + {I_{d,h}}

For low attitude also, the reflected component was taken into account [55]: 39Pr=ρrIhsin2λ{P_r} = {\rho _r}{I_h}{\sin ^2}\lambda where ρ_r – ground reflection factor (assumed 0,05). The total received irradiance is: 40Psa=Pb+Pd+Pr{P_{sa}}\;{\rm{ = }}\;{P_b}{\rm{ + }}{P_d}{\rm{ + }}{P_r}

The solar radiation model was implemented in MATLAB/SIMULINK and integrated with the numerical flight simulation of the aircraft (Fig. 13).

Figure 13.

Solar radiation model implemented in SIMULINK

Moreover, the presented model of the Sun’s position was validated by comparing the results obtained from SIMULINK with those from the website https://www.suncalc.org. Next, the model was used to predict the amount of solar radiation and prepare the flight scenarios for the real experiments.

During research, it was required to estimate the power demand and energy consumption of the individual subsystems (electric engine, autopilot, telemetry, and sensors). The total energy consumption of the system is recorded in the flight logs. However, the power module used to monitor the battery state must be properly calibrated.

A series of ground tests was conducted to obtain the amount of power required by the components of the electric systems. Voltage and current were measured manually using a voltmeter (connected in parallel) and an ammeter (placed in series) between the main battery and autopilot.

It should be noted that a secondary power source powers the solar radiation measurement system. The obtained results are presented in Table 2.

The energy consumption by the electronic components is practically constant. The power required by electronics is much smaller (below 3% of the total maximum value) than the power necessary to drive the electric motor. Increasing the angular speed of the motor results in a significant drop in the battery voltage.

A series of flight tests was conducted at different times of the day to collect the data. The experiments were evaluated at two locations: Zalesie (52,266874°, 20,751108°) and Przasnysz (53,009749°, 20,929911°) airfields (Fig. 14).

Figure 14.

Aircraft ready to start at the runway in Przasnysz airfield.

The primary objectives of the experiments were to collect data on flight parameters and solar radiation, which are essential for validating the developed mathematical models.

The weather conditions were good, with low cloud coverage (in both locations). In Zalesie, the wind gusts’ maximum speed was no more than 3 m/s, but in Przasnysz, it was significantly higher (up to 6 m/s). The aircraft was controlled remotely by an experienced drone pilot. The experiments were realized in two flight modes: manual and semi-automatic. Additionally, the second operator observed the flight parameters (e.g., speed, altitude, battery state of charge) using the Ground Control Station. In that way, the aircraft ensured that it would correctly and safely complete the predefined mission scenario. The flight logs obtained from the autopilot were analyzed using Mission Planner software (https://ardupilot.org/planner/) and UAV Log Viewer (https://plot.ardupilot.org/#/).

Next, the developed models were validated using the obtained data. The input data to the model (geographic coordinates, date, and time of the trials) were set according to the measured data. The flight trials were documented with photographs, so it was possible to obtain the level of cloudiness. The flight logs were aligned in time with data from the solar radiation measurement system. In order to calculate the cross-correlation of two discrete-time signals (aileron deflection), the MATLAB function “xcorr” was used.

In the subsequent paragraphs, selected results from the experiments are presented.

This test was performed in Zalesie on 2023-09-08 (year-month-day). The data logging process started at 16:38:56 (HH:MM:SS, hours, minutes, and seconds) Central European Summer Time (CEST).

The aircraft took off at 16:45:54. In that case, the UAV was piloted only in manual mode. The aircraft flight path is presented in Fig.15a (data has been cropped to the time range of interest). As can be observed, the pilot has frequently repeated trajectories similar to airfield traffic patterns (in the counterclockwise direction when looking from the top). Small yellow arrows in Fig.15a mean the direction of motion. The white arrow indicates the solar azimuth (the direction of the sun’s rays). At the start, the solar azimuth was 249,57° and the elevation was 22,04°. Azimuth is measured from North in a clockwise direction, and elevation from the horizontal plane (please see the convention in Fig.12). The roll, pitch, and yaw angles are presented in Fig.15b. In Fig.15c, the battery voltage, current, and the total amount of energy consumed by the aircraft are shown.

Figure 15.

Aircraft flight parameters (a) position (b) roll, pitch, yaw angles (c) battery voltage, current, and consumed energy (case 1)

The measured roll angle was between -48,16° and +20,74° (the plot is asymmetrical because the aircraft performed left turns). From the yaw angle plot, it can be concluded that approximately 8 full loops were evaluated in flight.

The voltage decreases slowly with time. The current has the highest amplitude at the beginning of the mission because then the aircraft performs the start and the ascending flight. To realize the mission 7.84 Wh of energy was spent.

In Fig.16, the measured solar illuminance is presented. The measured solar illuminance was converted into solar radiation intensity using calibration coefficients that were obtained experimentally. Then, the obtained results were integrated in time to estimate the energy that could be collected (per 1 square meter).

Figure 16.

Solar illuminance, radiation intensity, and energy (case 1)

The maximum possible amount of solar radiation intensity is approximately 554,6 W/m². The amplitude varies significantly due to changes in altitude. The radiation received by sensors strongly correlates with the roll and pitch angles. The aircraft performed a set of maneuvers. The intensity depends strongly on the angle between the incoming ray and the sensor’s surface.

Moreover, it can be observed that the data from both sensors differ a bit. The main reason for that phenomenon is the wing’s dihedral angle and the sensors’ minor mounting imperfections.

The amount of energy increases monotonically with time. The data from sensor 1 indicate that collecting 14.78 Wh of energy (from 1 square meter) is possible in these conditions. Using data from sensor 2, it was found that this result is slightly higher, 15.41 Wh.

The experimental radiation intensity values are in close agreement with the model predictions. The model overestimated the amount of collected energy (17,35 Wh).

The typical overall efficiency of photovoltaic systems that are used on solar-powered UAVs is around 20% [7]. This means that theoretically, the 2.956 Wh of energy can be obtained (using sensor 1 data).

The Multiplex Funcub NG is unable to achieve perpetual flight using only solar energy because the energy required for propulsion is significantly higher than the amount of energy that can be harvested from solar radiation.

It is worthwhile to take a closer look at the results to understand the relationship between the Euler angles and the amount of solar radiation received. The selected data portion (from 20 s to 50 s) is shown in Fig. 17 and Fig. 18 to make the details more visible.

Figure 17.

Roll, pitch, yaw angles (case 1, enlarged view)

Figure 18.

Solar illuminance, radiation intensity, and energy (case 1, enlarged view)

One of the peaks of the solar radiation occurred approximately at time 31.6 s (please see the datatip in Fig. 17). Then, the aircraft was banked on the leftwing (roll angle is negative and equal to around –32°) and moved at the yaw angle of 350° (from South to North). That means the solar radiation sensors were pointed in the direction of the oncoming radiation. The sketch of the aircraft location on the analyzed portion of the flight path is presented in Fig. 19.

Figure 19.

Approximate aircraft location at energy peak (case 1, enlarged view)

On the other hand, in 50 s, the amount of received solar radiation is relatively small because the aircraft moved from North to South, and the roll angle is ≈4.5° (the sensors were not pointed directly to the Sun).

Scenario 2 was also completed in Zalesie on 2023-09-12. The data logging in autopilot was initiated at 13:07:02 CEST. The plane took off at 13:14:42. In this experiment, the aircraft was piloted in both manual and semi-automatic modes. The resulting flight path is presented in Fig. 20a. The plane moved in a counterclockwise direction. The solar position (azimuth and elevation) was 193,19° and 41,23°, respectively.

Figure 20.

Aircraft flight parameters (a) position (b) roll, pitch, yaw angles (c) battery voltage, current, and consumed energy (case 2)

The roll angle (Fig. 20b) was in the range from – 56,25° to +29,17°. Also, the pitch angle was between –41.31° to +31.09°. The yaw angle indicates that the aircraft realized more than 11 loops.

The current (Fig. 20c) drops rapidly after 320 s because the plane landed and then was taxiing on the runway. The aircraft spent 10,83 Wh of energy during the mission.

In Fig. 21, the solar radiation measured in the second flight trial is shown.

Figure 21.

Solar illuminance, radiation intensity, and energy (case 2)

The solar radiation intensity exhibits a sinusoidal pattern corresponding to changes in the orientation angles. The maximum value of solar radiation intensity is 829,3 W/m² (significantly higher than in case 1). From sensor 1, it was estimated that collecting 45,21 Wh of solar energy is possible (per 1 square meter). Using sensor 2, this value was estimated to be 45,26 Wh.

The model predicted that the amount of energy that can be harvested is 45,65 Wh, which is close to the real values.

Again, assuming the aircraft can be equipped with solar panels (1 square meter in area) and the energy can be harvested with an efficiency of 20%, the obtained value is 9,042 Wh.

In this case, the amount of harvested energy is still insufficient to cover the instantaneous energy demand. Moreover, from a practical point of view, the hypothetical system cannot accumulate energy for the onboard batteries (because the energy consumption of 10,83 Wh was higher than production). As a matter of fact, it would not be possible to achieve perpetual flight functionality over a longer period.

In Fig. 22 and Fig. 23, the selected data portion (from 220 s to 260 s) is shown in an enlarged view.

Figure 22.

Roll, pitch, yaw angles (case 2, enlarged view)

Figure 23.

Solar illuminance, radiation intensity, and energy (case 2, enlarged view)

One of the energy peaks occurred at time 245,6 s. The aircraft performed a left turn and was banked on the leftwing (roll angle is –27,25°, please see Fig. 22). This means the solar radiation sensors were pointed toward the incoming sunlight. The situation scheme with the selected portion of the flight path is presented in Fig. 24.

Figure 24.

Approximate aircraft location at energy peak (case 2, enlarged view)

In case number 3, the experiment took place in Przasnysz on 2023-09-21. The data logging system was initiated at 15:41:26, and the plane taxied for a while. Next, the aircraft took off at 15:45:17. The plane started and landed in manual mode, but most of the mission was realized in loiter mode. The radius of the circular trajectory (Fig. 25a) was larger when compared to case 2. Solar azimuth and elevation angles were 234,53° and 24,41°.

Figure 25.

Aircraft flight parameters (a) position (b) roll, pitch, yaw angles (c) battery voltage, current, and consumed energy (case 3)

The roll angle (Fig. 25b) varied significantly from –45,75° up to 73,12°. To fulfill the mission, 12,28 Wh of the energy was spent by the onboard electronic system (Fig. 25c).

The light radiation intensity is presented in Fig. 26.

Figure 26.

Solar illuminance, radiation intensity, and energy (case 3)

The maximum measured value of the radiation intensity is 701,6 W/m².

Using data from sensor 1, it was calculated that 20,93 Wh of energy can be collected (from 1 square meter of surface). Data from sensor 2 indicates a bit lower value of 19,01 Wh. The model predictions (20,72 Wh) are between the values from both sensors.

The enlarged view of the selected data portion is shown in Fig. 27 and Fig. 28.

Figure 27.

Roll, pitch, yaw angles (case 3, enlarged view)

Figure 28.

Solar illuminance, radiation intensity, and energy (case 3, enlarged view)

For the selected part of the trajectory, the aircraft performed two rapid turns (please see Fig. 29). In the first turn, it was banked to the left, and as a result, the wing shaded the solar radiation sensors. Between 40 s and 50 s, the irradiance was smaller than 60 W/m² (the curve is practically flat for this time range). In the second turn, the aircraft banked on the right wing. The maximum roll angle occurred at 55.7 s and is equal to 36,06°.

Figure 29.

Approximate aircraft location at energy peak (case 3, enlarged view)

As a result, the sensors were pointed directly at the Sun. In Fig. 28, at 55,7 s, the local maximum of the signal was observed.