#### Abstract

The size of aircraft models that can be tested in icing wind tunnels is limited by the dimensions of the facilities in present; it is an effective method to replace the large model with a hybrid airfoil to carry out the experiment. A design method of multiple control points for hybrid airfoil based on the similarity of flow field in the leading edge of airfoil is proposed. Aiming at generating the full-scale flow field and ice accretion on the leading edge, multiobjective genetic optimization algorithm is used to design the hybrid airfoil under different conditions by combining the airfoil parameterization and solution of spatial constraint. Pressure tests of hybrid airfoils are carried out and compared with the leading edge pressure of the corresponding full-scale airfoils. The design and experimental results show that the pressure coefficient deviation between the hybrid airfoils designed and the corresponding full-scale airfoil in the 15% chord length range of the leading edge is within 4%. Finally, the vortex distribution and ice accretion process of the two airfoils were simulated by the unsteady Reynolds-averaged-Navier–Stokes (URANS) equations and multistep ice numerical method; it is shown that the hybrid airfoil can provide the same vortex distribution and ice accretion with the full-scale airfoil.

#### 1. Introduction

Icing encounters in aviation are common; when the aircraft operates in a cold and humid environment, supercooled droplets collide with the aircraft surface, resulting in ice being accreted on the surfaces of the aircraft [1]. The ice depletes the overall performance of the aircraft and significantly threatens flight safety if not properly mitigated [2]. To improve flight safety, one of the main tasks is to evaluate the performance of iced airfoils and aircraft, ensuring they can operate safely even in the worst icing conditions as set by the FAA [3]. There are several methods to assess iced aerodynamic performance of airfoils: conducting simulations, flying under natural conditions, and icing wind tunnel test. Among these methods, icing wind tunnel test is cost-effective, safe, and reliable. An icing wind tunnel can create an environment similar to the icing conditions in nature, and then, research experiments can be conducted and the icing process studied.

Theoretically, we can obtain the most accurate experimental results if a full-scale model is used for icing tests. However, icing wind tunnels are too small to test full-scale airfoils or wings; even for a moderate- or small-sized model, it may cause tunnel blockage [4] or wind tunnel wall interference when the model’s angle of attack (AoA) is adjusted. To weaken the blockage and acquire ice characteristics of the full-scale model under existing conditions, scaled model is used to substitute the full-scale model. There are two methods to scale models [5]: first, geometrically scaling, a full-scale airfoil or wing is scaled to a certain scale and then scaling the icing condition under the criterions of AEDC [6] and ONERA [7] to get icing similitude in existing facilities. But the usefulness of this method is limited and only valid with moderately sized models and scale factors. The second scaling method called “hybrid airfoils” is useful in the icing wind tunnel tests. By this method, an airfoil is designed with full-scale leading edges and a redesigned aft section to reduce the chord and provide full-scale icing condition in the leading edge region. A hybrid airfoil has a shorter chord length and can reduce the model size and tunnel blockage.

Hybrid airfoil research began in the 1950s by Glahn [8] and then the first systematic hybrid airfoil design method established and improved by Saeed et al. [9–11]. This was based on the fact that leading edge ice accretion will be the same between the hybrid and full-scale airfoil if cloud properties, droplet impingement, local leading-edge flow field, and model surface geometry are kept the same. In addition, ice usually only accretes on the leading edge of airfoils and so hybrid airfoils are designed such that it has a similar geometry to that of the leading edge of the full-scale airfoil. The aft section is redesigned to provide a similar flow field on the leading edge with the full-scale airfoil. With this approach, icing is in good conformity between the full-scale airfoil and hybrid airfoil; the greatest benefit of the hybrid airfoil design method is that it can obtain a similar icing process of the full-scale model with a hybrid airfoil without the parameters of the experiment changing, and this is very convenient for experiments. Fujiwara et al. [12, 13] completed research based on the work of Saeed, which employed an icing accretion code to perform comparisons of the ice shape between full-scale and hybrid airfoils; the sensitivity of different design parameters on hybrid airfoil performance as well as how circulation changes from the full-scale to the hybrid airfoils was analyzed [14]. Viscous and tunnel wall effects were checked with the use of 2D Reynolds-averaged-Navier–Stokes (RANS). Some useful conclusions have been drawn in his studies for hybrid design: a similar flow field, aerodynamic, and heat transfer characteristics between hybrid and full-scale airfoil are necessary; matching stagnation point is the first-order parameter, while suction peak magnitude is of second order; some geometric parameters such as nose droop angle, leading-edge extent, and scale factor have been evaluated for hybrid airfoil design [15]; he also pointed that compromises and trade-offs are required in the design process. See Refs. [12–15] for details.

Through the previous studies, Saeed proved the feasibility of the hybrid airfoil. Fujiwara confirmed that the similar flow field of hybrid airfoils resulted in a similar ice accretion to that of the full-scale airfoil. Also, hybrid airfoil design is based on matching the aerodynamic flow field between hybrid and full-scale airfoils; in this design method, take stagnation point location and suction peak magnitude as the indicators that the ice accretions between full-scale and hybrid airfoils will be similar.

However, in icing processes, the droplets are affected by the flow field in the near-field region and impinged at the leading edge of airfoil. Thus, in order to accurately design a hybrid airfoil, a new hybrid airfoil design method is presented, which is based on the previous work by Fujiwara and Bragg et al., and focuses on matching the full-scale aerodynamic flow field of hybrid and full-scale airfoil. In this method, a reference area is built and used; also, a multicontrol point optimization method based on MOGA (multiobjective genetic algorithm) is established. A high-fidelity flow solver based on 2D RANS methodology is used to calculate the flow field in the whole optimization. Also, a parameterization hybrid airfoil expression method is presented for the convenience of representing the different hybrid airfoils. Finally, the effect of design method was proved through numerical calculation and pressure tests.

#### 2. Parametric Hybrid Airfoil Design Method

##### 2.1. Definition of Hybrid Airfoil

An ideal hybrid airfoil design for icing wind tunnel testing will produce the same effect as full-scale icing; based on the results obtained from previous papers, a different hybrid design method is proposed; for different full-scale airfoils, corresponding hybrid airfoil can be designed under the specified leading edge range and scale factor. Scale factor (SF) of hybrid airfoil is defined as the full-scale chord divided by the hybrid airfoil chord. As shown in Figure 1, dashed line represents the full-scale airfoil of NACA0012, with a dimensionless chord length 1, and the solid line is the hybrid airfoil with a chord length of 0.5; thus, the SF of hybrid airfoil is 2 in here.

In this paper, 25% chord length of full-scale airfoil is taken for both upper and lower leading edge extents; as shown in Figure 1, hybrid airfoils have the same leading edge with full-scale airfoil, but a deformable aft section that is controlled by several points with the parameters of these control points varies in the optimization process, and the shape of hybrid airfoil also changes. Also, this leading extent can be changed by the designer according to different design conditions.

##### 2.2. Hybrid Airfoil Optimization

In this section, a Sum of Sine parameterization method for hybrid airfoil is introduced, and multiobjective genetic algorithm process for hybrid airfoil design is established.

###### 2.2.1. Optimization Algorithm

A conceptual illustration of the hybrid airfoil design approach is shown in Figure 2. There are two optimization loops: the first loop (method of Saeed) is based on XFOIL and GA, which is used to obtain an initial hybrid airfoil. The second loop is based on the initial hybrid airfoil get from the first loop, and Pointwise, 2D RANS methodology, and MOGA are employed to obtain the final hybrid airfoil. The specific design method is as follows. (1) Considering that usually the icing area of the wing is 10% of the chord length, take the 25% leading edge chord length of the full-scale airfoil as the leading edge of the hybrid airfoil, called nose section; this part of the full-scale airfoil is fixed in the hybrid airfoil optimization. (2) The aft section is designed to provide full-scale flow field on the nose section of the hybrid airfoil; it is generated by defining the control points arbitrarily. Then, a primitive airfoil is combined by nose section and aft section with the SS (Sum of Sine) and PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) methods. The designed aft section and primitive airfoil are controlled by several constraints, such as the scale of the hybrid airfoil and the continuity condition of second-order derivative. (3) Start the first loop; Cps (pressure coefficient) are calculated by XFOIL in view of its quick and easy features. (4) Compare the Cp of full-scale and hybrid airfoils; genetic algorithm optimization is used to update the aft section of the hybrid airfoil until the Cp of the nose section matches that of the full-scale airfoil. (5) Start the second optimization loop. It is the same procedure as the first optimization loop, but the difference is that it is more accurate because of the use of automatic grid generation based on Pointwise Glyph Script and CFD. (6) Compare and output the hybrid airfoil; ensure that the flow field of the hybrid and full-scale airfoils within the reference area meet the requirements; otherwise, continue the optimization with GA. The ultimate hybrid airfoil is obtained after completion of the second optimization. By changing the parameters, the hybrid airfoil at different AoA can be obtained.

###### 2.2.2. Parameterization for Hybrid Airfoil

The hybrid airfoil designed in this paper is different from common airfoils, as it has a fixed leading edge, a flat upper surface, and a trailing edge with large curvature. An inverse design method based on potential theory is not applicable to hybrid airfoils [16]. Thus, genetic algorithm is employed. In order to implement an automatic process of hybrid airfoil optimization, a parametric hybrid airfoil expression method is needed, selecting the appropriate parameters and transforming the hybrid airfoil into a function containing these parameters; the hybrid airfoil will deform by adjusting these parameters.

There are many parametric expressions to describe airfoils. (1) Analytic function methods, such as the Joukowsky [17] airfoil, are expressed by mean camber and thickness functions. This method is more suitable for static airfoils. (2) Bézier curve approximation algorithm is a method using many control parameters and requires more computing resources. (3) The Hicks-Henne [18] shape function is widely used in airfoil parameterization; this method is highly dependent on the original airfoil and has little ability to change the trailing edge angle.

According to the characteristics of hybrid airfoil, three parametric functions, Spline, Bezier, and Sum of Sine, are selected to fit the upper and lower surfaces of the hybrid; the fitting error is shown in Table 1; it can be seen that the fitting error is 10^{-4} by Spline function, 10^{-3} by Bezier method, and 10^{-6} with Sum of Sine method. Fitting shapes for different methods are presented in Figure 3. Obviously, the Sum of Sine method has the highest fitting precision. So it is selected in this paper to parameterize the hybrid airfoils.

**(a)**

**(b)**

**(c)**

**(d)**

###### 2.2.3. Reference Area Selection

Aerodynamic characteristics determine the motion of droplets, and impingement characteristics are responsible for the interaction of the droplets with the airfoil and ultimately affect the ice shape and range. The same aerodynamic and impingement characteristics of full-scale and hybrid airfoils are important principles for hybrid airfoil design. The same Cp at the leading edge of hybrid and full-scale airfoil was used as figures of merit for optimization in previous studies. However, in our opinion, icing is the result of water droplets moving along with flow field, impinged on the airfoil, heat transfer, and condensation in this region. Therefore, in order to design the hybrid airfoil accurately, the leading edge area should be included as the reference area in optimization process; the hybrid airfoil should be designed so that the aerodynamic characteristics of this region are the same as those of the full-scale airfoil. The reference area adapted is shown in Figure 4; it is C-H shape, and it consists of a semicircle with a radius of 0.1 c and a rectangle with a chord direction of 0.15 c and a normal direction of 0.2 c. The black points are the reference points, they are used to discretize the reference area, and the information value of these points represents the values of the whole reference area, which will be introduced into the genetic optimization to establish the fitness function and find the airfoil with the maximum fitness in the optimization.

###### 2.2.4. Definition of the Optimization Problem

Genetic algorithm (GA) is a global search algorithm according to Darwin’s theory of survival of the fitness and natural selection; it is widely used in computer science, engineering optimization, airfoil design, and other fields. At first, the optimization problem was simple and single objective; GA can get good results; however, with the increasing complexity of optimization problems and multiobjective, the simple GA cannot reach a satisfactory result. In this case, MOGA is more suitable; in this paper, the hybrid airfoil is optimized with the use of MOGA [19].

Seven control points and fourteen design variables are used to form a hybrid airfoil profile. Thus, the profile of hybrid airfoil has more degrees of freedom [20]. At first, the leading edge of the full-scale airfoil is taken as the hybrid airfoil’s leading edge and connected with an arbitrary aft section. Then, the flow field of hybrid and full-scale airfoil is obtained by using a RANS solver; compare the Cps of the two airfoils; the difference is shown in Figure 5. The design objectives are to minimize the Cp difference and SP (stagnation point) difference between full-scale and hybrid airfoil in the reference area. The shape of hybrid airfoil is changed in the MOGA process to get the individual with the maximum fitness. The genetic algorithm source code is built by authors, which is based on MATLAB and runs on PC.

The definition of optimization can be summarized as follows:

Objective 1:

Objective 2:

Meanwhile, fitness function is defined as where is reference area, subscript of f and h represent full-scale and hybrid airfoil, respectively, represent position of upper and lower surface, and represent position of upper and lower surface, respectively.

The purpose of genetic optimization is to find the individual who has the maximum fitness in evolution process by genetic optimization; in this genetic algorithm optimization, set the population size 30, the evolution generation number 30, the crossover rate 0.85, and the mutation rate 0.05.

#### 3. Results and Discussion

##### 3.1. Calculation Validation

An appropriately sized grid can increase the efficiency of CFD as well as the optimization. Grid topology of airfoils used is shown in Figure 6; the computational domain is 10 c above and below the airfoil and 10 c wide upstream, the far field extends 20 c downstream of the airfoil, and the initial grid spacing satisfies [21].

With this grid, the N-S equations are solved by numerical method, here is ANSYS FLUENT. The governing equations of continuity, momentum, are solved by the fluid solver, expressed as

The Spalart-Allmaras turbulence model is also used. The S-A model [22] is a one-equation model: where is turbulent eddy viscosity, is molecular viscosity, is destruction term, is production term, and where .

The NACA0012 airfoil of chord length 1 is used to validate the numerical method, where airfoils of NACA0012 at different AoA are compared to those of the published experimental data [23]. The parameters for verification calculation are the same as the experiment. Figures 7(a) and 7(b) show the pressure coefficient at 0° and 4° AoA, respectively. Figure 8 demonstrates the Cl and polar curve obtained from numerical method and experiment. As we can see, the Cp distributions have good conformity as well as the Cl and polar curve.

**(a)**

**(b)**

**(a)**

**(b)**

##### 3.2. Hybrid Airfoils on Design Conditions

With the method described above, different hybrid airfoils from 0° to 10° AoA were designed. As genetic algorithm is used in hybrid airfoil design, the population size is 30 and evolutionary generation is 30. The PC’s configuration is as follows: CPU is Intel® Core™ i5-8004, 16 GB of memory. It takes about 7 to 9 hours to get a hybrid airfoil with this optimization process; the time will be reduced if you do it with a high-performance computer or improve the algorithm.

The profile and pressure distribution comparison of hybrid airfoils at different AoA is illustrated in Figure 9; the magnified box of the same figure shows the reference area of interest (suction peak and stagnation point) in more detail; as we can see, the hybrid airfoil shows a good consistency of pressure coefficient with that of the full-scale airfoil. It is clear that the hybrid airfoil designed here has a perfect match in stagnation points as well as pressure coefficient with the full-scale airfoil. And this match proves that the aerodynamic characteristics of the hybrid airfoil are the same to the full-scale airfoil, which is the important criterion in the optimization process.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

#### 4. Test on Aerodynamic Characteristics of Airfoils

An experiment was conducted to test the aerodynamic characteristics of the hybrid airfoils obtained by the design method presented in this paper and also to check the effectiveness of the method.

##### 4.1. Quantitative Comparison of Hybrid and Full-Scale Airfoil

In this optimization, the reference area is shown in Section 2.2.3; the flow field characteristics in the reference region of hybrid and full-scale airfoil are compared quantitatively to show the design performance of this method. The flow field was calculated by CFD, and the reference points are list in Figure 4; the mean velocity errors of direction and direction of hybrid and full-scale airfoil of different reference points are 4.12% and 2.93%, respectively. So, we can know that, the hybrid airfoil designed by this method could provide the same aerodynamic characteristics with the corresponding full-scale airfoils.

##### 4.2. Pressure Coefficient Measurement

As described in Section 4.1, the performance of hybrid airfoil was compared by CFD; in order to check the aerodynamic characteristic of hybrid and full-scale airfoils, pressure measurement test was also employed. In this test, 10 different hybrid airfoils are selected and divided into two groups according to their corresponding full-scale airfoils: (1)Hybrid airfoils based on the full-scale airfoil NACA0012 at 0°, 2°, 4°, 6°, 8°, and 10° AoA, named Hy00, Hy02, Hy04, Hy06, Hy08, and Hy010(2)Hybrid airfoils based on the full-scale airfoil NACA2412 at 0°, 2°, 4°, and 6° AoA, named Hy20, Hy22, Hy24, and Hy26

In this test, the NACA0012 and NACA2412 airfoils both have a chord length of 300 mm and the different hybrid airfoils with a chord length of 150 mm. Each airfoil has 31 pressure taps, which are arranged in three rows and located on the upper and lower surfaces. The pressure taps are concentrated at the leading edge and gradually dispersed to the trailing edge with a broken line distribution. Figure 10 shows the location of pressure taps on the upper surface of the hybrid airfoil model. The pressure measurement test was carried out in the LTWT (Low Turbulence Wind Tunnel) of Northwest Polytechnic University. This wind tunnel is a low-speed and extremely low turbulence direct current research wind tunnel with the international advanced level and a minimum turbulence of 0.03%. It can provide an experimental environment with excellent flow field quality for the research of aircraft and other industrial aerodynamics. The layout of the wind tunnel is shown in Figure 11.

As shown in Figure 12, the airfoil was horizontally installed in the wind tunnel, and a pitot tube was installed above the trailing edge to facilitate subsequent data processing and pressure correction. End plates were used to make the flow field as two-dimensional as possible. The piezometers were led out from the pressure taps and connected to DSY104 electronic scanning micropressure measuring system. There are 192 pressure measuring channels in the system, the measuring ranges are and , the measuring precision is less than ±0.1% full scale, and the scanning rate is 20,000 points/sec. Measuring range of ±7.5kPa was selected in this experiment.

##### 4.3. Pressure Measurement Results

In order to validate the design performance of the aforementioned hybrid airfoils, the leading edge pressure measurement experiment was designed. For different airfoils, experimental conditions are shown in Table 2.

After data processing, the comparison of pressure coefficient between hybrid and full-scale airfoils is shown in Figure 13; and represent the error between the pressure coefficient distribution on the upper and lower surfaces of the hybrid airfoil and the corresponding full-scale airfoil, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

The errors of these 10 cases are listed in Table 3; obviously, the upper and lower Cp of the hybrid airfoils and full-scale airfoils at the leading edge are essentially the same; the average error is less than 4%.

#### 5. Ice Accretion and Ice-Induced Vortices of Hybrid and Full-Scale Airfoils

##### 5.1. Ice Accretion

To obtain a hybrid airfoil which generates the same ice accretion as a full-scale airfoil for icing wind tunnel testing is the purpose of this hybrid airfoil design method. Pressure tests have confirmed that the Cp of the two airfoils at the leading edge is the same. On this basis, it is necessary to analyze the icing process of hybrid and full-scale airfoil to check whether the icing accretion and aerodynamic characteristic changes are the same under the same icing condition. An improved icing calculation method [24] based on the finite volume method (FVM) is used to carry out multistep ice numerical calculation of hybrid and full-scale airfoils.

Generally, ice can be divided into glaze ice and rime ice according to ice temperature and liquid water content. For rime ice, as the temperature is low, water droplets will freeze immediately when they hit the airfoil [25] and ice accretion gradually increases along the flow direction at the leading edge; the ice does not radically alter the shape of the airfoil; however, it does increase the chord length of the airfoil. So rime ice is a good test case to compare if the droplet trajectory impact locations and icing limits of hybrid and full-scale airfoil were the same [26]. Here, we choose a rime ice condition to check the ice accretion for hybrid and full-scale airfoil. The ice simulation parameters are shown in Table 4.

Figure 14 shows the icing on the surface of the hybrid airfoil and the full-scale airfoil at the end of several time steps. It can be seen that the ice shapes coincide at s1 and s5 during the beginning of ice accretion, with no obvious differences between them. In later time steps, at the end of s8 and s14 iterations, the ice shapes of the lower part still have excellent agreement, but the ice shape of the upper part has a small difference. However, the upper and lower limit positions, the ice thickness limit, and the ice angle thickness [27] are close and thus still have a familiar resemblance. Comparing the ice shape in different iteration steps in Figure 14, the chord length of the airfoil after icing gradually increases in the negative direction compared with the clean airfoil.

Figure 15 shows the Cp comparison between the hybrid airfoil and the full-scale airfoil at the end of each time step. The ice shapes at the end of s1, s5, and s8 are in good agreement, especially in the area near the front edge, and there is no obvious difference in Cp. At the end of s14, the ice shape difference is caused by the icing accumulation error. At this time, the Cp oscillates and the suction peak of the hybrid airfoil weakens to a greater extent than that of the full-scale airfoil. It can also be seen that although ice accretion leads to the separation of air flow, which affects the upper and lower pressure distribution of the leading edge, the influenced region is small. Most of the oscillations occur in the 10% chord length of the leading edge, and thereafter, the Cp gradually recovers to the same as that of the clean airfoil. Only at the end of s14, because of the complex shape of ice accretion, does the separation of air flow seriously affects the Cp of the whole airfoil which makes the Cp of the iced airfoil smaller than that of the clean airfoil.

**(a)**

**(b)**

**(c)**

**(d)**

##### 5.2. Ice-Induced Vortices

The shape of the hybrid airfoil is such that the chord length is short and the aft section has a large curvature. Therefore, with the ice accretion on the leading edge of airfoils, flow separation and ice-induced vortices occur and grow gradually. Then, droplets will be affected and change their trajectories rapidly in this area. In order to check the vortices in the ice accretion process of hybrid and full-scale airfoils, vortex distribution of two airfoils (S5 and S14 of hybrid and full-scale airfoils) was simulated by the 2D URANS equations. When calculation is complete, the instantaneous vortex structures of these airfoils were compared by the -criterion [28] of the vorticity equipotential surface, which is defined as follows: where and are the symmetric and antisymmetric components of the velocity gradient tensor.

Instantaneous vortex diagrams of full-scale and hybrid airfoils are shown in Figure 16. The left picture shows the vortices of airfoils as the ice accretion begins, S5, and the right shows the end of accretion. As we can see, the vortex distribution is complex and concentrated in the leading edge; for the full-scale and hybrid airfoil, the vortices share the same distribution location and strength. As for S14, with more ice on the leading edge, the irregular ice induced more complex vortices and separation. Compared to S5, the vortices have larger and stronger distribution, especially in the lower surface area. As shown in the circles, vortices are similar for full-scale and hybrid airfoils at the same stations. These results illustrate the same flow and vortices for full-scale and hybrid airfoils in the ice accretion process and then the same droplet trajectories and ice accretion.

#### 6. Conclusions

Hybrid airfoils based on different full-scale airfoils and different AoA are obtained; leading edge pressure measurement experiment was conducted. The design and experimental results demonstrate that, by the same aerodynamic characteristics and Cp in the reference area, hybrid airfoil with good performance can be designed using the optimization method in this paper.

Compared with the parameterized expression method for general airfoils, the SS function is simple and clear because of its special function and use, which can better describe the shape characteristics of the hybrid airfoil, accurately express the shape of the hybrid airfoil, and easily change the airfoil profile through parameter adjustment.

The multiobjective genetic algorithm and two optimization loops improved the optimization accuracy of multicontrol point hybrid airfoils. Taking the same aerodynamic characteristics of the reference area as the optimization criterion, more accurate simulation of pressure distribution of full-scale airfoil can be realized and the same stagnation point and pressure peak value can be obtained.

The pressure measurement experiment proved the excellent performance of hybrid airfoils designed by this method; the instantaneous vortex structure comparison indicated that the hybrid airfoil shares the same vorticity distribution and variation with the full-scale airfoil, so ice accretion on the hybrid airfoil is the same as that on the full-scale airfoil.

#### Data Availability

All the data supporting the results were shown in the paper and can be applied from the corresponding author.

#### Disclosure

Initial idea and technical route verification began in 2018; it was introduced in Journal of Physics: Conference Series, 1786: 012023 as this paper's start.

#### Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

This research was supported by the Open Fund of Key Laboratory of Icing and Anti/De-icing (IADL20190106) and the National Project of Equipment’s General Technology (41406020501).