Aerodynamic Design of Airfoil Shape for Gust Generation in a Transonic Wind Tunnel

This article presents the aerodynamic design of the airfoil of the gust generator system being developed in the GUDGET project and conceived to generate high-amplitude gusts in a transonic wind tunnel. The system is made of vanes creating a flow deviation in turn by flapping around a rotational axis or by blowing air though a suitable sonic jet located close to the vane trailing edge. The airfoil shape optimization has been carried out using a design of experiment technique (DOE) and response surface optimization along with URANS CFD. The computational model has been preliminarily validated using data provided by ONERA for the baseline design at a lower Mach number (M=0.73\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {M}=0.73$$\end{document}) and then compared with the one actually required by GUDGET in the test chamber (M=0.82\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {M}=0.82$$\end{document}). All the cases have been optimized at a frequency of 40 Hz and then investigated at a frequency of 80Hz.


Background and Motivation
Aerodynamic loads on an aircraft can significantly and suddenly change due to atmospheric disturbances in severe turbulence conditions. Among the main phenomena, there are natural instabilities of air masses, water-land discontinuities, thunderstorms. Those impact all phases of flight. Close to the ground, significant turbulence levels can be encountered due to the flow variations associated with the atmospheric boundary layer. Sustained velocities of 3-5 m/s (which can go up to 10 m/s) have been reported at high altitudes [1]. Notwithstanding these quite challenging conditions, they are used to size the airframe and compute the structure fatigue life. Discrete gust and continuous turbulence models to be used for loads evaluation and structural design of large airplanes are specified by airworthiness standards such as the EASA CS-25, they are actually the results of tens of decades of research and statistics on gusts' severity and related effects on flying aircraft [2].
The investigation of the aircraft responses to high-amplitude gusts has become one of the most important investigated topics in aeroelastic experiments and simulations in the last years. Due to the sudden occurring load changes during a gust event, the loads on different structural components, like wings or horizontal tail planes, can reach very high values, becoming very often the most critical load condition for the structural sizing. Therefore, an accurate prediction of such time dependent loads is very important for modern aircraft design to ensure structural integrity and to avoid unstable flight conditions [3].
To handle such strong load changes, experimental and numerical investigations must be carried out first for a 1 3 better understanding of the physical phenomena and latter to improve methods for an increased load prediction reliability. This may lead to the design of better active load control to reduce undesired vibrations and associated loads, by increasing passenger comfort and enhancing the fatigue life of load-carrying structures. The development of improved aeroelastic models requires anyway an extensive experimental database which is currently difficult to obtain, especially in the transonic regime and for high-amplitude gusts. One of the earliest and simplest generators was constructed by Hakkinen and Richardson [4] and involved a plunging plate in a wind tunnel. The apparatus was used in an attempt to experimentally validate the Sears function, but provided unsatisfactory performance due to sensor noise. Subsequent gust-generation systems, understood that a single actuated plate was not sufficient for producing gusts in a controlled way, became more and more complex. In [5], Bennett and Gilman used four plates, mounted in pairs on the walls of a wind tunnel and actuated together by a series of linkages, to produce sinusoidal gusts for experiments with scale models of aircraft. Ham et al. [6] and Jancauskas and Melbourne [7] generated gusts using a pair of controlled-circulation airfoils. The concept was extended by Tang et al. [8] that considered an array of four airfoils. Approaches using arrays of six or more vanes (e.g. [9,10]) have been tested, though the wakes of the vanes introduced turbulent fluctuations into the downstream flow conditions. Simpler generation mechanisms involving two pitching plates (e.g. [11,12]) avoid wake effects by construction, but tend to be limited in both the amplitude of the gusts produced and reduced frequency they can achieve. However, this design has been proven to be effective in transonic wind tunnels ( [13]). Moreover, most of the studies cited here have considered traditional shapes of airfoils belonging to the NACA family.
The present work has been carried out in the frame of the H2020/Clean Sky 2 project called GUDGET (grant agreement No 831802, [14]). The objective of the project was to design, manufacture, calibrate, verify and finally install in the ONERA S3Ch Wind Tunnel (WT) [15] an enhanced gust generator system able to generate gusts with amplitudes up to +/-1 deg of flow deviation, in the frequency bandwidth [0-100 Hz], for a transonic flow field with WT airspeed Mach number equals to 0.82 and an aeroelastic half-model connected to the WT side wall. Within GUDGET project, a preliminary trade-off analysis has been performed to find the best configuration of the gust generator, by considering innovative configurations of tilting airfoils moved by mechanical actuators as well as blowing slots fed by fluidic actuators or a combination of both. In this framework, starting from the baseline shape of a NACA 0012, the present paper focuses on the shape optimization of the airfoil section of the gust generator in transonic conditions by using a design of experiment (DoE) procedure. The apparatus consists in a tandem of airfoils installed horizontally right upstream of the test section. The generator is set so as to deliver deterministic vertical and harmonic gusts in a flow dynamically scaled in Mach number. The objective of these tests is to qualify the gusts in the empty test section, in terms of amplitude, spatial and temporal form. Mach number 0.73 has been considered to validate CFD approach in pre-design phase using data provided by ONERA [16] and a Mach number of 0.82 in the test chamber has been considered during the design/ optimization phase for the aerodynamic shape design of the two vanes. Two different strategies have been considered for the generation of the desired gusts: -Rotating vanes -Fixed vanes with fluidic actuators generating a sonic jet at the vane trailing edge leading to two different optimum airfoil shapes. The characteristics of the fluidic actuators, the oscillation amplitude, frequency and maximum flow rate are not subject to the optimization process. The paper is organized as follows. In Sect. 2 the overall computational model is presented, including the geometries object of this study. The details of CFD formulation are given in Sect. 3 along with the approach used to validate it. The final results are presented in Sects. 4 and 5; finally, some concluding remarks are reported in Sect. 6.

Design of Experiment (DOE)
Design of experiment (DOE) is a procedure for choosing a set of samples in the design space, with the general goal of maximizing the amount of information gained from a limited number of samples. This representative set of data is then used to compute a response surface (Response surface method) to perform optimization of a baseline design. This technique is commonly used for the aerodynamic design ( [17,18]), and is useful for analysing the complex correlations between the geometrical parameters and the performance of the fluid machinery.
Central composite design (CCD) has been used as DOE technique for the optimization of vanes' airfoil .
A complete overview of the methodology can be found in [19][20][21] and [22]. Figure 1 shows the computational domain considered for the CFD calculations needed to properly estimate the baseline performance and to collect the aerodynamic data to perform the optimization. The gust vane device is installed right at the beginning of the wind tunnel test section and is composed of two movable 2D wings fixed on the walls of the wind tunnel, whose chord length is 200 mm. Available investigations, reported in [13,23], suggest that the choice of a 2D model gives a good agreement with the experimental data and a full description of the response of the system. The baseline airfoil shape of the two vanes is a NACA 0012. The airfoil, together with the points in the test section where the gust angle is evaluated (indicated as x 1 , x 2 , and x 3 ), are shown in Fig. 1. All the units in Fig. 1 are in mm.

Baseline Geometry
The gust angle, which from now on will be indicated as , is the main performance parameter. It is defined as follows: where v is the vertical velocity profile and u is the horizontal velocity profile at the control point. The wind tunnel layout is not subject to the optimization process. Figure 2 shows the baseline airfoil layout that, for the rotating vane design, rotates around the point at one quarter of the chord. For the configurations with fluidic actuators, two jet slots are present at 95% of the chord. The pitch motion type is a harmonic sinusoidal waveform with frequency 40 Hz and maximum amplitude ±2 deg, the double injector (upper and lower side) has harmonic synchronous alternating half-sine shape with frequency 40 Hz in terms of mass flow rate with maximum equal to 120 g/s/m.
More details about the choice of the baseline parameters can be found in [16]. Figure 3 gives a better understanding of pitching and blowing jet law. The pitch angle is considered positive when the trailing edge goes actually up (nose down).
In Fig. 3 the values of pitch angle must be read on the left and mass flow values for upper and lower slot on the axis on the right. If the vane rotation effect and the jet provided by fluidic actuators are combined, when the maximum positive pitch angle is attained, the upper slot reaches  It has to be noted that the values of the slot position, oscillation amplitude, frequency and maximum flow rate are given as an input by ONERA. For this reason, the only objective of this paper is to show how to improve the performances of the system by designing a new shape of the gust generator by considering a classical configuration. A different design process has been performed to optimize these parameters. To simulate the jet flow, the advanced capabilities of the CFD solver employed (ANSYS Fluent) have been used. In particular using a user-defined function, or UDF, the boundary condition at the bottom of the slot visible in Fig. 2 is iteratively changed from mass flow inlet (when the jet is active) to wall (when the jet is not active) according to the jet flow profile already presented in Fig. 3.

Computational Grid
Figures 4 and 5 show the computational grid developed for the baseline performance calculations, and created using ANSYS Meshing along with all the other computational grids used in this work. For the performance evaluation of the system with moving vanes, a sliding mesh approach has been used and two circular rotating sub-domains considered, which are clearly visible in Fig. 4. For the system with non-moving vanes and fluidic actuators, the sliding mesh approach is not necessary. The mesh type considered is a hybrid mesh with prismatic layers to properly simulate the All the computational grids used for the optimization are consistent with the ones presented in this section (which represents the baseline mesh), and all the quality parameters (skewness, aspect ratio) are within the required range.
To ensure the chosen mesh is able to describe the physical phenomena under study a sensitivity analysis has been performed. The mesh investigation carried out in this work is based on Grid Convergence Index from Celik et al. [25]. Basically, four different grids with different grid refinement have been created: coarse, medium, fine and ultrafine mesh (see Table 1).
For the calculations a Mach number equal to 0.82 has been imposed inside the empty test chamber (see Fig. 1), the Reynolds number referred to a reference length equal to the unity (the test chamber is actually empty) is 19 * 10 6 . The initial turbulent intensity and viscosity ratio at the inlet are 0.1 per cent and 0.001, respectively.
Considering the performance output for this engineering system is the gust angle, and consequently the vertical and horizontal velocity component profiles (evaluated at the check points x 1 , x 2 , and x 3 ), to perform the procedure described in [25] a series of steady calculation has been considered at first to evaluate the velocity magnitude at control points.
In the steady state calculations, the gust generator system is, of course, not active. Results are available in Table 2.
Indicating medium, fine and ultra-fine mesh as grid 3, 2 and 1 respectively. The Grid Convergence Index has been calculated according to the following expression: where Fs is a safety factor (considered equal to 1.25 in this work), e is the error between the two grids and r is the ratio of refinement (equal to 2 in this particular case). p is evaluated using: in which f i is the chosen output parameter (velocity magnitude) and the subscript refers to the specific mesh considered, for instance, i = 1 refers to what we called ultra-fine mesh while i = 3 to the medium mesh. Eventually, Eq. 3 leads to the value of 2.11. A further verification consists in evaluating of the ratio: In Eq. 4 the subscripts i and j refer to the specific mesh considered, the single index goes from 1 to 3. For instance, if i = 1 and j = 2 we are referring to the mesh called ultafine (mesh 1) and fine mesh (mesh 2), and GCI ij is calculate according to the output values obtained using the computational grids the subscripts refer to (ultra-fine and fine mesh in this example). Table 3 shows the Grid Convergence Index calculated starting from the velocity magnitude for each one of the three control points. Based on velocity at control points, it appears the asymptotic range has been reached (Table 3).
A further verification has been made considering an actual unsteady calculation, medium and fine meshes were considered and the vertical velocity profile extracted from control point x 2 . Figure 6 shows the time-dependent vertical velocity v profile when both rotation and jet are considered. For reasons of clarity, the results are reported only for medium and fine meshes since fine and ultra-fine mesh show almost coincident results and coarse grid produces unreliable results. It can be seen that there is a slight underestimation of the ,  maximum peak of velocity for the medium mesh (below 3 per cent) after the initial transition has ended up (say t > 0.04s).
Considering the general good behaviour of the medium mesh, it has been used for the optimization calculations to have a good compromise between numerical efficiency and accuracy, as a high amount of design points had to be analysed. Since the gust angle peak is the main objective, a slight underestimation of the peak has been considered as not affecting the optimization as long as all the calculations refer to the same flow and movement condition.

Moving Vane with Fluidic Actuators
Following indications by ONERA, the CFD methodology has been validated using CFD data available in [16] for Mach number 0.73 inside the test chamber obtained with a high accuracy structured mesh. Considering as the most challenging case, we decided to validate our calculation using the data for the case with effect of rotation and fluidic actuators combined. Anyway, further verifications have been made using experimental data available for a simplified case where only the effect of rotation was considered [13]. In particular, percentage differences below 3% have been found with respect to the angular excursion of the moving airfoils.
As already pointed out, the pitching frequency is 40 Hz with a pitch deflection range of ±2 deg and a harmonic pitch motion type, the jet profile is harmonic synchronous alternating half-sine with maximum flow rate 120 g/s/m (see Fig. 3).
Full unsteady RANS approach has been used for this calculation alongside with the sliding mesh methodology available within the capabilities of ANSYS Fluent solver for this particular case. A suitable time step equal to t = 1 * 10 −5 s has been considered and the k-SST turbulence model has been used. In terms of numerical approach, a second-order discretization has been used both in time and space. The solver used for all calculation is a pressure-based coupled solver. Pressure inlet boundary condition is specified at inlet while pressure outlet is specified at outlet. The actual values specified allows the flow to reach the required Mach number inside the test chamber.
Results in terms of gust angle extracted at control points are given in Fig. 7, and show that the initial transient regime is damped after about two cycles leading to a perfectly periodical signal.
For comparison with the reference, data extracted at the control point x 2 have been used and the initial transient has not been considered in this case. One can see the general good agreement of numerical results especially in terms of signal shape and average maximum gust angle peak. Figures 8 and 9 refer both to a condition when the maximum jet mass flow is attained at upper jet slot. From those figures, it can be seen how a sonic condition for the jet corresponds to the maximum mass flow.

Design Process
The first step of the design process was the choice of the airfoil family to generate the samples. The starting point (initial airfoil shape) has been represented by the NACA 0012 airfoil typically adopted for the design of gust generators. The result of the process was the generation of NACA 4-Digit Modified symmetric airfoil.
Geometry has been parametrized considering three major design parameters: designation of the leading edge radius, I , strictly related to leading edge radius ( r le ), chord-wise position of maximum thickness (T -in tenths of chord) and the maximum thickness itself ( t k ). The output (performance) parameter was the maximum gust angle (gust angle peak). The corresponding equations for leading edge radius and maximum thickness for NACA 4-Digit are: where, in Eq. 5, I is the designation of the leading edge radius, valid as long as I ≤ 8.
The complete procedure aimed at solving these equations and evaluating the coefficients to generate all the airfoil geometries is fully described in [26] and for the sake of brevity it is not repeated in this paper.
Once the geometry has been parametrized, a suitable design matrix has been defined considering the chosen DOE methodology (CCD) to run the numerical experiments. Considering the number of design parameters that have been considered (I, T, t k ), some extra points have been added to the initial matrix (which should have a minimum of 15 points) to refine the response, leading to the actual number of points employed during the process and described in the following sections of this work. For all design processes carried out in this work the design boundaries are summarized in Table 4.

Moving Vane
This section is dedicated to the shape design considering rotating vanes and inactive jet fluidic actuators. A series of samples have been generated according to the design space shown in Table 4 and the Central composite design DOE scheme. It has to be note that in considering the boundaries of the design space, aspects related to the mechanical design, which is not considered in this work, have led to use a design space not symmetrical with respect to the thickness t k . The design objective is to maximize the performance of the gust generators in terms of gust angle inside the test chamber where control points are placed, as already mentioned vertical and horizontal components are extracted to evaluate the temporal profile of the gust angle. Average peak gust angle is considered as output parameter (overall performance parameter). Inside the test chamber, a Mach number equal to 0.82 is imposed.
To generate the model, 22 design points have been analysed, results in terms of response are presented in Fig. 10 where the red dot identifies the baseline.
A zone which corresponds to the high-performance designs can be clearly identified, which leads to designs with a gain in terms of maximum gust angle higher than 15 per cent.
Among the design parameters under consideration, a strong influence is exerted by the position of the maximum thickness (see Fig. 10b). It has been found that the leading edge radius has considerable influence as well, although moving the maximum thickness position towards the leading edge seems to be the key aspect to increase the performance in terms of gust angle inside the test chamber.
In terms of influence of the maximum thickness t k the data analysed suggests that the best performance corresponds to intermediate thickness with respect to the design space considered for this parameter (see Fig. 10b).
Following the outcome presented in Fig. 10, a possible candidate point has been analysed using the fine mesh developed during the sensitivity analysis phase ( 280 * 10 3 elements).  The results of this analysis is presented in Figs. 11, 12, 13, 14 in terms of gust angle contour inside the test chamber compared with results obtained for the baseline design. Four different time instants have been chosen t = 0.06 s , t = 0.072 s , t = 0.085 s and t = 0.097 s which correspond to gust angle peaks (positive and negative) as seen by control points. One can see how there is a considerable improvement both in gust angle magnitude and extension of the high gust angle zone inside the test chamber. Fig. 15 shows the airfoil geometry corresponding to the proposed potential candidate point compared with the initial baseline (NACA 0012). The main difference between the two geometries proposed is actually the maximum thickness position which is moved forward for the candidate point with respect to the baseline.
To complete the analysis about the sensitiveness of the maximum thickness position parameter in terms of gust angle performance, an additional simulation has been performed by moving the maximum thickness position of the candidate best design point towards the leading edge, i.e. T = 0.15 instead of T = 0.2 . In this configuration, the parameters are I = 6 , T = 0.15 and t k = 0.121 . Even if this configuration has shown higher values of maximum gust angles compared to the one having T = 0.2 , a large number of disturbances in the gust angle signal has been observed: Fig. 16 shows the comparison in terms of pressure signal between the T = 0.15 and T = 0.25 by considering the reference time that produces the maximum positive gust angle in the test chamber.
It is interesting to analyse the pressure coefficient distributions on the airfoils located below the centerline of the wind tunnel when the maximum positive (green curves) and zero (red curves) gust angle conditions are reached for, respectively, T = 0.15 and T = 0.25 inside the test chamber  Fig. 17 for details). Of course, there is a phase delay between the time instant when the maximum gust angle is measured at control points inside the test chamber and the time instant corresponding at the pressure distribution on the airfoils that actually generated such a gust. This time difference has been considered when generating the pressure coefficient plots presented in Fig. 17. For T = 0.25 , the comparison between red and green curves shows that the shock wave on the upper part of the airfoil oscillates around X = −0.6 m going from the max positive to the zero gust angle instants. For T = 0.15 , it is possible to observe that the shock wave has a different position towards the leading edge of the airfoil and it covers a greater distance than T = 0.25 during the oscillation of the airfoil in the same time interval. Moreover T = 0.15 configuration shows a second shock wave on the lower part of the airfoil, that oscillates too. This situation can justify the noisy signal of the pressure and the higher values of the maximum gust angle observed.

Fixed vane with fluidic actuators
The second aerodynamic design process takes in consideration fixed vanes while the gust is generated by alternate jets of fluid. In this case a double injector is present (upper and lower side) and the pitch angle is set to zero (aifoil perfectly horizontal). The upper and lower bounds of the design space are unchanged with respect to the optimization of rotating vanes and central composite design (CCD) has been used as well. The number of design points analysed is comparable with the design performed for the rotating vane and presented in the previous section, although more points have been added in this case (29 points). Note the law of variation of mass flow rate for the jet slots is the one provided in Fig. 3.
It is known for the case with fixed vane and pulsed jet provided by fluidic actuator that the output in terms of gust angle is proportional to frequency, chord length and the derivative of the lift coefficient with respect to mass flow rate [16]. More in details, the following expression is valid: where f is the frequency, c is the chord length, V ∞ is the stream velocity inside the test chamber and C l m is the derivative of lift coefficient with respect to mass flow rate.
Considering C l ≈ C y in a condition of fixed vane position (perfectly horizontal vane), the design with the best performance in terms of gust angle should provide a higher C y as well when plotted versus time.
Two of the best design points on the best performance zones in Fig. 18 have been analysed using the same approach already presented for the moving vane (fine mesh used for better accuracy). We will refer to them as BDP 1 ( I = 4 , T = 0.4 , t k = 0.12 ) and BDP 2 ( I = 7 , T = 0.3 , t k = 0.125 ), the results in terms of gust angle contour inside the test chamber are presented in Figs. 19, 20 and 21, 22 for two time instants which correspond to maximum and minimum gust as seen by control points.
The upper side of Figs. 19, 20 and 21, 22 which corresponds to the performance of BDP 1 and BDP 2, respectively, clearly shows a good deal of improvement with respect to the baseline which is represented on the lower side of the same picture.
Comparing the variation of the vertical force coefficient versus time for the two best design points to the values obtained for the baseline design (Fig. 23) it can be seen how especially BDP 1 shows significantly higher C y having consequently a higher C l m . To conclude this section, the geometries of BDP 1 and 2 are compared in Fig. 24 to the baseline. It is evident how the geometry BDP 1 and BDP 2 provide more room at the TE, On the other hand, BDP 2 offers the advantage of the bigger thickness which impacts favourably on the mechanical design as well.

Performance at 80 Hz
Considering the simple model proposed in [13] for a moving vane gust generator, the performance in terms of gust angle can be expressed as follows: where A g is the pitching amplitude, dC L0 d is the slope of the lift characteristic of the airfoil, C(S t ) is a correction to the static lift which is called the lift-deficiency function, S t is Strouhal number and the gust maximum angle is indicated with A.
Looking at (9) it is intuitive how the Strouhal number and thus the frequency plays an important role for the gust generator performance. Although the model itself is relatively simple to be rigorously applicable to the systems under study, a considerable improving of the performance is expected increasing the frequency of either the moving vane or the jet or both for the last system under investigation.
To verify the behaviour of the candidate points produced during the design phase, their performance in terms of gust angle produced inside the wind tunnel test chamber has been evaluated at a frequency of 80 Hz. Results are presented in terms of gust angle contour inside the test chamber.

Candidate point-moving vanes
As expected, the improvement of the system is significant at a frequency of 80 Hz. In Figs. 25, 26 and 27, it should be noted that the time instants considered refer to minimum and maximum peak as seen by control points. That is the reason why the contour plot of candidate point and baseline are not taken exactly at the same time. Moreover, due to the practical difficulty of generating a movement of the vanes of ±2 • at 80 Hz, this working condition turned out to be not feasible into the WT. However, the optimization showed its effectiveness also at higher frequencies, as foreseen by the theory. Figure 28 shows the effect of doubling up the frequency on what in the corresponding section has been called BDP 1 ( I = 4 , T = 0.4 , t k = 0.12 ). As one can see, this design ensures a maximum gust angle order of magnitude about 0.8 degrees inside the test chamber. This result turned out to be very useful for the project, as high-frequency gusts are going to be generated inside the WT using fluidic actuators only (non-moving vanes).

Conclusions
Two different aero-shape designs have been performed considering different vane configurations. Although a good deal of improvement can be obtained re-designing the aero-shape it has to be taken into account the fact the aero-shape alone is not capable of drastically changing the magnitude of the gust angle peak. As a consequence, if a considerably higher performance is required (say 1 • as gust angle peak), other design aspect, in addition to the aero-shape changes must be taken into account such as jet inclination and position for non-moving vanes at high frequencies or the possibility to include a flap to each vane for the moving vanes at low frequencies.  The possibility to reach such a high gust are increased considering a vane which combines both rotation and jet effect, but of course the high difficulty to design such a systems can make it not feasible or so expensive and complex to be considered not attractive.

Disclaimer
The results shown in this paper reflects only the author's view and the EU Commission and Clean Sky 2 Joint Undertaking are not responsible for any use that may be made of the information the paper contains.
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/.