# CO2 footprint minimization of solar-powered HALE using MDO and eco-material selection

Edouard Duriez, Víctor Manuel Guadaño Martín, Joseph Morlier

<figure><img src="image_1.png" /></figure>

## To cite this version:

Edouard Duriez, Víctor Manuel Guadaño Martín, Joseph Morlier. CO2 footprint minimization of solar-powered HALE using MDO and eco-material selection. Scientific Reports, 2023, 13 (1), pp.11994. 10.1038/s41598-023-39221-3. hal-04171232

HAL Id: hal-04171232

https://hal.science/hal-04171232v1

Submitted on 24 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# CO₂ footprint minimization of solar-powered HALE using MDO and eco-material selection

Edouard Duriez¹,³, Víctor Manuel Guadaño Martín²,³ & Joseph Morlier ID¹,³✉

Multidisciplinary Design Optimization (MDO) enables one to reach a better solution than by optimizing each discipline independently. In particular, the optimal structure of a drone varies depending on the selected material. The CO₂ footprint of a solar-powered High Altitude Long Endurance (HALE) drone is optimized here, where the structural materials used is one of the design variables. Optimization is performed using a modified version of OpenAeroStruct, a framework based on OpenMDAO. Our EcoHale framework is validated on a classical HALE testcase in the MDO community (FBhale) constructed using high-fidelity codes compared to our low-fidelity approach. The originality of our work is to include two specific disciplines (energy and environment) to adapt to a new problem of CO₂ minimization. The choice of eco-materials is performed in the global MDO loop from a choice of discrete materials. This is achieved through a variable relaxation, enabling the use of continuous optimization algorithms inspired from multimaterial topology optimization. Our results show that, in our specific case of electric drone, the optimal material in terms of CO₂ footprint is also the optimal material in terms of weight. It opens the door to new researches on digital microarchitected materials that will decrease the CO₂ footprint of the drone.

High Altitude Long Endurance (HALE) drones driven by solar power could be an alternative to satellites for some missions. Their solar power cells and batteries enable them to fly for a few years, which added to their high altitude (above 20 km) make them fit for missions similar to those of satellites¹. HALE drones can offer permanent coverage of a point or be re-positioned, and are repairable, unlike satellites, which are in orbit. Their lower altitude can offer better resolution for earth observation, but also results in smaller coverage. Their biggest advantage is their lower cost compared to satellites. HALE drones could also be more environmentally friendly as they do not require a high energy consuming launcher. This advantage can be enhanced if a special attention is paid to their environmental impact. For a fully electric HALE drone, most of this impact comes from the materials used and the manufacturing of the drone. Multidisciplinary Design Analysis and Optimization (MDAO) often compressed to MDO enables a better solution to be reached than by optimizing each discipline independently by developing efficient architecture². The most advanced numerical framework is developed by NASA³. MDO has been successfully applied to the ecodesign of commercial planes⁴. However, HALE drone mission requirements are very different from commercial planes. In addition to the aforementioned features, flight speed is not a requirement as they fly in closed trajectories. Moreover, a large wing surface area is required in order to gather sufficient solar power. These requirements lead to very different designs, with high aspect ratios. The design of HALE drones has already been studied extensively. Global configuration optimization has been made to an ultralight 3.2 m span solar powered drone using analytical methodology for the conceptual design of such an aircraft⁵. Another related paper focused on the development of a multidisciplinary tool for analysis, design, and optimization of HALE UAVs which act as an "atmospheric satellites" with an extreme aspect ratio concept (500 ft wingspan) using standard approaches, ranging from conceptual design and mission analysis, to preliminary aircraft design methods⁶. In Montagnier et al.⁷ the mass of the flexible wing was minimized using composite materials. The cruise speed versus lift coefficient diagram revealed an optimal solution with a payload of approximately 4 percent of the total mass of 817 kg for a 69 m wing span. From an experimental point of view: scale-sized prototypes have been manufactured, for mechanical testing⁸,⁹ or aeroelastic studies¹⁰,¹¹. More recently, a multi-fidelity design vehicle framework has been developed¹². Botero et al.¹³ applied this framework to a Solar UAV design, highlighting the fact that the energy required for the larger vehicle is approximately twice that of

¹ICA, Université de Toulouse, ISAE-SUPAERO, MINES ALBI, UPS, INSA, CNRS, 3 Rue Caroline Aigle, 31400 Toulouse, France. ²Université de Toulouse, ISAE-SUPAERO, 10 avenue edouard Belin, 31400 Toulouse, France. ³These authors contributed equally: Edouard Duriez, Víctor Manuel Guadaño Martín and Joseph Morlier. ✉email: joseph.morlier@isae-supaero.fr
the smaller one (same 10 m wing span but different payload). The vehicles are not optimized for an objective. The optimizer is used to find a feasible vehicle using the mostly energy related constraints. Optimizations seeking higher fidelity were carried out¹⁴, and a design framework was built¹⁵. However, all these studies focused solely on mass optimization and did not consider the environmental footprint of the drones. The aim of this work is to close this gap in a global, low-fidelity but fast, environmental impact optimization using open source tools in order to produce reproducible research.

This work is therefore based on OpenAeroStruct (OAS), a global low-fidelity aerostructural optimization framework<sup>16,17</sup>. This framework is itself based on OpenMDAO<sup>3,18</sup>. The Fig. 1 presents the MDA/MDO process of a typical coupled aerodynamic and structural design problem (OAS), while the next chapter highlights the main differences with the original OAS tool and thus guide the reader to understand the HALE design process starting from a standard aircraft design process.

For a conventional aircraft, most of the CO₂ emissions come from the fuel burned. For a solar powered HALE drone, however, no fuel is burned and most of the CO₂ emissions come from the materials utilized and the manufacturing of the drone. If the structural material used is fixed, optimizing the environmental impact of the HALE drone is equivalent to optimizing the mass of the drone. Therefore, we decided to include the choice of structural material in the optimization. A recurring problem in optimal material selection is that the material handled as a discrete variable. Another problem is that optimal properties are usually antagonistic. For example a material with a higher Young's modulus will also have a higher density.

In order to address these difficulties, Ashby indexes can be used¹⁹. These allow materials to be ordered according to a unique index. However, this method only works when choosing a material for a given simple part under a given loading. This method has been extended for simultaneous material selection and geometry design²⁰ and more recently solved using Variational Autoencoder²¹. A new approach has also been developed by the authors to solve material and process selection and thus link ecodesign and topology optimization²²,²³. This method, however, is limited to a structure and cannot take into account the impact of this structure on a bigger system, along with the resulting feedback loops. Making the material variable continuous is a suitable solution for the gradient-based algorithm used in this paper and was first proposed to find the optimal material²⁴ and developed in the form of an interpolation in topology optimization²⁵. However, the use of this approach in a multidisciplinary design optimization framework is unique in the state of the art and thus constitutes a major contribution of this paper.

The paper is organized as follows. We present in section “HALE drone MDO formulation” the HALE drone model and the optimization framework. Then, we introduce in section “Results and discussion on CO2 footprint minimization” the numerical results and also give some hints on how to interpret them. Finally, we give some concluding remarks in section “Conclusion”.

## HALE drone MDO formulation

**OpenAeroStruct to Eco-HALE.** We choose to derive our work from the OpenAeroStruct package with the wingbox full design (¹⁷). OpenAeroStruct is mainly utilized to study commercial aircraft which we modified to model a HALE drone. To validate our code, we used the single-boom HALE drone test case (data from FBhale publication¹⁴).

First, as there is no need for a 2.5 G maneuver on a HALE drone, we changed the design points used. We used one cruise flight design point where the lift and power constraints must be satisfied, and one gust design point where the structural constraints must be satisfied.

The OpenAeroStruct constraint on performance, based on the Bréguet equation²⁶ cannot be used for a solar-powered HALE aircraft, as the total weight of the aircraft does not change during flight and the range is not limited. The Bréguet range is a simple yet mathematically sound relationship between the main physics (aero-dynamic, structural, and propulsion) of the aircraft. Instead, it is a power equilibrium that gives a performance constraint. Actually, the power used by the propulsion and the payload must be produced by the solar panels or stored in the batteries. The power needed for propulsion $P_{prop}$ is derived from the 1D equilibrium between thrust

<figure><img src="image_4.png" /><figcaption>**Figure 1.** An overview of MDO applied to an aerostructural wing design problem¹⁶,¹⁷.</figcaption></figure>
(T) and drag (D). Thrust and drag are expressed in Eqs. (1) and (2), where W is the total weight of the drone, C_d is its drag coefficient, C_l is its lift coefficient, and v its speed.

$$D = W \cdot \frac{C_d}{C_l} \qquad (1)$$

$$P_{prop} = T \cdot v \qquad (2)$$

These give a link between the power needed for propulsion and the total weight of the drone (Eq. 3).

$$P_{prop} = W \cdot v \cdot \frac{C_d}{C_l} \qquad (3)$$

The relationship we actually use is Eq. (4), in order to account for the propulsion efficiency η and the power needed by the payload and the avionics $P_{payload}$.

$$P_{needed} = \frac{W \cdot v \cdot C_d}{C_l \cdot \eta} + P_{payload} \qquad (4)$$

During the day, power is supplied by the solar cells. Therefore a minimum wing surface area ($S_{wing}$) is needed in order to place a sufficient number of solar cells. This is expressed in Ineq. (5).

$$S_{wing} > \frac{P_{needed}}{A_{PV}} \qquad (5)$$

$A_{PV}$ is the power produced per unit of area of solar cells. Ineq. (5) is the constraint on performance we use. $A_{PV}$ is chosen carefully in order to be in the same case as FBhale (¹⁴). Figure 2 comes from¹⁴, it was used to extract the mean solar power harvested during the 24h power equilibrium considered. This mean power was divided by the area of solar cells in¹⁴, and leads to $A_{PV} = 54W/m^2$.

The computed surface area of solar cells necessary to produce the needed power is taken into account in the total weight computation, meaning that weight and power are coupled, as can be seen in Fig. 3.

The mass of the solar cells added to the total mass is computed using a surface density from¹⁴.

The mass of the batteries is also taken into account. In order to compute this mass, an assumption that the batteries power the drone during a 13 h night period ($t_{night}$) is made. This assumption is also based on Fig. 2, and enables a great simplification of the problem by not modeling the flight and gliding phases. This time multiplied by the power needed gives the size of the battery in kWh, which in turn gives the mass of the batteries ($M_{bat}$) when divided by the energy density ($d_{bat}$ in kWh/kg). The energy density was taken equal to that of FBhale. The mass of the batteries is therefore computed with Eq. (6). Taking the battery mass into account adds to the weight and power coupling, as on Fig. 3.

$$M_{bat} = \frac{P_{needed} \cdot t_{night}}{d_{bat}} \qquad (6)$$

These battery and solar panel masses are distributed along each of the beam elements of the wing, as was done in OpenAeroStruct for the fuel mass. This is closer to reality than a concentrated mass, and relieves the wing, allowing thinner structural skins. Other masses are also taken into account but not distributed, such as the avionics, payload, and the propulsion system mass considered as fixed during optimization.

<figure><img src="image_5.png" /><figcaption>Figure 2. Solar and battery power for FBhale (¹⁴). The nighttime we considered starts at the end of the gliding phase.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Power and weight coupling: if the total weight changes, the power needed to thrust that weight changes; if the needed power changes, the battery weight, the photo-voltaic cells (PV) weight and possibly the wing weight, change, resulting in a total weight change. This coupling is solved by the same Gauss-Seidel solver as the aerostructural coupling.</figcaption></figure>

The propulsion mass ($M_{prop}$) models both the mass of the engine and the propeller. It is obtained as the product of $P_{prop}$ (Eq. 3) and the propulsion density ($d_{prop}$ in kg/W), as in Eq. (7).

$$M_{prop} = P_{prop} \cdot d_{prop} \qquad (7)$$

In order for the data to be as close as possible to FBhale to validate our framework, without modelling the propeller, this propulsion density is estimated based on¹⁴. This estimate is obtained by dividing the propulsion mass (extracted from Supplementary materials) by the power used for propulsion (extracted from Fig. 2). These motors and propellers are added as point masses onto the wing structure. As in Eq. (8), the mass of each motor ($M_{mot}$) is obtained by dividing the propulsion mass by the number of motors ($n_{mot}$).

$$M_{mot} = \frac{M_{prop}}{n_{mot}} \qquad (8)$$

The power system (solar panels, maximum power point tracker and batteries) is sized for a worst case corresponding to the winter solstice where the least solar power is available. By choosing a more favorable launch day, much more power is available during climb. Therefore, the climb phase does not influence the sizing of the power system. As a result, computing the propulsion mass as in Eq. (7) also avoids the need to model the climb phase of the drone. This phase is taken into account through the propulsion mass data from¹⁴. A shear gust wall is added to the framework, in order to compute the loads used for the sizing of the wings. It has the same magnitude as in¹⁴. Our model only takes into account the weight of the wing, batteries, solar panels, propulsion, maximum power point tracker (MPPT), avionics (only the fixed mass) and payload mass. In order to account for the extra mass of the harness, the landing gear, the interfaces, the horizontal and vertical tail, the boom and the pod, the total weight of the drone is increased by 10%. This value comes from Supplementary materials (¹⁴).

Finally, the operational conditions of the HALE, such as the high altitude and the low speed, imply a low Reynolds number that lies in the range of 150,000–200,000. For this reason, high values of lift coefficient ($C_l$) are needed and a low-drag configuration is mandatory in order to minimize the power needed for propulsion (Eq. 3). Thus, as proposed in²⁷, profile NACA 63412 is selected, which is a laminar airfoil with a high maximum lift coefficient ($C_{lmax}$), a low value of the moment coefficient ($C_m$) and a “drag bucket” covering the desired $C_l$ range (see Fig. 4).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Polar curve for NACA 63412 airfoil.</figcaption></figure>

**MDA framework summary.** The multidisciplinary design analysis (MDA) framework is summarized in Table 1. The method modeling each discipline is shown in this table, along with the implementation used and a reference.

MDA converges using the Gauss-Seidel method (also known as fixed-point iterations). Each discipline analysis computes its own set of coupling variables that is passed to other disciplines (analysis). At the end of the MDA process, each discipline returns the final set of coupling variables computed at convergence. Additionally, the chosen optimizer will solve the constrained optimization problem with all the interactions shown in Figure 5.

<table><caption>**Table 1.** Summary of multidisciplinary design analysis framework.</caption><thead><tr><td>Discipline</td><td>Method</td><td>Implementation</td><td>References</td></tr></thead><tbody><tr><td>Aerodynamics</td><td>VLM</td><td>OAS</td><td>28</td></tr><tr><td>Structure</td><td>Wingbox beams</td><td>OAS</td><td>17</td></tr><tr><td>Energy</td><td>Simple in-house method</td><td>Section “OpenAeroStruct to Eco-HALE”</td><td>Data from<sup>14</sup></td></tr><tr><td>Environmental</td><td>Proportional to mass</td><td>Section “MDO framework summary”</td><td>Data from<sup>29,30</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Component layout of the performance group adapted from OpenAeroStruct<sup>16,17</sup>. The coupling between power and weight is visible as a “PV\_surface” variable from which the weight of the batteries and solar panels are derived.</figcaption></figure>
The energy and environmental models are described in other sections of this paper, listed in Table 1. We briefly present hereafter the aerodynamics and structural models used. These are kept unchanged from the OpenAeroStruct¹⁶ implementation, apart from the changes to the structure model described in section “MDO framework summary”.

The aerodynamics model is the implementation of the VLM model²⁸ by OpenAeroStruct¹⁶. Each panel of the lifting surface is modeled as a horseshoe vortex. This horseshoe vortex is made of three vortex segments : one span-wise bound vortex located at 1/4 th of the panel's chord-wise length and two semi-infinite chord-wise vortices extending towards infinity in the free-stream direction. The total lifting surface is modeled as the superimposition of all these horseshoe vortices. Flow tangency conditions are imposed preventing normal flow through the panel. The aerodynamic forces on each panel are derived from the circulation strength of the vortices. Finally, lift and drag are derived from these forces. Skin friction drag is also estimated. The reader is referred to¹⁶ for further details concerning this implementation.

The structure model is the beam FEM with wing-box cross-sections implemented in OAS (OpenAeroStruct¹⁷). The span-wise discretization of this model (7 elements) is the same as for the VLM mesh. Each span-wise section of wing is modeled as a beam element. The cross section of this beam element is a hollow rectangle modeling the upper and lower skins of the wing and the two spars. This model is used to compute the worst case von Mises stress values that are used as a mechanical failure constraint. The reader is referred to¹⁷ for more detail on this implementation.

**MDO framework summary.** Adding more physical constraints in OAS. The aim of this part is to add a more detailed sizing process to the existing OAS code and to describe the objective and constraints of the HALE design problem. One has to notice that a mechanical failure constraint leads to insufficient skin thicknesses. Therefore, a new (linear) buckling constraint is added in the design process. This buckling constraint imposes that the stress in the top skin of each element stays below the buckling critical stress of that panel. In order to have only one constraint, a Kreisselmeier-Steinhauser (KS) aggregation is used (³¹,³²), as for the initial mechanical failure constraint. The buckling critical stresses are estimated by the case of a homogeneous rectangular curved plate subject to combined axial compression and shear. Therefore, two different buckling critical stresses are considered: an axial critical stress ($\sigma_c$) and a shear critical stress ($\tau_c$), which can be expressed in the form of Eqs. (9) and (10), respectively, where $k_c$ is the buckling coefficient for axial load, $k_s$ the buckling coefficient for shear, b the width of the panel, and D, the flexural stiffness. In our approximate buckling case, we take b to be the distance between the two spars (set at half the chord) and the values of $k_c$ and $k_s$ can be obtained graphically from the charts presented in³³, using the geometry of the panel as an input.

$$ \sigma_c = k_c \cdot \frac{\pi^2 \cdot D}{b^2} \qquad (9) $$

$$ \tau_c = k_s \cdot \frac{\pi^2 \cdot D}{b^2} \qquad (10) $$

Thus, the buckling constraint is approximated by a parabolic interaction, as in Ineq. (11), proposed in³³, $R_s$ and $R_c$ being the stress ratios for shear and axial compression, respectively. These stress ratios are defined as the ratio of the stress at buckling under combined loading to the buckling stress under simple loading.

$$ R_s + R_c < 1 \qquad (11) $$

We also add a constraint imposing that the skin thickness cannot be greater than half the wing thickness, in order to prevent the top and bottom skin from intersecting at the wing tips, where the chord is smaller and the thickness-to-chord ratio of the wing is also small. To simplify the MDO process, the lower bound is fixed at 1 mm for all materials such as homogeneous metals or homogenized composites materials. At the conceptual detailed design level this limitation can be easily treated.

**New objective function: CO2 footprint.** We choose the CO₂ emitted by the HALE drone during its life cycle (CO₂tot) as our objective function. As the drone does not burn any fuel during flight, it is assumed that the CO₂ is mainly emitted before the use of the drone and is in particular due to the materials used and their processing. As the aim is to have a simple and fast tool, we consider only the CO₂ emitted by the material used for the structure (CO₂struct), the CO₂ emitted by the solar panels (CO₂PV), and the CO₂ emitted by the batteries (CO₂bat), as these have the most influence on the total emitted CO₂. This is shown in Eq. (12).

$$ CO2_{tot} = CO2_{struct} + CO2_{PV} + CO2_{bat} \qquad (12) $$

The CO₂ emitted by the structure is computed as the product of the mass of the spars ($M_{spar}$) by the CO₂ footprint of the material used for the spars (CO₂mat1) and the product of the mass of the skins ($M_{skin}$) by the CO₂ footprint of the material used for the skins (CO₂mat2), as shown in Eq. (13). This CO₂ footprint of the materials is considered to be equal to a weighted sum of the CO₂ footprint of their primary production (CO₂p) and the CO₂ footprint of their recycling (CO₂r), depending on their recycled fraction in current supplies ($\eta_r$), as in Eq. (14).

$$ CO2_{struct} = M_{spar} \cdot CO2_{mat1} + M_{skin} \cdot CO2_{mat2} \qquad (13) $$
$$CO2_{mat} = \eta_r \cdot CO2_r + (1 - \eta_r) \cdot CO2_p \quad (14)$$

It can be seen in Eqs. (13) and (14), that the CO₂ footprint of the processing of the materials is not taken into account. We made this choice because CO₂ footprint of processing varies a lot depending on the process, and we have no information related to this stage. This choice is acceptable as the footprints of the most probable processes are much lower than the footprint of the material primary production, and this is especially true for the materials that were found to be optimal.

The CO₂ emitted by the solar panels is computed as the product of the power needed on board the HALE ($P_{needed}$ detailed in section “OpenAeroStruct to Eco-HALE”) and the emissions per power ($CO2/W$), as in Eq. (15). The value of $CO2/W$ is taken from²⁹.

$$CO2_{PV} = P_{needed} \cdot CO2/W \quad (15)$$

The CO₂ emitted by the batteries is computed as the product of the size of the battery ($P_{needed} \cdot t_{night}$ as in section “OpenAeroStruct to Eco-HALE”) and the emissions per Wh ($CO2/Wh$), as in Eq. (16). The batteries are assumed to be similar to electric car Li-ion batteries, and the value of the emissions per Wh is taken equal to the mean of the three batteries studied in³⁰.

$$CO2_{bat} = P_{needed} \cdot t_{night} \cdot CO2/Wh \quad (16)$$

**Design variables.** We used eight geometric design variables and a material design variable which we discuss in detail in section “MDO framework summary”. B-splines parameterized using four control points each were adopted in order to vary wing geometric design variables along the semi-span. The only angle design variable we used was the vector of control points for the twist, because it also enabled us to control the angle of attack. Therefore, our variable twist represents the sum of geometric twist and angle of attack at cruise conditions.

We kept the wingbox design variables of OpenAeroStruct, which are the control points for the skin thickness and the spar thickness. These are mainly determined by the buckling constraint and the mechanical failure constraint, respectively. We used the same thickness distribution for both the upper and lower skins, and the same thickness distribution for both the forward and rear spars.

In order to satisfy both the power constraint and the weight-equals-lift constraint, we kept the control points for wing geometry as thickness-to-chord ratio, span, chord and taper ratio as design variables. These variables can generate the particular geometries of HALE drones with very high aspect ratios.

Additionally, the motor spanwise location was added as a new design variable. Since we are considering a symmetrical twin-motor HALE, this new design variable is defined as the ratio of the distance between the plane of symmetry and the motor to the semi-span of the wing.

Some limitations must be imposed to the design space in order to compensate for some missing physics in the model. Wings cannot be too tapered because tip stall can happen otherwise. As stall is not modelled in our framework, we choose a lower limit of 0.3 for the taper ratio design variable. For the same reason, a lower and upper twist limits of -15 and 15 degrees are considered. In this way, the critical or stalling angle of attack is taken into account. Furthermore, wings cannot be too narrow because 1-cosine gusts limit their aspect ratio. These gusts are not modelled, but the rule exists in the FAA (Federal Aviation Administration) regulation, so we also added a lower limit of 1.4 m for the root chord.

**Material selection.** The battery and solar panel types are fixed in our case. Therefore, minimizing the CO₂ emissions due to the solar panels and batteries, consequently, minimizes the power needed, which means minimizing the mass of the drone to be thrusted. Similarly, minimizing the CO₂ emissions due to the structural materials means minimizing the mass of materials used, if the choice of materials is not a variable. Therefore, if the structural materials were fixed, optimizing the CO₂ impact of the drone would be equivalent to optimizing the mass of the drone. Therefore, we chose to add the choice of the structural materials as a design variable. This new additional design variable takes the form of a vector with two components: the density of the material used for the spars and the density of the material used for the skins. Distinguishing these two structural elements adds new interesting degrees of freedom to our problem.

One of the main simplifications of our approach, is the material selection through their density. For the sake of simplicity this choice is also linked to the MDO process that needs to select a unique material at the end. The material data we need (Young's modulus, shear modulus, yield strength, and CO₂ emissions) is therefore accessed as a function of the density design variable. This limitation was solved in one of our recent publications²².

In order to have only continuous variables, we decide to make the density design variable continuous. This is achieved by linearly interpolating each material property in the space between real materials. The resulting interpolation is shown on Fig. 6 for two isotropic homogenized materials: Carbon Fiber Reinforced Composite (CFRP), Glass Fiber Reinforced Composite (GFRP) and two homogeneous metals: aluminium and steel. This is one of the limits of our approach, always linked to the simplification of the MDO process. Of course, this will lead to sub-optimal composite materials compared to more detailed structural sizing process.

However, with this interpolation, it is possible to find an optimum in an interpolated material corresponding to no real material. In order to eliminate this possibility, the interpolation is penalized by the use of a power term in the interpolation, as in²⁵. An example of this interpolation for Young's modulus is given in Eq. (17), where E is the interpolated Young's modulus at density $ρ$, $ρ_i$ and $ρ_{i+1}$ are the densities of the real materials framing the one being interpolated, $E_i$ and $E_{i+1}$ are the respective Young's modulus of these real materials, and p is the power used for penalization.
<figure><img src="image_9.png" /><figcaption>Figure 6. Young’s modulus example of linear interpolation of materials.</figcaption></figure>

$$E(\rho) = A \cdot \rho^p + B \qquad (17)$$

$$\text{with } A = \frac{E_{i+1}-E_i}{\rho_{i+1}^p-\rho_i^p} \text{ and } B = E_i - A \cdot \rho_i^p$$

The resulting interpolation is shown on Fig. 7 for $p = 5$, and for the same real materials as on Fig. 6.

To the left of a real material, a small decrease in density brings for example a great decrease in Young’s modulus. Therefore, the slight advantage due to a slightly lower density is compensated by the Young’s modulus evolution. To the right of a real material, even a big increase in density brings a very small increase in Young’s modulus. Therefore, the slight advantage due to a slightly higher Young’s modulus is compensated by the higher density. As a result, choosing a material that is not real (non-physical) is not beneficial, and the final optimum will be a real material.

In order to not advantage these non-physical materials, the penalization isn’t used in the case where a lower density leads to a better property. This appears for example in the left-hand part of Fig. 7, between CFRP and GFRP, where a decrease in density leads to an increase in Young’s modulus. For material properties for which a smaller value is better, such as CO₂ emissions, the penalization is inverted, as in Fig. 8.

This penalization must be applied to the interpolation only at the end of the optimization. In fact, as we use a gradient optimization, once the penalization is given to the interpolation, the solution is stuck in the local minimum corresponding to the closest material. The optimal material being an easily identified angular point (for example in Fig. 8), it is often also the optimum found even without any penalization in the interpolation. However, adding the penalization is more robust. Table 2 summarizes the MDO framework.

It is important to highlight that the SLSQP (Sequential Least Squares Programming) optimization algorithm is used in this work. The chosen optimizer options (stopping criteria) are $10^{-3}$ for the convergence accuracy and 250 for the maximum number of iterations. We checked the analytical derivatives using the complex step method³⁴.

<figure><img src="image_10.png" /><figcaption>Figure 7. Young’s modulus example of penalized interpolation of materials.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** CO₂ emissions example of penalized interpolation.</figcaption></figure>

<table><caption>**Table 2.** Summary of multidisciplinary design optimization framework.</caption><thead><tr><th>Objective function</th><th>Dimension</th><th>Bounds</th></tr></thead><tbody><tr><td>CO<sub>2tot</sub></td><td>ℝ</td><td></td></tr><tr><td>Design variables</td><td></td><td></td></tr><tr><td>Density</td><td>ℝ<sup>2</sup></td><td>[400, 8000] kg/m<sup>3</sup></td></tr><tr><td>Twist control points</td><td>ℝ<sup>4</sup></td><td>[-15, 15] deg</td></tr><tr><td>Skin thickness (t<sub>skin</sub>) control points</td><td>ℝ<sup>4</sup></td><td>[0.001, 0.1] m</td></tr><tr><td>Spar thickness control points</td><td>ℝ<sup>4</sup></td><td>[0.001, 0.1] m</td></tr><tr><td>Thickness-to-chord ratio control points</td><td>ℝ<sup>4</sup></td><td>[0.01, 0.4]</td></tr><tr><td>Span</td><td>ℝ</td><td>[1, 1000] m</td></tr><tr><td>Root chord</td><td>ℝ</td><td>[1.4, 500] m</td></tr><tr><td>Taper ratio</td><td>ℝ</td><td>[0.3, 0.99]</td></tr><tr><td>Motor location over semi-span ratio</td><td>ℝ</td><td>[0, 1]</td></tr><tr><td>Constraints</td><td></td><td></td></tr><tr><td>Mechanical failure σ &lt; σ<sub>max</sub></td><td>ℝ<sup>7</sup></td><td></td></tr><tr><td>Buckling R<sub>s</sub><sup>2</sup> + R<sub>c</sub> &lt; 1</td><td>ℝ<sup>7</sup></td><td></td></tr><tr><td>Skin thickness 2t<sub>skin</sub> &lt; t<sub>wing</sub></td><td>ℝ<sup>4</sup></td><td></td></tr><tr><td>Power equilibrium P<sub>needed</sub>/A<sub>PV</sub> &lt; S<sub>wing</sub></td><td>ℝ</td><td></td></tr></tbody></table>

## Results and discussion on CO₂ footprint minimization

Once the framework was validated (see Supplementary materials), an optimization was run. It should be recalled that our goal was to optimize the CO₂ footprint of the HALE, therefore, the objective function is now this total CO₂ emitted. Nevertheless, gradient based solvers converge to a locally optimal point, so the search for a global optimum depends heavily on the starting point of the optimization. For this reason, the multi-start strategy described in Table 3 was used.

Every combination of values is used, meaning that a total of 2 · 3 · 3 · 3 · 4 = 216 optimizations are run. For each variable, the values are regularly spread between the highest and lowest. For example, for the skin thickness control points, the three values given at the start to that variable are [0.001 0.002 0.003 0.004], [0.0015 0.003 0.0045 0.006], and [0.002 0.004 0.006 0.008].

The optimization was run with seven homogeneous materials integrated into the material design variable. Their properties are shown in Table 4.

Those named material 1, 2 and 3 are roughly homogenized sandwiched panels. Their skins are of the same carbon fibre reinforced polymer (CFRP) as the CFRP material in Table 4. The expanded polystyrene (PS) foam core of material 1 is replaced by a balsa core in material 2 and a cork core in material 3. Those cores are heavier than the one of material 1 but emit less CO₂ per kilogram of material.

Two material indices are computed according to the method of¹⁹. Both indices correspond to minimizing the CO₂ emitted by the structure CO₂struct, (defined in Eq. 13), for a flat plate loaded with in-plane compression. The buckling index is useful in the case where the buckling constraint described in section “MDO framework summary” is active. In this case, the index to be maximized can simply be established in 18 as:

$$E^{1/3} / \rho / \text{CO}_2\text{mat} \qquad (18)$$
<table><caption>**Table 3.** Design variable starting values for multi-start strategy.</caption><thead><tr><td>Design variable</td><td>Lowest starting value</td><td>Highest starting value</td><td>Number of starting values</td><td>Unit</td></tr></thead><tbody><tr><td>Density</td><td>[500 500]</td><td>[600 600]</td><td>2</td><td>kg/m³</td></tr><tr><td>Twist control points</td><td>[10 15 15 15]</td><td>-</td><td>1</td><td>deg</td></tr><tr><td>Skin thickness control points</td><td>0.002·[0.5 1 1.5 2]</td><td>0.004·[0.5 1 1.5 2]</td><td>3</td><td>m</td></tr><tr><td>Spar thickness control points</td><td>0.001·[1 1 1 1]</td><td>0.003·[1 1 1 1]</td><td>3</td><td>m</td></tr><tr><td>Thickness-to-chord ratio control points</td><td>0.05·[0.75 1 1 1.25]</td><td>0.17·[0.75 1 1 1.25]</td><td>3</td><td>-</td></tr><tr><td>Span</td><td>25</td><td>100</td><td>4</td><td>m</td></tr><tr><td>Root chord</td><td>1.5</td><td>-</td><td>1</td><td>m</td></tr><tr><td>Taper ratio</td><td>0.3</td><td>-</td><td>1</td><td>-</td></tr><tr><td>Motor location over semi-span ratio</td><td>0.3</td><td>-</td><td>1</td><td>-</td></tr></tbody></table>

<table><caption>**Table 4.** Material properties.</caption><thead><tr><td>Property</td><td>Material 1</td><td>Material 2</td><td>Material 3</td><td>CFRP</td><td>GFRP</td><td>Aluminum</td><td>Steel</td><td>Unit</td></tr></thead><tbody><tr><td>Density</td><td>504.5</td><td>529</td><td>560.5</td><td>1565</td><td>1860</td><td>2800</td><td>7750</td><td>kg/m³</td></tr><tr><td>CO₂ emissions</td><td>44.9</td><td>42.8</td><td>40.3</td><td>48.1</td><td>6.18</td><td>8.66</td><td>3.28</td><td>kgCO₂/kg</td></tr><tr><td>Young's modulus</td><td>42.5</td><td>42.5</td><td>42.5</td><td>54.9</td><td>21.4</td><td>72.5</td><td>200</td><td>GPa</td></tr><tr><td>Shear modulus</td><td>16.3</td><td>16.3</td><td>16.3</td><td>21</td><td>8.14</td><td>27</td><td>78.5</td><td>GPa</td></tr><tr><td>Failure strength</td><td>587</td><td>237</td><td>587</td><td>670</td><td>255</td><td>445</td><td>562</td><td>MPa</td></tr><tr><td>Buckling index</td><td>0.1539</td><td>0.1543</td><td>0.1544</td><td>0.05049</td><td>0.24153</td><td>0.1720</td><td>0.2302</td><td>N¹·m³/kgCO₂</td></tr><tr><td>Strength index</td><td>25885</td><td>10484</td><td>25959</td><td>8901</td><td>22184</td><td>18331</td><td>22119</td><td>N·m·kgCO₂</td></tr></tbody></table>

where E is the Young's modulus of the material, $\rho$ its density, and $CO_{2mat}$, the $CO_2$ footprint of the material (defined in Eq. 14). The strength index is useful in the case where the mechanical failure constraint is active. In this case, the index to be maximized is defined in (19):

$$ \sigma_f / \rho / CO_{2mat} \qquad (19) $$

where $\sigma_f$ is the failure strength (yield strength or tensile strength depending on the material). These material indices appear in Table 4.

The $CO_2$ emissions and density convergence graphs obtained during optimization are shown on Figs. 9 and 10, respectively.

The convergence graphs demonstrate that the lightest material is selected for both spars and skins.

Additionally, the CAD model of the resulting HALE is shown in Fig. 11 (with generic boom and tail shapes since only the wing structure is optimized). The optimization converges towards Material 1 for both spars and skins. This material is the one with the lowest density among those accessible through the material design variable. However, it is not the one with the best $CO_2$ emission indices. Indeed, the density indices for Materials 2 and 3 are both higher than those of Material 1. The optimization was re-run with fixed-material (1, 2 or 3) and confirmed that Material 1 performs better. This material is also the optimal material in case the total weight of the

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Objective function: total $CO_2$ emitted by the drone.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Material density convergence graph.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Rendering of the CAD model (using CATIA) for the optimal HALE wing.</figcaption></figure>

drone is the objective function. This optimization clearly confirms that in our specific case of an electric HALE drone, minimizing the mass is equivalent to minimizing the CO2 footprint. As expected, a slight increase in the total weight of the drone leads to an increase in the weight of the battery and the solar panel in order to propel a heavier drone, and also to an increase in the weight of the wing structure that induces a more important lift to compensate. These weight increases contribute to a further increase in the overall weight of the drone. This is called the “snowball” effect.

This result can be explained by the fact that the CO₂ emissions taken into account in the objective function are not only due to the structure (Eq. 12) explains this result. Although the structure composed of Material 3 emits less CO₂ than the one with Material 1, it is also heavier. The structure being heavier with Material 3, it needs more power to be thrusted, and thus, more CO₂ is emitted by the batteries and solar panels. This compensates the lower emission due to the structure. As a result, material indices as in¹⁹ can't be used in this case. A direct ratio between the weight of the structure and the CO₂ emitted by the batteries and solar panels would be needed in order to have a truly useful material index. However, this ratio is not accessible as it evolves during optimization. Therefore, a method such as the gradient optimization proposed in this work is necessary to choose the best material, if all material candidates can't be tested individually.

An informal test was carried out to compare the relative influence of the density and the CO₂ emissions of a material on its being optimal. For this test, the CO₂ emissions of Material 3 are lowered until this material becomes optimal. In this test, Material 1 and Material 3 differ only by their densities and CO₂ emissions. The density of Material 3 is 11% higher than that of Material 1. The results are presented in Fig. 12. We can note that the CO₂ emissions of Material 3 must be approximately 20% lower than those of Material 1, in order to compensate for its 11% higher density. It is reminded that Material 3 is the eco-material surrogate for Material 1, the expanded PS foam core of Material 1 being replaced by a cork core in Material 3. Therefore, in this case, in order to lower the total CO₂ emissions of the drone, an eco-material surrogate must also be almost as good as the initial material in terms of weight.

The code (see section Data Availability) is not optimized in terms of computational time. This will be addressed in future work. However, the speed of the optimization method presented in this article does not depend on the number of materials in the material catalogue. It is therefore particularly well-suited to large catalogs in comparison to a brute force method consisting in optimizing the structure for every material in the catalog. We also wish to emphasize that the results of our low fidelity tools are in very good agreement with the medium fidelity tools used for the design of FBhale¹⁴ as demonstrated in the mass minimization problem applied to FBhale (see Supplementary materials).
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Optimal material depending on changes to the CO₂ emissions of Material 3. CO₂ ratio is the ratio of the CO₂ emissions of Material 3 to the CO₂ emissions of Material 1.</figcaption></figure>

## Conclusion

The authors have developed a low fidelity MDO framework for preliminary design of solar powered HALE. The code is opensource and available for the community. Material choice was successfully integrated in a continuous multidisciplinary design optimization regarding CO₂ emissions. The simple HALE model developed showed acceptable agreement with more complex models such as FBhale, which was designed using much more sophisticated and costly tools.

Our modified version of OpenAeroStruct includes certain improvements with respect to the original version of the tool. For instance, it is adapted to HALE drones and contains more physics: material choice from among a discrete catalogue, batteries, solar panels, buckling constraint, etc.

A key result is that an optimal material in terms of drone total weight is also optimal in terms of drone total CO₂ emissions, even for an electrical drone. In order to be competitive in terms of total CO₂ emissions, an eco-material substitute must be almost as good as the initial in terms of drone total weight. This is due to a “snowball” effect on weight. One solution for tackling this problem is to generate this eco-material substitutes by creating digital materials based on unit cell optimization and apply this multiscale topology scheme to ribs and spars²,³,⁵.

It is true that HALE are highly flexible aerostructures. Thus more extensive sizing such as aeroelastic analysis and gust response should be conducted in a near future.

Finally, a good feature that could be addressed in the future would be a multi-objective optimization between CO₂ footprint and cost.

## Data availability

The optimization presented in the results can be obtained by following the instruction of the file howToStart at: https://github.com/mid2SUPAERO/ecoHALE/tree/downloadEcohale. The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Received: 10 February 2023; Accepted: 21 July 2023

Published online: 25 July 2023

## References

1. Dinc, A. et al. A performance review of a high altitude long endurance drone. *Int. J. Progress. Sci. Technol.* **24**, 32–39. https://doi.org/10.52155/ijpsat.v24.1.2461 (2021).

2. Martins, J. R. & Lambe, A. B. Multidisciplinary design optimization: A survey of architectures. *AIAA J.* **51**, 2049–2075 (2013).

3. Gray, J. S., Hwang, J. T., Martins, J. R., Moore, K. T. & Naylor, B. A. Openmdao: An open-source framework for multidisciplinary design, analysis, and optimization. *Struct. Multidiscipl. Optim.* **59**, 1075–1104 (2019).

4. Antoine, N. & Kroo, I. Framework for aircraft conceptual design and environmental performance studies. *AIAA J.* **43**, 2100–2109. https://doi.org/10.2514/1.13017 (2005).

5. Noth, A. *Design of Solar Powered Airplanes for Continuous Flight*. Doctoral Thesis, ETH Zurich (2008). https://doi.org/10.3929/ethz-a-005745622.

6. Morrisey, B. & McDonald, R. Multidisciplinary design optimization of an extreme aspect ratio HALE UAV. In *9th AIAA Aviation Technology, Integration, and Operations Conference (ATIO)* (American Institute of Aeronautics and Astronautics, 2009). https://doi.org/10.2514/6.2009-6949.

7. Montagnier, O. & Bovet, L. Optimisation of a solar-powered high altitude long endurance uav. In *Proceedings of the 25th International Congress of the Aeronautical Sciences* (Nice, 2010).

8. Romeo, G., Frulla, G., Cestino, E. & Corsino, G. HELIPLAT: Design, aerodynamic, structural analysis of long-endurance solar-powered stratospheric platform. *J. Aircr.* **41**, 1505–1520. https://doi.org/10.2514/1.2723 (2004).

9. Frulla, G. & Cestino, E. Design, manufacturing and testing of a HALE-UAV structural demonstrator. *Compos. Struct.* **83**, 143–153. https://doi.org/10.1016/j.compstruct.2007.04.008 (2008).

10. Cesnik, C. E. S., Senatore, P. J., Su, W., Atkins, E. M. & Shearer, C. M. X-HALE: A very flexible unmanned aerial vehicle for nonlinear aeroelastic tests. *AIAA J.* **50**, 2820–2833. https://doi.org/10.2514/1.J051392 (2012).

11. Jones, J. R. & Cesnik, C. E. S. Preliminary flight test correlations of the X-HALE aeroelastic experiment. *Aeronaut. J.* **119**, 855–870. https://doi.org/10.1017/S0001924000010952 (2015).
12. MacDonald, T., Clarke, M., Botero, E. M., Vegh, J. M. & Alonso, J. J. SUAVE: An open-source environment enabling multi-fidelity vehicle optimization. In *18th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference* (American Institute of Aeronautics and Astronautics, 2017). https://doi.org/10.2514/6.2017-4437.

13. Botero, E. M. et al. Suave: An open-source environment for conceptual vehicle design and optimization. In *54th AIAA aerospace sciences meeting* 1275 (2016).

14. Colas, D., Roberts, N. H. & Suryakumar, V. S. HALE multidisciplinary design optimization part I: Solar-powered single and multiple-boom aircraft. In *2018 Aviation Technology, Integration, and Operations Conference* (American Institute of Aeronautics and Astronautics, 2018). https://doi.org/10.2514/6.2018-3028.

15. Hwang, H., Cha, J. & Ahn, J. Solar UAV design framework for a HALE flight. *Aircr. Eng. Aerospace Technol.* **91**, 927–937. https://doi.org/10.1108/AEAT-03-2017-0093 (2019).

16. Jasa, J. P., Hwang, J. T. & Martins, J. R. R. A. Open-source coupled aerostructural optimization using Python. *Struct. Multidiscip. Optim.* **57**, 1815–1827. https://doi.org/10.1007/s00158-018-1912-8 (2018).

17. Chauhan, S. S. & Martins, J. R. R. A. Low-fidelity aerostructural optimization of aircraft wings with a simplified Wingbox model using OpenAeroStruct. In *EngOpt 2018 Proceedings of the 6th International Conference on Engineering Optimization* (eds. Rodrigues, H. et al.) 418–431 (Springer International Publishing, 2019). https://doi.org/10.1007/978-3-319-97773-7_38.

18. Gray, J., Moore, K. & Naylor, B. OpenMDAO: An open source framework for multidisciplinary analysis and optimization. In *13th AIAA/ISSMO Multidisciplinary Analysis Optimization Conference* (American Institute of Aeronautics and Astronautics, 2010). https://doi.org/10.2514/6.2010-9101.

19. Ashby, M. F. *Materials Selection in Mechanical Design*, 2 edition (Elsevier Butterworth-Heinemann, 2004).

20. Rakshit, S. & Ananthasuresh, G. K. Simultaneous material selection and geometry design of statically determinate trusses using continuous optimization. *Struct. Multidiscipl. Optim.* **35**, 55–68. https://doi.org/10.1007/s00158-007-0116-4 (2007).

21. Chandrasekhar, A., Sridhara, S. & Suresh, K. Integrating material selection with design optimization via neural networks. *Eng. Comput.* **38**, 4715–4730 (2022).

22. Duriez, E., Azzaro-Pantel, C., Morlier, J. & Charlotte, M. A fast method of material, design and process eco-selection via topology optimization, for additive manufactured structures. *Clean. Environ. Syst.* **2023**, 100114 (2023).

23. Duriez, E., Morlier, J., Azzaro-Pantel, C. & Charlotte, M. Ecodesign with topology optimization. *Procedia CIRP* **109**, 454–459 (2022).

24. Bendsoe, M. P., Guedes, J. M., Haber, R. B., Pedersen, P. & Taylor, J. E. An analytical model to predict optimal material properties in the context of optimal structural design. *J. Appl. Mech.* **61**, 930–937. https://doi.org/10.1115/1.2901581 (1994).

25. Zuo, W. & Saitou, K. Multi-material topology optimization using ordered SIMP interpolation. *Struct. Multidiscipl. Optim.* **55**, 477–491. https://doi.org/10.1007/s00158-016-1513-3 (2017).

26. Cavcar, M. Bréguet range equation?. *J. Aircr.* **43**, 1542–1544 (2006).

27. Mattos, B. S. D., Secco, N. R. & Salles, E. F. Optimal design of a high-altitude solar-powered unmanned airplane. *J. Aerosp. Technol. Manage.* **5**, 349–361 (2013).

28. Anderson, J. *Fundamentals of Aerodynamics* (McGraw-Hill, 1991).

29. Wetzel, T. & Borchers, S. Update of energy payback time and greenhouse gas emission data for crystalline silicon photovoltaic modules. *Prog. Photovolt. Res. Appl.* **23**, 1429–1435. https://doi.org/10.1002/pip.2548 (2015).

30. Hao, H., Mu, Z., Jiang, S., Liu, Z. & Zhao, F. GHG emissions from the production of lithium-ion batteries for electric vehicles in China. *Sustainability* **9**, 504. https://doi.org/10.3390/su9040504 (2017).

31. Kreisselmeier, G. & Steinhauser, R. Systematic control design by optimizing a vector performance index. *IFAC Proc. Vol.* **12**, 113–117. https://doi.org/10.1016/S1474-6670(17)65584-8 (1979).

32. Lambe, A. B., Kennedy, G. J. & Martins, J. R. R. A. An evaluation of constraint aggregation strategies for wing box mass minimization. *Struct. Multidiscipl. Optim.* **55**, 257–277. https://doi.org/10.1007/s00158-016-1495-1 (2017).

33. Gerard, G. & Becker, H. Handbook of structural stability. Part 3. Buckling of curved plates and shells. In *Tech. Rep., NATIONAL AERONAUTICS AND SPACE ADMINISTRATION WASHINGTON DC* (1957).

34. Martins, J. R. R. A., Sturdza, P. & Alonso, J. J. The complex-step derivative approximation. *ACM Trans. Math. Softw.* **29**, 245–262. https://doi.org/10.1145/838250.838251 (2003).

35. Duriez, E., Morlier, J., Charlotte, M. & Azzaro-Pantel, C. A well connected, locally-oriented and efficient multi-scale topology optimization (emto) strategy. *Struct. Multidiscipl. Optim.* **64**, 3705–3728 (2021).

## Acknowledgements

The authors would like to thank Mauricio Dwek and Nicolas Martin from Granta Design for their support in using CES Selector. The authors would also like to thank Fondation ISAE-SUPAERO for funding and enabling the Research Assistant position that lead this article to fruition.

## Author contributions

Must include all authors, identified by initials, for example: J.M. conceived the methodology and the numerical experiment(s), V.M.G.M and E.D. conducted the numerical experiment(s) and analyzed the results. All authors reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-023-39221-3.

**Correspondence** and requests for materials should be addressed to J.M.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_16.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023