TOMARES

In aerospace industry, saving mass on spacecrafts always remain in large demand to save launch costs or increase the available payload mass. A case study is carried out designing a first concept of an additive manufactured flywheel of a reaction wheel, as it is one of the heaviest parts of wheel systems. As an objective the mass is minimized, while obtaining an angular momentum suitable according to mission requirements and maintaining recent performances. As references the SeaSAT mission and a commercial reaction wheel are used. The work includes a preliminary dimension of the flywheels design space by MATLAB calculations, where in total 15 shapes are analyzed and compared. The most promising design space is afterwards analyzed via the finite-element tool ANSYS and is defined as the reference flywheel. The reference flywheel is used for topology optimizations (ANSYS Topology Optimization), where different boundary conditions are considered. The final designed flywheel obtains 16% higher energy density than the reference flywheel and withstands the mission loads. It can be concluded that it was possible to design a flywheel obtaining less mass while keeping the expected performance.


Introduction
A common way to control the attitude of a satellite is the use of one or more wheel systems such as reaction wheels, which are part of the attitude and orbital control system (AOCS). One of the key element of a reaction wheel is the flywheel, which generates an angular momentum to adjust the satellite's orientation. The flywheel is next to the torque motor one of the heaviest parts of the attitude control component. From a satellite system's perspective each subsystem should have the lowest possible mass so that the payload mass could be maximized. Therefore heavy components are analysed to identify potential mass reduction.
Currently, almost all flywheels in use are constructed as solid or hollow cylinders [1]. Since the performance of a flywheel depends on its angular momentum, which is controlled by the moment of inertia and therefore mass distribution, a sufficiently large enough angular momentum can be created by an effective mass distribution; removing unnecessary mass and increasing useful mass. Topology optimization offers great potential to reduce the overall mass while obtaining desired structural properties and therefore to create bionic lightweight structures. Additive manufacturing (AM) enables to produce such filigree structures.
An analytical approach for the optimization of simple flywheel shapes using MATLAB has already been developed [2], as well as an approach based on SIMP, the solid isotropic method with penalization [3]. The development of a flywheel with the highest possible energy storage capacity was attempted by combining different materials and geometries in simulations by [4]. [5] followed the same goal by using parametric geometry modeling and a shape optimization method to optimize the flywheel rotor geometry. Utilizing topology optimization based on the variable density method, [6] also pursued the same goal. [7] investigated different flywheel 1 3 geometries and their energy storage capacity per unit mass. Additionally, the energy storage capacity was optimized by the Injection Island Genetic Algorithm [8] and by an evolutive system method [9]. Alternatively, [9] also considered unsymmetrical shapes. [10] and [11] followed a different approach. While [10] adapted the flywheel shape to achieve an even stress distribution by varying the thickness along the radius, [11] altered the flywheel thickness as well as the material properties (e.g. by using functionally graded materials) to optimize the flywheel concerning the stress to strength ratio. In summary, at present less work dealt with mass reduction of additive manufactured flywheels, which are used in satellites and significantly smaller, lighter and different in rotational speed than industrial flywheels.
Therefore a project has been initiated to further investigate on this matter. The project title is "Topology Optimization of an Additive Manufactured Reaction Flywheel Designed for an Earth-Observation Satellite" which can be abbreviated to the project name TOMARES. The project's objective is to create a first concept of an additive manufactured lightweight flywheel of a reaction wheel to decrease the overall mass of an AOCS subsystem by maintaining the demanded performances. For this project, a case study is performed using the SeaSAT mission and the Honeywell HR14 reaction wheel with an angular momentum of 25 Nms as reference. Due to the increased interest of PEEK for space applications [12] by reasons of its advantages in terms of lightweight, thermal expansion and elongation properties compared to traditionally used materials (Aluminium 7075, Titanium TA6V) [13], the high performance polymer PEEK (Polyether ether ketone) is used for this project. As magnetic bearings shall be applied in this project, PEEK additionally has the advantage that it obtains no magnetic permeability and therefore does not influence with e.g. the current of the motor.
The paper is structured as follows: Sect. 2 summarizes the relevant flywheel characteristics. In Sect. 3 different flywheel shapes that are currently in use are optimized according to the mission requirements as well as additional objectives. The optimized shapes are afterwards compared to find the most suitable overall shape of the flywheel, which is analysed according to on-board loads using the finite-element method (FEM) in Sect. 4. The best shape of Sect. 3 serves as basis for the topology optimization (Sect. 5). The characteristics of the final topology optimized model are afterwards analysed and compared to the findings of the preliminary flywheel in Sect. 6. The results are discussed in Sect. 7, while the conclusion is given in Sect. 8.

Flywheel characteristics
According to [7] the overall performance of a flywheel is dependent on (1) material strength, (2) flywheel geometry (flywheel cross-section) as well as (3) rotational speed. The material and geometry are altered hereinafter to find the optimal flywheel performance, while the rotational speed is specified by the set mission requirements. Following flywheel characteristics are relevant subsequently to identify the optimum flywheel.
Flywheels are currently designed to obtain a high stored kinetic energy E k , which is equal to the rotation energy. The kinetic energy of a flywheel depends on the rotational speed (respectively the angular velocity ) and the moment of inertia I about the rotation axis, which in turn depends on the mass m and the shape of the flywheel and consequently on the geometry of the cross-section. Mathematically E k is represented in Eq. 1 [14].
To show the effect of the flywheel geometry on its energy storage capability per mass, [7] defined the energy density e: Additionally, the mass specific angular momentum L mass is considered, putting the angular momentum L in relation to the flywheel mass.
During operation, the stress distribution within the flywheel is mainly driven by centrifugal forces, putting the material strength as a limiting factor. Thus, the Rankine criterion with its equivalent stress v provides statements regarding the durability of the flywheel, whereas the form factor K enables to compare different flywheels according to their internal equivalent stress distribution. If a uniform distribution is present, the form factor obtains a maximum value of one [15]. As insignificantly small shear stress is expected, the last term in the calculation of the Rankine criterion is neglected subsequently [16,17]. r : radial stress, t : tangential stress, : density.
According to [18], the Rankine criterion is chosen for two reasons. First, a high strength material will be used which show brittle fracture behaviour. Second, an equivalent stress is more significant than the maximum stress, since the present multiaxial stress state is reduced to an uniaxial equivalent stress. It can therefore be seen as more meaningful in relation to the yield stress of the material.

Outer contour definition
The overall dimensions of the reference flywheel Honeywell HR14, restricted by the dimensions of the motor, serves as basis for the design space to determine the overall flywheel shape. Resulting characteristics are summarized in Table 1.
Note that the HR14 achieves 25 Nms, which is also the target value for the new flywheel in this work. Since the development consists of two consecutive phases (optimizing outer contour and topology optimization of the outer contour), the angular momentum for the optimization of the outer contour must be greater than 25 Nms, due to the physical relationship between angular momentum and moment of inertia (which is proportional to the mass, see Eq. 3) and due to the flywheel mass reduction in the context of topology optimization. Otherwise the result of each topology optimization would be a flywheel with less than 25 Nms in angular momentum. For this reason an angular momentum of 50 Nms is requested for the optimization of the outer contour.
Basic shapes such as discs with and without central hole as well as with a constant, curved or conical profile are considered ( Fig. 1). Additionally, the shapes are extended either with a rim to generate a higher angular momentum or a shaft (and a hub) to achieve a better connection between the flywheel and the rest of the wheel system. The rim and shaft geometries are also considered in combination.
The aim of this work is to develop a flywheel that is as (1) light as possible, obtaining the set (2) angular momentum (see Table 1) as well as (3) not exceeding the maximum stress that is allowed due to the material properties. Therefore, 15 different shapes are analysed in MATLAB in several nested loops taking the above written objectives into account, whereas the analysis is done by changing variable sizes of the flywheel geometry.
The analysis is done as follows: first the material properties and maximum flywheel dimensions are imported enabling the calculation of the maximum possible flywheel mass m max , whereas all variable dimensions are set to its maximum values. Afterwards the actual analysis process starts by reducing the variable dimensions in sections of 1 mm until the minimum values are reached. For each dimensional combination of the sections, the system checks whether the mass is less than the original flywheel, the angular momentum sufficiently high compared to mission properties and the maximum stresses (radial and tangential stress, as well as the Rankine criterion) less than the yield stress acc . If all these conditions are fulfilled, the current dimension of the flywheel and the corresponding mass is stored and used as a comparison value for the next reduction iteration. In this way, when the minimum dimensions are reached, the best dimension can be obtained, for which further evaluation criteria and the flywheel profile are then calculated and exported. Considered as variable are the thickness and radius of rim and plate ( y rim and r rim , respectively y plate and r plate ). Contrary, the radius of the shaft r shaft has a constant value, as well as the thickness y hub and radius r hub of the hub, since both values depend on the radius of the shaft. The thickness of the shaft y shaft is partially variable, as it depends on the variable thickness of the plate (or hub, if existent) and the length needed to hold the bearings and the motor. During the analysis process each variable is changed independently until the requirements are fulfilled. A workflow-chart of the analysis done in MATLAB is shown in Fig. 2.
All geometries show a rotational symmetry around their z-axis and a horizontal plane of symmetry in the xy-plane. Except from the conical and curved plate, all parts of the Table 1 Requested flywheel characteristics, based on the Honeywell HR14 according to [19] a For the outer contour definition, b For the topology optimization and the final flywheel design  The hub is always placed between the shaft and plate and its thickness is always higher than the plate's, while the rim always obtains a value higher than the plate, but smaller than the shaft. If possible in regard to other constraints, every variable has a minimum thickness and radius of 1 mm. For a detailed description of the corresponding formulas and calculations of radii and heights of each flywheel section, it is referred to the work of [15][16][17][18]. For the sake of simplicity, only a few basic aspects of the calculation, as well as the assumed simplifications and assumptions are described below.
For the calculation of stress and deformation in flywheels with a constant angular velocity, a static load case can be assumed, since the wheel rotates with a constant angular velocity and centrifugal forces are applied as external force to the mass elements [18]. Further, bending, torsion and oscillation is neglected, as bending and torsion can be seen as minor compared to the appearing radial and tangential stress [16]. This is due to the fact that the wheel spins around the vertical, center-of-gravity crossing, rotational axis at constant speed, so that no moment is generated in axial direction [15]. In addition, an axis symmetrical plane stress condition is assumed, since no or insignificantly small shear stress appears in radial direction and the thickness is small with respect to the radius [15][16][17]. Moreover, the entire flywheel is homogeneous and made of the same material. Loads at the inner (if existing) and outer surface are assumed as equably distributed. The radial stress at the maximum radius is set to zero. Gravity is neglected and the temperature in the entire flywheel is assumed as constant.
All stress values are multiplied by a safety factor of s = 2 as common safety factors for conventional materials lie between 1.25 and 2 [20].
Based on the work of [15] and [17] the differential system of equation is used to calculate the stress of the flywheel, whereat the height is described as a function of the radius Fig. 2 Overview of the workflow chart of the analysis process for the pre-dimensioning of the flywheel for different flywheel geometries conducted in MATLAB. The left side of the figure shows the overall process and the right side a detailed view of the analysis process itself (visualized through the light red rectangle). During each iteration the variable values are reduced for every flywheel sections until the minimum values are reached and for each combination, the system checks whether the mass is less than the original flywheel, the angular momentum sufficiently high compared to mission requirements and the maximum stresses less than the yield stress. If a value or condition is in parentheses, the value or condition only counts for flywheel shapes that have these properties. Otherwise, it is neglected. Rectangles in light blue symbolizes processes, while hexagons stand for conditions. The orange conditions belong to the properties of the rim and the yellow ones to the plate properties. Green hexagons represent general conditions, which apply for all shapes. Circles visualize the start and end of the process and dark grey tetragons import and export procedures. m max : maximum possible mass of the flywheel; y: thickness of the corresponding flywheel section; r: radius of the corresponding flywheel section; r : radial stress; t : tangential stress; v : equivalent stress (Rankine criterion); acc : yield stress; K: form factor y = y(r) . Depending on the considered flywheel geometry, the differential equation system can be solved in different ways. For the sake of simplicity only the source for different solutions are mentioned. [16] describes the solution for a wheel which is made of one part with a constant height, similar to a flat disc. If the wheel is estimated by multiple parts [18] provides the required equations. Also [18] discussed the required equations if the wheel is made of one part with a central hole. If at least one part of the flywheel has a non-constant height, a conical or curved shape is chosen, [17] presents the equations for the conical and [16] the curved shapes.
To solve the equation system an approximation adapted from [17] by [15] is used, where the differential equations are converted to difference equations and the stress is calculated step by step from in-to outside of the wheel. The flywheel shape obtaining best results is improved by applying a correction curve according to [15] in order to further decrease stress peaks.
In summary, for the selection of a outer contour 15 flywheel shapes are considered and multiple variants for each of these 15 shapes are created based on a parameter analysis. By changing a variable parameter (see Table 2) by one millimeter per iteration, several hundred shapes per geometry are compared to each other to determine the shape with the lowest mass, the lowest stresses and an angular moment of 50Nms within a single shape. Afterwards, the best variants of every geometry are compared to each other to define the best flywheel shape.
Mathematically, the correlation can be described as y total : Thickness of the entire flywheel shape along all sections, r total : Radius of the entire flywheel shape along all sections.
Since the connection and influence of the individual parameters on the function is unknown, a quantity of possible parameter settings per shape are calculated and compared to define one parameter set that fits the objective the most.
Comparison parameters All flywheel shapes are compared according to mass, mass specific angular momentum (Eq. 3), maximum Rankine criterion (Eq. 5) and the form factor (Eq. 6). In addition, an average value of the Rankine criterion is calculated over the entire flywheel shape, so that not only stress peaks but the entire stress curve is considered in the evaluation.
Concerning the scoring, the points are awarded as follows: The mass, the maximum and the average Rankine criterion should be as small as possible while the mass specific angular momentum and the form factor should be as high as possible.
(7) f optimal (y total , r total , v,max , r,max , t,max , m, L) Results Comparing the 15 different shapes, the flywheel shapes combining plate, shaft, hub and rim provide best properties. Within this group, the "curved plate with shaft and rim" obtains best characteristics concerning the above mentioned parameters. Figure 3 and Table 3 show an extract of the evaluation.
The "curved plate with shaft and rim" has best results, which is why the correction curve is additionally applied on this shape. Its results are summarized in Fig. 4 and Table 4. As this flywheel shape obtains best results for this project, it is considered as the preliminary design for the following topology optimization.
Based on the results, several general statements can be made: (1) the curved and conical shape are similar and exhibit a lower and more constant stress flow, (2) every thickness difference between parts lead to a jump in the stress flow, (3) a hub decreases the stress flow jumps, (4) a curving of junctions lead to a decrease of the maximum stress, disappearance of stress jumps and improvement of all considered attributes and (5) a rim supports an optimal mass to angular momentum ratio.

FEM analysis of reference flywheel
The "curved plate with shaft and rim" containing the correction curve is used as a reference flywheel for the topology optimization and is tested using FEM (ANSYS). As boundary conditions, fixed displacements are placed on the upper and lower part of the shaft, while a rotational velocity is applied at the z-axis. Also static and dynamic imbalances and oscillations due to roll-yaw oscillations caused by reaction wheels themselves are taken into account. The former is modelled by applying a force at one side of the flywheel, while the latter is simulated by a pair of centrifugal forces as well as transverse moments that is applied in the xz-plane, as suggested by [21]. The responses of a reaction flywheel in the SeaSAT mission to roll-yaw oscillations was investigated by [22] by applying a white-noise disturbance. The same method is taken into account subsequently, creating the white noise in MATLAB and implementing it in the finite element (FE) model in the xy-and xz-plane for 100 s.
Due to missing data during launch, only on-board loads are considered. The corresponding data are taken from Honeywell HR14 and are summarized in Table 5. As the geometry is axially symmetric, only a quarter of the flywheel is evaluated, reducing the computational time significantly.
Static structural, modal and harmonic sine load cases are applied. For the static analysis, the impact of the rotational velocity is transferred to a radial acceleration load that is applied on the outer rim of the flywheel. Internal analyses show the accordance of the Rankine stress values for the application of rotational velocity and adequate radial acceleration loads. Only the maximum value appearing in the white noise is applied for the static structural analysis. Concerning the harmonic sine analysis, base-point excitation in three different directions (x, y, z), using an amplitude of 1 g, are performed and a damping coefficient of 2% is used as recommended by Calvi [23]. The used material properties for PEEK for the analysis are summarized in Table 6.
Results A mesh convergence study shows sufficient results for a mesh size of 4 mm using solid hexahedral elements. Figure 5 gives an overview of the static structural analysis. Due to the axially symmetric geometry, the results are exploited in a cylindrical coordinate system with the origin of coordinates in the middle of the flywheel. Highest stress values appear at the transition between the plate   Table 7.

Topology optimization of reference flywheel
Using ANSYS Topology Optimization (design constraint: mass -mass ratio: 0.4; objective function: minimize compliance) five different boundary conditions are tested to find the most suitable topology. The shaft as well as the outer rim of the rim is always taken as non-design space due to manufacturing constraints as well as to ensure a connection to the motor. Figure 6 summarizes the different design spaces, which are chosen to maximize the angular momentum (Non-Design 1 and 5), to provide minimal boundary conditions (Non-Design 2) and to maintain the transition (Non-Design 3) as well as the overall shape of the flywheel (Non-Design 4). The same FE model as described in the previous chapter is taken into account as basis for the topology optimization. Five optimizations are performed taking the five different design spaces into account. The outcome of each optimization is afterwards reconstructed in CAD (Computer-Aided Design) and static structural analyses including only rotational loads are performed, comparing the different shapes according to the highest stress values as well as mass and angular momentum. Static and dynamic imbalances are not considered within the optimization in order to obtain a symmetrical flywheel. As the raw topology optimization results obtain sharp edges and consequently high stress peaks appear within the static analysis, all five models that are obtained from the optimization are adapted to reduce stress peaks in several iterations by smoothing the model and implementing roundings.The models are adapted accordingly until for all optimized flywheels, stress values similar of the reference flywheel appear (max. aberration 50% of the highest stress values of reference flywheel). The models are afterwards compared according to the mass specific angular momentum. Only the flywheel model with highest mass specific angular momentum is considered further. An additional topology optimization is afterwards performed on the chosen model in order to obtain the correct angular momentum for the given mission requirements. The same optimization parameters are chosen (design constraint: mass, objective function: minimize compliance), while the mass ratio is decreased in 0.05 steps starting from 1.0. The mass is reduced until mission requirements concerning the angular momentum are fulfilled. Additionally, manufacturing constraints regarding selective laser sintering (SLS) are taken into account to prevent enclosed powder within the flywheel.
Results Figure 8 summarizes the different topology optimization results concerning the design space alterations. The different topologies are adapted to obtain similar stress values as present in the reference flywheel. The changes mainly include the smoothing of sharp edges, where highest stress peaks appear. Exemplary the results of the first and the last iteration with adapted topology are shown for Non-Design 1 in Fig. 7. Mass and angular momentum of the five optimized flywheels after the adaptations are given in Table 8.
The topologies of Non-Design 1 and 2 of the last generation obtain stress values most similar to the reference flywheel. As highest mass specific angular momentum appears for the Non-Design 1 model, this model is used for further optimizations.
The mass ratio of Non-Design 1 is further increased until mission requirements as presented in Table 1 are  fulfilled. An additional mass ratio of 0.65 obtains the required angular momentum. The structural changes of the topology of the second optimization compared to the above described flywheel can mainly be found in the connection between the upper and lower surface. While the previous Non-Design 1 obtains solid connections, the second topology optimization exchanges these by struts. As the flywheel is constructed and printed as a closed part, four holes are placed symmetrical on the outer side of the rim to allow powder release for the post-processing of the additive manufactured part while still keeping a symmetrical geometry. Figure 9 shows the final design. The properties are summarized in Table 9.

FEM analysis of optimized flywheel
The same FEM analyses are performed as described in Sect. 4, using an equivalent mesh size. Figure 10 shows the results of the static structural analysis. Additionally, the main properties of the analysis are summarized in Table 10.
Comparison Following main statements can be drawn between the reference flywheel and the final optimized flywheel: Static-structural analysis The final optimized flywheel obtains an stress increase of 74% compared to the reference flywheel. In addition, the areas where highest stress values appear shift. While in the reference flywheel the maximum stress values appear at the transition plate-rim, they are found in the final optimized model at the transition shaft-plate as well as the areas where the struts are placed. Modal analysis Higher 1st, but lower 2nd and 3rd eigenfrequency can be stated for the final optimized flywheel.
Harmonic sine analysis A higher response is observed in the optimized model, as the first mode shows 15% higher values, while at the second mode an increase of 21% is seen. Additionally, a shift from an excitation in x/y direction to a z-direction for the main mode appears.
The optimized flywheel shows higher values in all analysis. The mass however is reduced by 56% but the energy density (respectively the mass specific angular momentum) is increased by 16% compared to the reference model while maintaining mission properties. Table 11 gives an overview, inter alia, of the mechanical properties that are considered as important for this project.

Discussion
Non-Design 1 obtains the required properties of low mass while keeping mission requirements as a high angular momentum exists due to the high mass that is placed at the outer ring of the flywheel. The latter is less present for Non-Design 2 and 3. Inspired by the results of the second topology optimization the adaptation of the ring towards mission requirements was achieved in this project by reducing the ring material by implementing several holes in the outer ring of the flywheel until the corresponding angular momentum was reached. Alternatively, a decrease of the ring diameter for the adjustment leads to similar results. The design therefore shows high potential to easily be adjusted. Non-Design 2, 3 and 4 obtained high stress values as the flywheel start to flatten during operation resulting in high deformations. The flattening of Non-Design 1 is less present due to the two supporting connections in the plate and transition of plate-rim preventing high stress values. Hence, the highest stress values still appear at the transition as observed in the reference flywheel. Only small stress values are observed for the adapted Non-Design 1 as the areas with high stress (especially transition plate-rim) are supported with enough material. Non-Design 5 obtains a high angular momentum, however also high stress peaks appeared at the connection between the plate and rim and a rather low mass specific angular momentum due to the high mass in comparison to the other designs. It can therefore be concluded that Non-Design 1 shows most suitable properties for the set objective when comparing the different obtained topologies, as it obtains a high angular momentum, low mass and a good stress distribution, also due to the preliminary set overall shape of the flywheel.
By comparing the results of the static structural analysis higher stress values appear in comparison to the reference flywheel. Even though higher stress values were obtained for the optimized model, the maximum stress shows a value lower than 2 3 of the yield stress of PEEK and can therefore be seen as applicable. In addition, the main mode shifts to lower eigenfrequencies, but is still higher than the critical frequency range that is excited during launch. A total mass reduction of 56% compared to the set reference flywheel can be stated. However, as both flywheels obtain different properties concerning the angular momentum (see Sect. 3), the energy density e is taken into account, as the mass or angular momentum in particular does not allow a meaningful comparison. The values are summed up in Table 11. An increase of about 16% for the energy density can be stated, which means that the final optimized model can store more energy per mass than the reference flywheel.
It has to be noted that the actual shape of the Honeywell HR14 is representative for the "constant plate with shaft and rim", which was ranked second after the "curved plate with shaft and rim". Due to lower stress values, the latter was considered as having higher potential. Similar conclusions concerning the stress distribution and energy density were drawn by [6] and [9]. The here developed flywheel therefore enables to combine good properties concerning high and adjustable angular momentum, low stress and low mass.
As the best flywheel shape for this project should obtain a high angular momentum while keeping low mass, the optimization's objective should be chosen accordingly. However due to software limitations, this was not possible and the objective to obtain maximum stiffness for the design was selected as increased stiffness leads to higher eigenfrequencies which is an important property for lightweight structures in aerospace application. The design of the presented flywheel was therefore advanced to previous flywheels by the MATLAB parameter study as well as the topology optimization for various parameters to not only optimize for high angular momentum but rather also including its stress, stiffness and eigenfrequency properties. Nevertheless additional optimizations with a different software is recommended to further improve the design on this matter.

Conclusion
It can be concluded that it was possible to design a new flywheel topology that increases the energy density compared to nowadays flywheels by focusing on the reduction of the stress distribution within the flywheel by the preliminary design of the flywheel. The angular momentum and mass was afterwards adapted and decreased by topology optimizations enabling an increase of the energy density about 16% compared to the reference flywheel while fulfilling mission requirements and withstanding the set loads. Further investigations, including sensitivity analyses, further optimizations as well as including additional load cases should be done to further specify the design.
Availability of data and materials Further data can be provided if a request is sent to the corresponding author.

Compliance with ethical standards
Conflict of interest The authors declare that there is no conflict of interest.
Funding Open Access funding enabled and organized by Projekt DEAL. The author would like to acknowledge the support of ZARM Technik AG.
Code availability Code can be provided if a request is sent to the corresponding author.
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://creat iveco mmons .org/licen ses/by/4.0/.