**Modeling and Investigations on Surface Colors of Wings on the Performance of Albatross-Inspired Mars Drones and Thermoelectric Generation Capabilities**

#### **Devyn Rice 1, Samah Ben Ayed 2, Stephen Johnstone <sup>3</sup> and Abdessattar Abdelkefi 1,\***


Received: 28 July 2020; Accepted: 13 August 2020; Published: 16 August 2020

**Abstract:** Thermal effects of wing color for Albatross-inspired drones performing in the Martian atmosphere are investigated during the summer and winter seasons. This study focuses on two useful consequences of the thermal effects of wing color: the drag reduction and the thermoelectric generation of power. According to its color, each wing side has a certain temperature affecting the drag. Investigations of various configurations have shown that the thermal effect on the wing boundary layer skin drag is insignificant because of the low atmospheric pressure. However, the total drag varies as much as 12.8% between the highest performing wing color configuration and the lowest performing configuration. Additionally, the large temperature differences between the top and the bottom wing surfaces show great potential for thermoelectric power generation. The maximum temperature differences between the top and bottom surfaces for the summer and winter seasons are, respectively, 65 K and 30 K. The drag reduction and the power generation via thermoelectric generators both contribute to enhancing the endurance of drones. Future drone designs will benefit from increased endurance through optimizing the wing color configuration.

**Keywords:** drones; Martian atmosphere; thermal analysis; aerodynamics; thermoelectric generators

#### **1. Introduction**

Over the past few decades, various government and private programs have been developed for space exploration. Among the high interest exploration targets is the planet Mars. NASA has performed many missions in the past and will continue with new mission in the future including the upcoming Mars 2020 mission [1]. The private company, SpaceX, hopes to establish a supply depot on Mars by the year 2024 as a cornerstone for future manned missions to the planet [2]. As Mars begins to generate greater interest in the near future, planetary exploration will become much more important. The current methods used to explore the Martian surface include rover, lander, and satellite [3]. In the past few years, the concept of using drones for planetary exploration has gained traction and many design concepts have been developed [3]. To design a drone suitable for Mars exploration, it is important to understand the Martian atmosphere and other properties for flight. Compared to Earth, the atmosphere of Mars is much thinner, and the temperature is also colder [4]. Atmospheric properties such as density, pressure, and temperature were obtained at varying altitudes from the high-speed entry of a probe into the atmosphere [5].

Due to the extreme conditions and challenges associated with flying in the Martian atmosphere, maximizing the efficiency of a drone is essential [3,4]. Drones require high endurance to be able to

perform many defined missions for planetary exploration. The endurance of a drone is primarily limited by the available energy [6]. There are various propulsion systems and fuel/oxidizers that provide energy for drones, and there are additional methods, such as solar cells, batteries, etc., that provide energy to the drone [7,8]. Solar cells have been gaining popularity for use in space exploration for several systems. Orbitals and drones have had studies performed for Solar Electric Propulsion (SEP) applications [9,10]. When considering enhancing the endurance of a drone, there are various methods to take into consideration. Optimizing the propulsion system to increase the total available energy is one option. In the case of an electric propulsion system, recharging the battery by various means can enhance the endurance of the drone. Solar panel methods have been researched as a possible option to power or partially power drones [11–14]. Additionally, thermoelectric generators are another possible method to harness energy via converting the thermal energy from the temperature difference between two surfaces to electrical power [15]. Another option for enhancing the endurance of a drone is to optimize the structure so that energy can be best conserved which will result in a reduction of necessary power [16–21]. Control algorithms have also a great impact on the energy consumption [22,23].

Researchers have explored various techniques to decrease and manage the required power for drones. An example is to propose new design methods for drones, which improve their aerodynamic and mechanical efficiencies and consequently reduce the power consumption. Therefore, design methods enhance the performance and efficiency of drones which will extend their endurance and improve their utility in complex environments. Reduction of drag can be considered one of the main factors during the design process to conserve energy for the mission. In other words, drag reduction of fixed wing drones is crucial to flight efficiency on Mars [24].

The structure of a fixed wing drone can come in many different shapes and sizes. Often, when trying to optimize a drone, inspiration can be sought from nature. As seen from nature, on Earth there exists different modes of flight depending on the species of bird including flapping, hovering, gliding, and soaring [25]. Specifically, soaring and gliding are flight modes that are used by fixed wing drones. The nature of a bird's mode of flight was used in a study to design a drone capable of wing morphing. Two different flight conditions were achieved through either sweeping or spanning the wing [26]. Another study was previously performed to analyze the thermal effects of wing color on the resulting heated boundary layer flow over the wing of an Albatross [27]. The study analyzed different wing color configurations to identify the most efficient configuration. The results verify the drag efficiency of the natural wing color configuration of the Albatross (which is black on the top surface and white on the bottom surface) compared to other color configurations. To design an efficient drone capable of flying in the Martian atmosphere, a similar study will be performed to analyze the thermal effects of wing color on the heated boundary layer of an Albatross wing shape within the Martian atmosphere. To the authors' knowledge, a thermal boundary layer analysis of wing color configurations for drag efficiency calculations of a drone based on wing color on Mars has not been previously performed. This study will provide information on the surface temperature of a wing during different seasons and for different color configurations. An optimal wing color configuration will be determined to achieve maximum drag efficiency. Additionally, due to a possible temperature difference between the top and bottom wing surfaces, an investigation into the application of thermoelectric generators (TEGs) will be performed. Thermoelectric generators are composed of semiconductor material that produces thermoelectric modules. These modules are between ceramic wafers and then connected to an electrical load resistance. Electrical power is generated from the temperature difference as a result that can be explained through the Seebeck effect. The specific type of TEG considered for this study is the BiTe which is considered suitable for low-temperature heat recovery applications [28–30].

The paper is organized as follows: the thermal model describing the energy balance, wing color configurations, and the surrounding temperatures and heat fluxes is discussed in Section 2. Section 3 considers both the wing surfaces to investigate the skin drag forces and the effects of the angle of attack. A 3D aerodynamic analysis is performed to determine the total drag forces in Section 4. Section 5

studies the application of thermoelectric generators for drones in the Martian atmosphere. A summary and conclusions are drawn in Section 6.

#### **2. Mathematical Model**

An analytical model is developed using an energy balance to determine the wing surface temperature and consequently the skin friction drag. Atmospheric properties, such as sky, ground, and ambient temperatures are determined and applied to the energy balance. Each section of the energy balance is discussed in detail for the following sections: solar energy absorption, convective heat transfer, and radiative heat transfer. The conductive heat transfer will be defined in Section 4.

#### *2.1. Energy Balance*

To understand the thermal effects of different wing colors, a thermal analysis is performed using an energy balance and including the influences of solar energy absorption, convection heat loss, and radiative heat loss. The energy balance on a drone fixed wing can be expressed as:

$$\mathbf{P}\_{\text{solar}} = \mathbf{P}\_{\text{convective}} + \mathbf{P}\_{\text{radiative}} + \mathbf{P}\_{\text{conductive}} \tag{1}$$

The various sources of energy acting on the wing can be seen in Figure 1. The absorbed solar radiation is balanced by the convective, radiative and conductive transfer. Determining the various fluxes in the energy balance requires information from the Martian atmosphere. In this study, the ambient temperature, pressure, wind speed, and other atmospheric properties are obtained from the dataset collected by Viking Lander 1 and the entry probe (VL1), which is located at 22.37 N latitude and 47.97 W longitude on the surface of Mars [31].

**Figure 1.** Heat flux on a wing.

The top and the bottom wing surfaces experience different heat fluxes, which can be determined by the view factors, surface temperatures, surface absorptance, and surface emissivity. The total energy balance equations for both the top and the bottom surfaces of the wing are given in Equations (2) and (3), respectively [32].

$$\begin{aligned} \alpha\_{\text{t}} \cdot \mathbf{G}\_{\text{b}} \cdot \cos \theta &+ \alpha\_{\text{t}} \cdot \mathbf{F}\_{\text{d}}^{\text{t}} \cdot \mathbf{G}\_{\text{dh}} + \alpha\_{\text{t}} \cdot \text{al} \cdot \mathbf{F}\_{\text{al}}^{\text{t}} \cdot (\mathbf{G}\_{\text{b}h} + \mathbf{G}\_{\text{dh}}) - \mathbf{h} \cdot (\mathbf{T}\_{\text{s} \rightarrow \text{top}} - \mathbf{T}\_{\text{amb}}) - \varepsilon\_{\text{t}} \cdot \sigma \\\ \mathbf{F}\_{\text{d}}^{\text{t}} \cdot (\mathbf{T}\_{\text{s} \leftarrow \text{top}}^{4} - \mathbf{T}\_{\text{sky}}^{4}) - \varepsilon\_{\text{t}} \cdot \sigma \cdot \mathbf{F}\_{\text{al}}^{\text{t}} \cdot (\mathbf{T}\_{\text{s} \leftarrow \text{top}}^{4} - \mathbf{T}\_{\text{g} \leftarrow \text{fl}}^{4}) - \frac{1}{\mathcal{R}\_{\text{eff}}} \cdot (\mathbf{T}\_{\text{s} \leftarrow \text{top}} - \mathbf{T}\_{\text{s} \leftarrow \text{bottom}}) = 0 \end{aligned} \tag{2}$$

$$\begin{array}{c} \text{c}\_{\text{b}}\text{-}\text{F}\_{\text{d}}^{\text{b}}\text{-}\text{G}\_{\text{dh}} \quad + \text{c}\_{\text{b}}\text{-}\text{al}\text{-}\text{F}\_{\text{al}^{\text{l}}}^{\text{b}}(\text{C}\_{\text{bh}} + \text{G}\_{\text{dh}}) - \text{h}\cdot(\text{T}\_{\text{s-bottom}} - \text{T}\_{\text{amb}})\\ \qquad \qquad \qquad \qquad \qquad - \text{c}\_{\text{b}}\cdot\sigma\cdot\text{F}\_{\text{d}}^{\text{b}}\cdot(\text{T}\_{\text{s-bottom}}^{4} - \text{T}\_{\text{sky}}^{4})\\ \qquad \qquad \qquad \qquad \qquad - \text{c}\_{\text{b}}\cdot\sigma\cdot\text{F}\_{\text{al}^{\text{l}}}^{\text{b}}(\text{T}\_{\text{s-bottom}}^{4} - \text{T}\_{\text{grd}}^{4})\\ \qquad \qquad \qquad \qquad - \frac{1}{\text{R}\_{\text{vfl}}}\cdot(\text{T}\_{\text{s-bottom}} - \text{T}\_{\text{s-top}}) = 0 \end{array} \tag{3}$$

In the energy balance equations, the meanings of the different symbols are detailed in the nomenclature, using subscripts t and b to denote top and bottom wing surface, respectively. The convective heat transfer coefficient, denoted h is obtained from calculating the Nusselt Number on a flat plate in laminar regime as detailed in [33] using Mars atmospheric properties and assuming that the dominant gas is *CO*2. These properties were obtained from [34]. The atmospheric density was determined by using ideal gas law at 8.5 mbar of atmospheric pressure.

#### *2.2. Sky, Ground, and Ambient Temperatures*

As mentioned previously, this study includes the summer and winter seasons. As such, the sky, ground, and ambient temperatures for both seasons are needed. The information for a summer day corresponds to a solar longitude of 140◦ and for a winter day corresponds to a solar longitude of 280◦. Figure 2 shows the ambient temperatures based on the season. The ambient temperature for summer is higher than winter and has a larger variation over the day. The summer daily variation has a range of almost 50 ◦C, while the winter only has a daily variation of nearly 15 ◦C [31]. Figure 3 shows the sky and ground temperatures as reported by Matz et al. [32] at the same location on the same day. The ground temperature follows a similar trend where the summer temperatures are higher and vary more throughout the day compared to the winter ground temperatures. The sky temperatures, however, have the opposite trend. The winter sky temperature is higher and varies more compared to the summer sky temperatures.

**Figure 2.** Ambient temperature for (**a**) summer and (**b**) winter in the Martian atmosphere [31].

**Figure 3.** Sky and ground temperatures for (**a**) summer and (**b**) winter in the Martian atmosphere [32].

#### *2.3. Solar Energy Absorption*

Being inspired by migrating birds' flight on Earth, different color combinations are investigated. It is assumed that the top and bottom surfaces have different colors and therefore different absorptivity values. The different color configurations considered during this study are shown in Figure 4. The four configurations are on the top and bottom surfaces, respectively, black–white, black–black, white–black, and white–white. The solar irradiance that acts on the top surface of the wings includes direct beam, diffuse, and albedo radiations. The wing surface exposure to different radiation fluxes is governed by the view factors to the sky and the ground. The absorbed energy on the two wing surfaces can be written as [32]:

$$\mathbf{P\_{solar-top}} = \boldsymbol{\alpha\_{t}} \cdot \mathbf{G\_{b}} \cdot \cos \theta + \boldsymbol{\alpha\_{t}} \cdot \mathbf{F\_{d}^{t}} \cdot \mathbf{G\_{dh}} + \boldsymbol{\alpha\_{t}} \cdot \mathbf{a\_{l}} \cdot \mathbf{F\_{u}^{t}} \cdot (\mathbf{G\_{b}} \mathbf{h} + \mathbf{G\_{dh}}) \tag{4}$$

$$\mathbf{P\_{solar-bottom}} = \boldsymbol{\alpha\_{b}} \cdot \mathbf{F\_{q}^{b}} \cdot \mathbf{G\_{dh}} + \boldsymbol{\alpha\_{b}} \cdot \mathbf{al} \cdot \mathbf{F\_{al}^{b}} \cdot (\mathbf{G\_{bh}} + \mathbf{G\_{dh}}) \tag{5}$$

**Figure 4.** Wing color configurations in this study.

The view factors for the top and bottom surfaces are needed and are expressed by Equations (6)–(9), where β is the angle of attack. The view factor for the top surface to the sky is *F<sup>t</sup> <sup>d</sup>* while the view factor for the top surface to the ground is F*<sup>t</sup> al*. For the bottom surface, the view factor for the bottom to the sky F*b <sup>d</sup>* is identical to <sup>F</sup>*<sup>t</sup> al*. Likewise, the view factor from the bottom surface to the ground <sup>F</sup>*<sup>b</sup> al* is identical to *Ft d* . In this case, the drone is considered to be flying very close to the ground because some of the atmospheric properties used were measured by VL1 at ground level. It is noted that at a zero-degree angle of attack Equations (2) and (3) are decoupled because the view factor from the top surface to the sky is 1 and the view factor for the top surface to the ground is zero.

$$F\_d^t = \frac{1 + \cos \beta}{2} \tag{6}$$

$$\mathbf{F}\_{al}^{t} = \frac{1-\cos\beta}{2} \tag{7}$$

$$\mathbf{F}\_d^b = \mathbf{F}\_{al}^t \tag{8}$$

$$F\_{al}^b = F\_d^t \tag{9}$$

The absorptivity values for α*<sup>t</sup>* and α*<sup>b</sup>* used in this study are 0.88 and 0.23, which correspond to the absorptance of anodize black paint and biphenyl white solid paint, respectively. Both are often used for space applications [35]. Depending on several variables, such as the solar angle of incidence shown in Figure 5, the optical depth of the Martian atmosphere, and the season, the solar irradiance changes. NASA researchers from Lewis Research Center in Ohio have derived a set of equations to determine the solar irradiance on Mars [36]. One of the most crucial elements for determining the solar irradiance is the position of the wing with respect to the sun. The solar incidence angle θ is used to describe this relation and can be determined from Equation (10), as follows [37]:

$$\begin{aligned} \cos\theta &= \sin\varphi \sin\delta \cos\beta - \sin\delta \cos\varphi \cos\lambda \sin\beta + \cos\varphi \cos\delta \cos\omega \cos\beta \\ &+ \cos\delta \sin\varphi \cos\lambda \cos\omega \sin\beta + \cos\delta \sin\lambda \sin\omega \sin\beta \end{aligned} \tag{10}$$

**Figure 5.** Solar angles.

To determine the solar incidence angle, positional information about the drone must be known, such as the latitude ϕ, declination angle δ, hour angle ω, angle of attack β, and the wing azimuth angle λ. The wing azimuth angle can be considered as the direction of flight. In this study, the drone is assumed to be flying directly southward. The hour angle only considers daylight hours from 7:00 to 17:00. Additionally, the hour angle refers to solar Martian time (0:00–24:00). The declination angle depends on the season which is described by the areocentric longitude *Ls*. The latitude selected for the study is the location of the rover which is where the atmospheric properties are obtained. The values for latitude, declination angle, tilt angle, and wing azimuth angle that are used in this study are shown in Table 1 for both summer and winter seasons.



After determining the solar incidence angle, the solar irradiance can then be determined. The direct beam irradiance Gb, diffuse irradiance on a horizontal surface *Gdh*, and the global irradiance on a horizontal surface Gh which are used to determine the beam, diffuse, and albedo irradiance, respectively, can be determined from Equations (11)–(13) [36]. The optical depth τ varies depending on the date which is represented by the areocentric longitude Ls. The optical depth is a result of the atmospheric dust particles in the atmosphere [38]. The values for the solar irradiance variables that are used for this study are shown in Table 2. The normalized net flux function f(θ, τ) values [36] used during this study that correspond to the solar incidence angle and optical depth are presented in Table 3.

$$\mathbf{G}\_{\mathsf{b}} = \mathbf{G}\_{\mathsf{cb}} \cdot \mathbf{e}^{\left(\frac{-\pi}{\cos \mathsf{B}}\right)} \tag{11}$$

$$G\_{dl\mathfrak{h}} = G\_{\mathfrak{h}} - G\_{\mathfrak{l}\mathfrak{h}} \tag{12}$$

$$\mathbf{G}\_{\rm h} = \frac{\mathbf{G}\_{\rm cb} \cdot \cos \theta \cdot \mathbf{f}(\theta, \tau)}{1 - \mathbf{al}} \tag{13}$$


**Table 2.** Solar irradiance variables [36].



**Table 3.** Normalized net flux function [36].

The solar irradiance at the top of the atmosphere Gob, which can be seen in Equation (14) [36], is needed along with the optical depth τ, solar incidence angle θ, and the normalized net flux function f(θ, τ) to determine the direct beam irradiance and the global irradiance. Mars eccentricity b is around 0.093377, as can be seen in Table 2. The diffuse irradiance can be determined by subtracting the direct beam irradiance on a horizontal surface *Gbh* from the global irradiance on a horizontal surface Gh. The equation for direct beam irradiance on a horizontal surface can be expressed as shown in Equation (15) [36].

$$\mathbf{G\_{ob}} = \frac{590 \cdot \left[1 + \mathbf{b} \cdot \cos\left(\mathbf{L\_s} - 248^\circ\right)\right]^2}{\left(1 - \mathbf{b}^2\right)^2} \tag{14}$$

$$\mathbf{G\_{bh}} = \mathbf{G\_{ab}} \cdot \cos \theta \cdot \mathbf{e^{\left(\frac{-\mathbf{r}}{\cos \theta}\right)}} \tag{15}$$

#### *2.4. Radiative Heat Loss*

Equations (16) and (17) represent the radiative heat loss for the top and bottom wing surfaces, respectively. The top and bottom wing surfaces will differ depending on the emissivity ε and the shape

factors. For a zero-degree angle of attack, the top wing surface only loses energy to the sky while the bottom wing surface only loses energy to the ground. The values used for emissivity are based on anodize black paint and biphenyl white solid paint, respectively. For both white and black colors, the emissivity is 0.88 [24]. The sky and ground temperatures can be found in Figure 3a,b, respectively.

$$P\_{\rm radiative-top} = \varepsilon\_{\rm t} \cdot \sigma \cdot \mathrm{F\_{dt}} \cdot \left( T\_{\rm s-top}^{4} - T\_{\rm sky}^{4} \right) + \varepsilon\_{\rm t} \cdot \sigma \cdot \mathrm{F\_{alt}} \cdot \left( T\_{\rm s-top}^{4} - T\_{\rm grid}^{4} \right) \tag{16}$$

$$P\_{\text{radiative-bottom}} = \varepsilon\_{\text{b}} \cdot \sigma \cdot \text{F}\_{\text{db}} \cdot (\text{T}\_{\text{s-bottom}}^4 - \text{T}\_{\text{sky}}^4) + \varepsilon\_{\text{b}} \cdot \sigma \cdot \text{F}\_{\text{alb}} \cdot (\text{T}\_{\text{s-bottom}}^4 - \text{T}\_{\text{grd}}^4) \tag{17}$$

#### *2.5. Thermal Boundary Layer Analysis*

For simplicity, the wing curvature is ignored, and a flat plate is considered for the boundary layer analysis. Figure 6 represents a schematic view of the boundary layer over a flat plat [25].

**Figure 6.** Considered boundary layer over a flat plate [25].

The thermal boundary layer effects on the skin friction drag are determined using Equation (18) [39], where the wingspan ws, is 3.5 m and the chord length, L, is 0.22 m [7]. To determine the viscosity μ, Sutherland's formula, which can be found in Equation (19), is used [40]. T0 and μ<sup>0</sup> are the reference temperature and viscosity, respectively, which are 293.14 K and 1.48 <sup>×</sup> <sup>10</sup>−<sup>5</sup> kg <sup>m</sup>·<sup>s</sup> for the Martian atmosphere [41]. The Sutherland constant C, for CO2 is 240 K [41] and Ts is the wing surface temperature. The density ρ is determined using ideal gas law considering CO2 in Equation (20). The ideal gas constant for the Martian atmosphere R, is 192.1 <sup>J</sup> kg·<sup>K</sup> [30]. The pressure P is found from the information collected by VL1.

$$\mathbf{D} = 0.664 \cdot \mathbf{w} \mathbf{s} \cdot \sqrt{\mathbf{p} \cdot \mathbf{u} \cdot \mathbf{u}^3 \cdot \mathbf{L}} \tag{18}$$

$$
\mu = \mu\_0 \cdot \left(\frac{T\_0 + C}{T\_s + C}\right) \cdot \left(\frac{T\_s}{T\_0}\right)^{3/2} \tag{19}
$$

$$\rho = \frac{\text{P}}{192.1 \text{R} \cdot \text{T}\_{\text{s}}} \tag{20}$$

#### **3. Wing Surface Analysis**

In this section, the whole wing (top and bottom surfaces) is considered for the study. The wing structure will be considered as a hollow wing composed of a hollow center with a top layer and a bottom layer of material. The model of the overall wing can be seen in Figure 7a,b. From Figure 7b, it can be observed that the conductive heat transfer will cover three layers of material; the first layer of material, the interior which is air (CO2 for the Martian atmosphere), and the third layer of material.

**Figure 7.** (**a**) Hollow airfoil and (**b**) conduction model.

The expression for heat conduction is given in Equation (21), where the conduction is determined from the temperature difference between the top surface Ts-top, and the bottom surface Ts-bottom divided by the effective resistance of all materials. The effective resistance is presented in Equation (22) where Lsurface is the depth of the wing material, Lair is the average depth of the hollow part of the wing, Ksurface denotes the thermal conductivity of the material, Kair represents the thermal conductivity of the Martian atmosphere, Asurface is the area of the material, and Aair denotes the area of the hollow part.

$$P\_{\text{conductive}} = \frac{T\_{\text{s-top}} - T\_{\text{s-bottom}}}{R\_{\text{eff}}} \tag{21}$$

$$\mathbf{R\_{eff}} = \frac{2\mathbf{L\_{surface}}}{\mathbf{K\_{surface}}\mathbf{A\_{surface}}} + \frac{\mathbf{L\_{air}}}{\mathbf{K\_{air}}\mathbf{A\_{air}}}\tag{22}$$

Air (CO2) acts as an insulator in this case which prevents conductive heat transfer from occurring between the top surface and the bottom surface. The surface temperature for the top and bottom surfaces with and without considering conductive heat fluxes are shown in Figure 8 for summer.

**Figure 8.** Surface temperature with and without conduction for (**a**) top surface and (**b**) bottom surface.

Air has a very low thermal conductivity, which will result in a very high thermal resistance. The heat conducted through the wing can be considered negligible compared to the total energy balance. This allows Equations (2) and (3) to be simplified into Equations (23) and (24). The two new energy balance equations become uncoupled without the conductive heat transfer and are used for the rest of this study.

$$\begin{aligned} \mathbf{^0a\_t}\mathbf{^0b\_b}\cos z + \mathbf{^0\_t}\mathbf{^0\_d}\mathbf{^0\_d}\mathbf{^0\_d} + \mathbf{^0\_t}\mathbf{^1d}\mathbf{^1\_{sl}}(\mathbf{^0\_{bh}} + \mathbf{^0\_{dh}}) - \mathbf{h^{\cdot}}(\mathbf{T\_{s\cdot t\text{op}}} - \mathbf{T\_{amb}}) - \boldsymbol{\varepsilon\_t}\cdot\boldsymbol{\sigma}\cdot\mathbf{F^t\_d} \\ \cdot (\mathbf{T^4\_{s\cdot t\text{op}}} - \mathbf{T^4\_{s\text{ky}}}) - \boldsymbol{\varepsilon\_t}\cdot\boldsymbol{\sigma}\cdot\mathbf{F^t\_{al}}\cdot(\mathbf{T^4\_{s\cdot t\text{op}}} - \mathbf{T^4\_{g\text{rd}}}) &= 0 \end{aligned} \tag{23}$$

$$\begin{aligned} \mathbf{a}\_{\rm b} \cdot \mathbf{F}\_{\rm d}^{\rm b} \cdot \mathbf{G}\_{\rm dh} + \boldsymbol{\alpha}\_{\rm b} \cdot \mathbf{al} \cdot \mathbf{F}\_{\rm al}^{\rm b} \cdot (\mathbf{G}\_{\rm bh} + \mathbf{G}\_{\rm dh}) - \mathbf{h} \cdot (\mathbf{T}\_{\rm s\text{-bottom}} - \mathbf{T}\_{\rm amb}) - \boldsymbol{\varepsilon}\_{\rm b} \cdot \boldsymbol{\sigma} \cdot \mathbf{F}\_{\rm d}^{\rm b} \\ \cdot (\mathbf{T}\_{\rm s\text{-bottom}}^{4} - \mathbf{T}\_{\rm sky}^{4}) - \boldsymbol{\varepsilon}\_{\rm b} \cdot \boldsymbol{\sigma} \cdot \mathbf{F}\_{\rm al}^{\rm b} (\mathbf{T}\_{\rm s\text{-bottom}}^{4} - \mathbf{T}\_{\rm grd}^{4}) = 0 \end{aligned} \tag{24}$$

When considering an angle of attack of zero degrees, the equations are completely decoupled and thus can be solved individually. Therefore, the effects of the conduction are not considered in further analysis. As can be seen in Figure 9a,b, the solar absorption for the top surface (both white and black colors) are higher than the bottom surface. This is due to the direct beam irradiance that is only experienced by the top surface. Additionally, the solar absorption during summer is higher than during the winter due to the high solar intensity that occurs during the summer period of the Martian solar cycle.

**Figure 9.** Surface solar absorption for (**a**) summer and (**b**) winter.

The two new energy balance Equations (23) and (24) can then be used to determine the surface temperature for both the bottom and top wing surfaces for black and white wing colors as presented in Figure 10. During the summer and winter seasons, the wing surface experiences the highest temperatures when a black color is applied to the top surface. The lowest wing temperatures occur when white paint is used on the top surface. As expected, due to the high solar absorption during summer, the overall temperatures on the wing surfaces are higher during summer compared to winter.

**Figure 10.** Surface temperature for (**a**) summer and (**b**) winter.

Using the surface temperature, Equations (18)–(22) can be solved to find the viscosity, density, and drag. The viscosity for both wing surfaces and for black and white colors are plotted in Figure 11a,b, respectively. The viscosity follows a similar trend to the surface temperature. The higher surface temperature wing color configuration results in a higher viscosity along the boundary. The density, which can be seen in Figure 12a, is lower along the boundary layer when the surface temperature is higher. The total drag of the combined wing color configurations is plotted in Figure 13a,b for summer and winter, respectively. The four different wing color combinations are as follows: black on top and black on bottom (black–black), black on top and white on bottom (black–white), white on top and black on bottom (white–black), and white on top and white on bottom (white–white). It follows from Figure 13a, during the summer, that the highest drag color configuration is the white–white, and the lowest drag is the black–black. During the winter, it is shown in Figure 13b that the highest drag is from the color configuration black–black and the lowest drag is from the color configuration white–white.

**Figure 11.** Surface viscosity for (**a**) summer and (**b**) winter.

**Figure 12.** Surface density for (**a**) summer and (**b**) winter.

**Figure 13.** Surface drag for (**a**) summer and (**b**) winter in the Martian atmosphere.

Comparing the relative change between the highest drag color configuration and the lowest drag color configuration results in Figure 13a,b. During the summer the highest drag color configuration is white–white and the lowest is the black–black color configuration. During the winter, the trend is reversed with the highest drag color configuration being black–black and the lowest color configuration white–white.

During the summer season, the optimal wing color configuration is black–black which results in a decrease of drag by a maximum of 0.125% compared to the white–white wing color configuration. During the winter season, the optimal wing color configuration is white–white which has a maximum drag reduction of 0.15% compared to the black–black configuration, as shown in Figure 14.

**Figure 14.** Relative change for (**a**) summer and (**b**) winter in the Martian atmosphere.

Until now, the results have been generated using a zero-degree angle of attack (AoA). To understand the effects of the angle of attack, the drag as a function of angle of attack at solar noon is plotted in Figure 15. The change in drag during the summer period is near constant over a range of 10 degrees AoA. During the winter, the maximum change in drag occurs for the white–black color configuration with a change percentage of 0.001%. It can be concluded that the change in angle of attack is negligible to determining the drag force. This led to using a constant zero-degree angle of attack for the study.

**Figure 15.** Drag as a function of angle of attack for (**a**) summer and (**b**) winter.

#### **4. 3D Aerodynamic Analysis**

The skin friction drag has already been determined based on the boundary layer analysis in Section 3. In this section, the aerodynamic performance for the whole wing will be analyzed by determining the total drag. The drag will be determined using an aerodynamic analysis software named XFLR5. The total drag includes the form drag, the induced drag, and the parasitic drag. The software used to analyze the 3D wing geometry can determine the drag coefficients. The study considers a constant zero-degree angle of attack under similar conditions as the boundary layer analysis from Section 3 and the two results are compared. Four analysis methods, namely, Lifting Line Theory (LLT), Horseshoe Vortex (VLM1), Ring Vortex (VLM2), and 3D panel method are used in order to estimate the total drag of an Albatross-inspired drone in the Martian atmosphere. All of these methods consider the assumption of inviscid flow.

The LLT method uses the assumption that the lift coefficient as a function of the angle of attack is linear. One major assumption for the LLT method is that the surfaces are in the X-Y plane and that the dihedral angles and the sweep are not needed for the lift distribution calculations. The lift of the wing for this theory is determined from the incremental vortices shed along the trail span of the wing along the freestream direction [42].

The VLM1 and VLM2 methods allow for more freedom in the choice of wing geometry. This includes winglets, high dihedral angles, low aspect ratio, and sweep. For the VLM methods, the general idea is to model the perturbation of the wing planform using a vortices' distribution sum. The lift coefficient is computed by integrated surface forces, such as the moment coefficients. The VLM methods calculate the lift distribution, induced angles, and induced drag. Because of this, the methods are independent of the wing speed. The VLM methods assume small angles of attack, so stall angles should be avoided when computing results. The main difference between the VLM1 and VLM2 methods is that the VLM1 method does not consider side slip [42].

The 3D panel method is able to account for the thickness of the wing unlike the VLM methods that only consider the mean camber line. Because of this, the 3D panel method can allow for an understanding of the center of pressure distribution over the top and bottom surfaces of the wing. The 3D panel method uses an approach that sum the doublets and sources distributed of the surface of the wing from the perturbations generated. The 3D panel method can be used to polish the results of the other three methods, so this is often considered to be the most accurate method. The 3D panel method will be considered as the standard for this study that the other three methods will be compared to [42].

As mentioned in the previous section, the 3D aerodynamic analysis methods all have inviscid assumptions. To account for the viscosity in the 3D analysis, a 2D viscous analysis is performed for the selected airfoil to obtain 2D viscous information. The 2D information is then interpolated and then incorporated into the 3D analysis and this is how the viscous effects are included in the analysis. The 2D boundary layer problem is solved using a software called XFoil by using an iterative method along with the inviscid velocity field as an input. The method provides a viscous velocity from the boundary layer solution which is used as an input for the potential flow solver. This iterative method is called the "Interactive Boundary Layer" [42]. The wing shape considered for the 3D analysis is that of an Albatross. Figure 16 shows the analyzed wing shape and Table 4 lists the geometric properties.

**Figure 16.** 3D wing shape.

**Table 4.** Wing geometry of the drone.


To solve the 2D and 3D analyses, certain properties are needed, namely, the Reynolds number and wing velocity. The wing velocity and the Reynolds number for the 3D analysis are the same as for the boundary layer analysis. The Reynolds number is calculated based on the viscosity and density found in Section 3. The total drag varies depending on the analysis method. As mentioned previously, the 3D panel method is considered to be the most accurate method and the other three methods are compared to the 3D panel method to determine their accuracy. The number of panels on the wing is determined in order for the different methods of analysis to converge. Figure 17 compares the total drag for the different analysis methods for both summer and winter. Inspecting the plotted curves in these figures, it is clear that VLM1 and VLM2 methods have comparable results to the 3D panel method. The LLT method varies greatly compared to the other three methods. In Tables 5 and 6, the percent difference between the different methods compared to 3D panel method is shown. The percent difference between VLM1 and VLM2 compared to the 3D panel method is 3~4%, which shows good agreement. However, the LLT method compared to the 3D panel method has a greater than 50% error, which shows that the LLT method is not capable of handling the complex geometry of the Albatross-like wing shape.

**Figure 17.** Total drag comparison of the different analysis method for black–black color configuration for (**a**) summer and (**b**) winter.


**Table 5.** Drag percent difference compared to 3D panel method during the summer.


**Table 6.** Drag percent difference compared to 3D panel method during the winter.

The total drag forces for the wing at a constant angle of attack and at different times throughout the day is calculated using the 3D panel method. The total drag includes the form drag, the induced drag, and the parasitic drag. Figure 18 shows the total drag dependent on the hour for the 3D panel method. Compared to the drag determined using the boundary layer analysis, the drag calculated using the 3D panel method is much larger, by a factor of 5 on average. During the winter, the drag is much higher than during the summer. This remains consistent with the trend found for the skin friction drag. Additionally, the black–black wing color configuration is the most optimal case that results in the least amount of drag for both summer and winter. The percent change between the most optimal wing color configuration (black–black) and the least optimal color configuration (white–white) is much more significant when considering the total drag instead of just the skin drag. For summer the percent change is 12.6% and for the winter the change is 6.8%. This demonstrates that the effects of the skin drag are significantly less compared to the total drag.

**Figure 18.** Total drag as a function of time using 3D panel method for (**a**) summer and (**b**) winter.

#### **5. Thermoelectric Energy Generators in Martian Atmosphere**

The basis of a thermoelectric generator is the ability to convert heat energy into electrical power when there exists a temperature difference between two surfaces. Most TEGs are created using pairs of thermocouples that are connected in series (electrical) and parallel (thermal) in between two plates. A representation of a typical TEG is presented in Figure 19. The power produced from the TEG is the difference between the heat supply Ph and heat removal Pc which is also equal to the product of the voltage across the external load resistance V and the current I. The energy balance can be expressed as:

$$\mathbf{P}\_{\rm TEG} = \mathbf{P}\_{\rm h} - \mathbf{P}\_{\rm c} = \mathbf{V} \cdot \mathbf{I} \tag{25}$$

**Figure 19.** Schematic of a thermoelectric generator.

The derived equation for the electrical harvested power is written as [15]:

$$\mathbf{P\_{TEG}} = \mathbf{R\_L} \cdot \frac{\mathbf{a^2 \cdot (T\_h - T\_c)}^2}{\left(\mathbf{R\_i - R\_L}\right)^2} \tag{26}$$

It follows from Equation (26) that the electrical power is dependent on the Seebeck coefficient α, the internal resistance Ri, the external resistance RL, and the temperature difference between the hot surface Th and the cold surface Tc. Depending on the specific TEG and the manufacturer, the coefficients will differ. It is often assumed that the internal properties of a TEG device are constant and independent of temperature. However, it has been demonstrated through experimental data by Hsu et al. [28] that the internal properties vary as a dependent on the hot side temperature. For the experiment, Hsu et al. considered a TEG module (TMH400302055, Wise Life Technology, Taiwan) and varied the hot side

temperature to determine the effects on the Seebeck coefficient and the internal resistance. Expressions for both the internal resistance and the Seebeck coefficient are determined by [15]. The expressions are validated to the experimental data by Hsu et al. [28]. The expressions for the Seebeck coefficient and the internal resistance are given by:

$$\mathbf{a}(\mathbf{T}) = \mathbf{a}\mathbf{1}^2 + \mathbf{a}\mathbf{2}\mathbf{T} + \mathbf{a}\mathbf{3} \tag{27}$$

$$\mathbf{R}\_{\mathbf{i}}(\mathbf{T}) = \mathbf{b}\_{1}\mathbf{T} + \mathbf{b}\_{2} \tag{28}$$

The temperature range that is experienced by the wing surfaces on Mars are 180–300 K. In Figure 20, the internal properties for the TEG under consideration are determined by using the expressions above for the Seebeck coefficient and the internal resistance.

**Figure 20.** Internal properties of a TEG device over a temperature range for (**a**) Seebeck coefficient and (**b**) internal resistance.

An analysis is performed to determine the harvested power over a range of values for the external resistance, Seebeck coefficient, and internal resistance. These coefficients were determined over a range of temperatures from 180–300 K. According to Figure 20, the Seebeck coefficient ranges from 0.064–0.072 (V/K) and the internal resistance ranges from 0.85–1.15 (Ω). However, the set of temperatures were obtained using data from one specific day each during the summer and winter. Consequently, to account for a certain amount of variation, a Seebeck coefficient ranging from 0.05–0.08 (V/K) and an internal resistance ranging from 0.8–1.3 (Ω) will be used for the parametric study and performance analysis. Figure 21 shows a 3D plot representing the harvested power as a function of the external resistance and the Seebeck coefficient with a fixed value for the internal resistance during summer for each color configuration. Figure 21a–d are for the color configurations black–white, black–black, white–black, and white–white, respectively. Each color configuration has a different time of day for which the maximum power occurs. This is because the maximum power is occurring when the temperature difference between the top and bottom surfaces is greatest. For the maximum harvested power, the time of day for each color configuration (a–d) is 10:00, 10:00, 15:00, and 15:00 solar time.

**Figure 21.** Harvested power as a function of the external resistance and Seebeck coefficient when Ri = 0.8 Ω for (**a**) Black–white (10:00), (**b**) Black–black (10:00), (**c**) White–black (15:00), and (**d**) White–white (15:00).

The optimal power generation for winter and summer is shown in Figure 22. The corresponding time of day for each color configuration is also shown. The black–white color configuration has the largest power generation while the white–white color configuration has the lowest. During the earlier part of the day, when the top surface of the wing is black, the TEG experiences higher power generation. We can derive the optimal power generation QP-Optimal from Equation (26) using the expression:

$$\frac{\partial P\_{TEG}}{\partial \mathcal{R}\_L} = 0 \tag{29}$$

**Figure 22.** Optimal harvested power as a function of internal resistance for (**a**) summer and (**b**) winter when the Seebeck coefficient is equal to 0.08.

The resulting optimal power generation equation is:

$$P\_{\text{TEG-Optimal}} = \frac{\alpha^2 \cdot (\Delta \text{T})^2}{4 \cdot \text{R}\_i} \tag{30}$$

In Table 7, the optimal power generation along with the maximum temperature difference and the time of day for each color configuration for both summer and winter is given. The optimal wing color configuration that generates the most harvested power is the black–white configuration for both summer and winter. This is due to having the highest temperature difference. The white–white configuration has the lowest power generation and the lowest temperature difference. Overall, the TEG performs better during the summer than during the winter.

**Table 7.** Optimal power comparison (Time, ΔTmax, Power).


#### **6. Conclusions**

To enhance the design of an efficient drone capable of flying in the Martian atmosphere, a parametric study was performed to analyze the thermal effects of wing color on the heated boundary layer of a wing. This study provided information on the optimal wing color configuration during different seasons on Mars considering the color configurations of black–black, black–white, white–black, and white–white. A model based on an energy balance, and including wing color configurations, sky, ground, and ambient temperatures, solar energy absorption, convective heat transfer, and radiative heat loss, was developed and discussed. The adiabatic case, which only includes the top wing surface and drag effects, was analyzed. The effects of conduction and angle of attack on the total energy balance were studied. They both had a negligible effect on the thermal energy balance. The drag effects considering the color configurations on both sides of the wing were investigated. The first method analyzed the boundary layer and the skin friction drag, while the second method included a 3D analysis that investigated the total drag. The skin drag did not significantly change as a result of the color configuration. However, the total drag experienced a change percent of 12.8% during summer and 6.8% during winter. This shows that the skin drag does not play a significant role in the total drag on Mars due to the low atmospheric pressure. The study performed on the thermoelectric generator yielded interesting results. Due to the high temperature differences on the top and bottom wing surfaces that occur on Mars, TEGs are capable of high power generation. The black–white configuration had the largest temperature difference and resulted in the greatest power generation for both respective seasons, summer and winter.

**Author Contributions:** Conceptualization, D.R., S.B.A., S.J., and A.A.; methodology, D.R., S.B.A., and A.A.; software, D.R.; formal analysis, D.R. and S.B.A.; investigation, D.R.; resources, A.A.; writing—original draft preparation, D.R.; writing—review and editing, S.B.A., S.J., and A.A.; supervision, S.B.A., S.J., and A.A.; project administration, A.A.; funding acquisition, S.J. and A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by New Mexico Consortium and Los Alamos National Laboratory.

**Acknowledgments:** The authors D. Rice and A. Abdelkefi would like to acknowledge the financial support from New Mexico Space Grant Consortium. Additionally, the authors would like to thank Mostafa Hassanalian for his fruitful discussions in this subject.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

A = surface area al = Albedo b = Mars eccentricity C = Sutherland- s constant D = Drag force f(z, τ)= Normalized net flux Fdt= View factor top-side to the sky Falt= View factor top-side to the ground Fdb= View factor bottom-side to the sky Falb= View factor bottom side to the ground Gb= Direct beam irradiance on surface Gbh= Direct beam irradiance horizontal surface Gdh= Diffuse irradiance on horizontal surface Gh= Global irradiance on horizontal surface Gob= Beam irradiance at top of the Atmosphere h = heat transfer coefficient I = Current K = Atmospheric thermal conductivity L = Chord length Ls= Areocentric Longitude Nu = Nusselt number P = Atmospheric Pressure Pr = Prandtl number Ph and Pc = Hot and cold sources PConductive = Conductive heat transfer Pconvective= Convective heat loss Pradiative= Radiative heat loss Psolar= Solar energy absorbed PTEG= Electrical Power R = Ideal gas constant R*<sup>i</sup>* = Internal resistance Re = Reynolds number Reff = Effective conductive resistance of the wing RL = External resistance T0= Reference temperature for Sutherlands Tamb= Ambient temperature Tgrd= Ground temperature Ts= Wing surface temperature Tsky= Sky temperature u = velocity V = Voltage across external load resistance ws = Wing span α = Seebeck coefficient αt= Top surface absorptance αb= Bottom surface absorptance β = Angle of attack δ = Declination angle εt= Top surface emissivity εb= Bottom surface emissivity θ = Solar incidence angle

λ = Wing azimuth angle (direction)

*Drones* **2020**, *4*, 43


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
