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

: Thermal e ﬀ ects 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 e ﬀ ects of wing color: the drag reduction and the thermoelectric generation of power. According to its color, each wing side has a certain temperature a ﬀ ecting the drag. Investigations of various conﬁgurations have shown that the thermal e ﬀ ect on the wing boundary layer skin drag is insigniﬁcant because of the low atmospheric pressure. However, the total drag varies as much as 12.8% between the highest performing wing color conﬁguration and the lowest performing conﬁguration. Additionally, the large temperature di ﬀ erences between the top and the bottom wing surfaces show great potential for thermoelectric power generation. The maximum temperature di ﬀ erences 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 beneﬁt from increased endurance through optimizing the wing color conﬁguration. The e ﬀ ects of conduction and angle of attack on the total energy balance were studied. They had a negligible e ﬀ ect on the thermal energy balance. The drag ﬀ the color conﬁgurations on both sides of the wing were investigated. The ﬁrst method analyzed the boundary layer and the skin friction drag, while the second method a 3D analysis that investigated the total drag. The skin drag did not signiﬁcantly change as a result of the color conﬁguration. 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 signiﬁcant 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 di ﬀ erences on the top and bottom wing surfaces that occur on Mars, TEGs are capable of high power generation. The black–white conﬁguration had the largest temperature di ﬀ erence and resulted in the greatest power generation for both respective seasons, summer and winter.


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

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.

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: 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]. 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].
α · G · cos ϴ + α · F · G + α · al · F · (G + G ) − h · T − T − ε · σ · F · T − T − ε · σ · F · T − T − · T − T =0 (2) 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].
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.

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.
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 . These properties were obtained from [34]. The atmospheric density was determined by using ideal gas law at 8.5 mbar of atmospheric pressure.

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.

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]: P = α · G · cos ϴ + α · F · G + α · al · F · (G + G ) (4) P = α · F · G + α · al · F · (G + G ) 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 while the view factor for the top surface to the ground is F . For the bottom surface, the view factor for the bottom to the sky F is identical to F . Likewise, the view factor from the bottom surface to the ground F is identical to . 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.

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]:

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]: P = α · G · cos ϴ + α · F · G + α · al · F · (G + G ) (4) P = α · F · G + α · al · F · (G + G ) 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 while the view factor for the top surface to the ground is F . For the bottom surface, the view factor for the bottom to the sky F is identical to F . Likewise, the view factor from the bottom surface to the ground F is identical to . 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. 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 t d while the view factor for the top surface to the ground is F t al . For the bottom surface, the view factor for the bottom to the sky F b d is identical to F t al . Likewise, the view factor from the bottom surface to the ground F b al is identical to F t 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.
Drones 2020, 4, 43 The absorptivity values for α t and α b 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]: cos θ = sin ϕ sin δ cos β − sin δ cos ϕ cos λ sin β + cos ϕ cos δ cos ω cos β + cos δ sin ϕ cos λ cos ω sin β + cos δ sin λ sin ω sin β (10) The absorptivity values for α and α 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]: 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 . 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. 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 L s . 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. Table 1. Location parameters [36]. After determining the solar incidence angle, the solar irradiance can then be determined. The direct beam irradiance G b , diffuse irradiance on a horizontal surface G dh , and the global irradiance on a horizontal surface G h 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 L s . 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. The solar irradiance at the top of the atmosphere G ob , 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 G bh from the global irradiance on a horizontal surface G h . The equation for direct beam irradiance on a horizontal surface can be expressed as shown in Equation (15) [36].

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.

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]. 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]. T 0 and µ 0 are the reference temperature and viscosity, respectively, which are 293.14 K and 1.48 × 10 −5 kg m·s for the Martian atmosphere [41]. The Sutherland constant C, for CO 2 is 240 K [41] and T s is the wing surface temperature. The density ρ is determined using ideal gas law considering CO 2 in Equation (20). The ideal gas constant for the Martian atmosphere R, is 192.1 J kg·K [30]. The pressure P is found from the information collected by VL1.

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 (CO 2 for the Martian atmosphere), and the third layer of material. 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. 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 Tsbottom 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.
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. The expression for heat conduction is given in Equation (21), where the conduction is determined from the temperature difference between the top surface T s-top, and the bottom surface T s-bottom divided by the effective resistance of all materials. The effective resistance is presented in Equation (22) where L surface is the depth of the wing material, L air is the average depth of the hollow part of the wing, K surface denotes the thermal conductivity of the material, K air represents the thermal conductivity of the Martian atmosphere, A surface is the area of the material, and A air denotes the area of the hollow part.
Air (CO 2 ) 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. 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.
When considering an angle of attack of zero degrees, the equations are completely decoupled 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.
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.
balance equations become uncoupled without the conductive heat transfer and are used for the rest of this study. α · G · cos + α · F · G + α · al · F · (G + G ) − h · T − T − ε · σ · F · T − T − ε · σ · F · T − T = 0 (23) 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 Figures 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. 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. 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.
(a) (b) 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 Figures  13a,b for summer and winter, respectively. The four different wing color combinations are as follows: 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 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.        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. Comparing the relative change between the highest drag color configuration and the lowest drag color configuration results in Figures 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. 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. 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.
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.

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.  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.  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.
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.  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. 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 (whitewhite) 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.

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.

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 P h and heat removal P c 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: A representation of a typical TEG is presented in Figure 19. The power produced from the TEG is the difference between the heat supply P and heat removal P 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: P = P − P = V · I (25) Figure 19. Schematic of a thermoelectric generator.
The derived equation for the electrical harvested power is written as [15]: It follows from equation (26) that the electrical power is dependent on the Seebeck coefficient α, the internal resistance R , the external resistance R , and the temperature difference between the hot surface T and the cold surface T . 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. [43] 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. [43]. The expressions for the Seebeck coefficient and the internal resistance are given by: The derived equation for the electrical harvested power is written as [15]: It follows from Equation (26) that the electrical power is dependent on the Seebeck coefficient α, the internal resistance R i , the external resistance R L , and the temperature difference between the hot surface T h and the cold surface T c . 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: 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. 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. Figures 21a-d are for the color configurations blackwhite, 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. 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.
of the external resistance and the Seebeck coefficient with a fixed value for the internal resistance during summer for each color configuration. Figures 21a-d are for the color configurations blackwhite, 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. 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 Q from Equation (26) using the expression: The resulting optimal power generation equation is: 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 Q P-Optimal from Equation (26) using the expression: 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 Q from Equation (26) using the expression: The resulting optimal power generation equation is:  The resulting optimal power generation equation is: 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.

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.  Nomenclature A = surface area al = Albedo b = Mars eccentricity C = Sutherland s constant D = Drag force f(z, τ)= Normalized net flux F dt = View factor top-side to the sky F alt = View factor top-side to the ground F db = View factor bottom-side to the sky F alb = View factor bottom side to the ground G b = Direct beam irradiance on surface G bh = Direct beam irradiance horizontal surface G dh = Diffuse irradiance on horizontal surface G h = Global irradiance on horizontal surface G ob = Beam irradiance at top of the Atmosphere h = heat transfer coefficient I = Current K = Atmospheric thermal conductivity L = Chord length L s = Areocentric Longitude Nu = Nusselt number P = Atmospheric Pressure Pr = Prandtl number P h and P c = Hot and cold sources P Conductive = Conductive heat transfer P convective = Convective heat loss P radiative = Radiative heat loss P solar = Solar energy absorbed P TEG = Electrical Power R = Ideal gas constant R i = Internal resistance Re = Reynolds number R eff = Effective conductive resistance of the wing R L = External resistance T 0 = Reference temperature for Sutherland s T amb = Ambient temperature T grd = Ground temperature T s = Wing surface temperature T sky = 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) µ 0 = Reference viscosity for Sutherland s Formula µ = Boundary layer viscosity µ a = Atmospheric viscosity ρ = Boundary layer density ρ a = Atmospheric Density σ = Stefan-Boltzmann Constant τ = Optical depth ϕ = Latitude ω = Hour angle