**About the Editors**

#### **Spiros Pantelakis**

Spiros Pantelakis is Director of the Laboratory of Technology and Strength of Materials of the Mechanical Engineering and Aeronautics Department at the University of Patras. He served as Chairman of the Department of Mechanical Engineering and Aeronautics from September 2009 to September 2013. In September 2007 to August 2010, he was an Executive Board member of the Research Committee of the University of Patras. Furthermore, between September 2010 and August 2013, he was Vice Chairman of the Executive Board of the Research Committee of the University of Patras. He is one of the founding members and, since 2008, has been Chairman of the European Aeronautics Science Network (EASN) Association. From 2006 to 2011, he was Representative of the European Aeronautics Academia in the plenary of ACARE, and from 2005 to 2011, he was Chairman of ACARE's Working Group on Human Resources and Chairman of the Board of Directors of the Greek Metallurgical Society from 2008 to 2011. Since 2015, he has been a member of the Board of Directors of the Hellenic Aerospace Industry. He has over 35 years' experience in the field of aerostructures and aeronautics materials, during which he has been involved in over 100 international aeronautics research projects—in many of them, as Scientific Coordinator. He is a member of the Editorial Board and Guest Editor-in-Chief for a number of international scientific journals and Reviewer for many international scientific journals. Furthermore, he has been author or served as Editor in numerous books published by various international publishing houses. He is the author of more than 250 scientific publications in international peer-reviewed journals and conference proceedings. He has been Chairman of some dozens of international conferences and workshops, and he is founder of the Conference Series ICEAF (International Conference of Engineering Against Failure).

#### **Andreas Strohmayer**

Andreas Strohmayer is professor of aircraft design at the Institute of Aircraft Design, University of Stuttgart, with a research focus on manned (hybrid) electric flight ("Icare 2" and "e-Genius") and ´ scaled UAS flight testing. He studied Aeronautical Engineering at TU Munich and graduated in 2001 under Dr.-Ing. in the field of conceptual aircraft design. From 2002 to 2008, he was director of Grob Aerospace in Mindelheim, Germany, responsible for the design, production and support of the Grob fleet of all-composite aircraft: the development of a four-seat aerobatic turboprop, a seven-seat turboprop, and the SPn business jet. From 2009 to 2013, he was program director for a 19-seater commuter aircraft project at Sky Aircraft in Metz, France, then VP Programs at SST Flugtechnik in Memmingen, Germany, setting up an EASA-approved design organization, holding the TC for a six-seater all-composite turbo-prop aircraft. He then left the industry to join University of Stuttgart in 2015, teaching aircraft design and promoting electric flight. In 2016, he became member of the Board of Directors of the European Aeronautics Science Network (EASN), and since 2019, he has been an EASN Chairman. In this position, he also represents academia in the ACARE General Assembly.

### *Editorial* **Special Issue "9th EASN International Conference on Innovation in Aviation & Space"**

**Spiros Pantelakis 1,\* and Andreas Strohmayer 2,\***


This Special Issue contains selected papers from works presented at the 9th EASN International Conference on Innovation in Aviation & Space, which was successfully held in Athens, Greece, between the 3rd and 6th of September 2019. The event included 9 keynote lectures and more than 360 technical presentations distributed in approximately 70 sessions. Furthermore, 40 HORIZON2020 projects disseminated their latest research results, as well as the future trends on the respective technological field. In total, more than 450 participants joined the 9th EASN International Conference.

In the present Special Issue, eight engaging articles are contained, with more than 1800 views each until now, related to aviation and space research. Slavik et al. [1] performed a multi-objective optimization in order to select a fixed propeller for a short take-off and landing (STOL) category aircraft. The aim was to achieve the highest possible performance with fixed propeller, i.e., high maximal horizontal and cruise speed, short take-off and high rate of climb; Pareto sets were implemented in order to select the optimal propeller. As a result, a high-performance power system with a low price (fixed pitch propeller) for STOL category aircraft could be designed. Plagianakos et al. [2] studied the effect of hot–wet storage aging on the mechanical response of a carbon fiber polyether ether ketone (PEEK)-matrix woven composite. A wide range of static loads and selected cyclic load tests on the interlaminar fatigue strength were performed, with the results providing a useful basis towards preliminary design with PEEK-based woven thermoplastic composites during service in aerospace applications. Capovilla et al. [3] designed a CFRP structural/battery array configuration in order to integrate the electrical power system with a spacecraft bus primary structure. The results indicated that, by implementing the designed configuration, more volume and mass was made available for the payload, compared with traditional, functionally separated structures employing aluminum alloys. The study of Khustochka et al. [4] proposed a novel method for a stable estimation of the engine performance parameters using a priori information about the engine, its mathematical model and expected performance, in view of fuzzy sets, and also about the measuring system and measuring procedure. A comparison of the proposed approach with traditional methods underlined its main advantage regarding its high stability of estimation in the parametric uncertainty conditions, while the proposed method can be implemented for matching thermodynamic models to experimental data, gas path analysis, as well as adapting dynamic models to the needs of the engine control system. In the work of Piscitelli et al. [5], a superhydrophobic coating for metallic substrates with a simplified and non-expensive method was developed, which could be employed as a usual paint able to prevent/reduce the formation of ice, especially on small aircraft. The surface properties and the wettability of the developed coating were investigated, with the authors concluding that the specific coating can be potentially employed as a passive anti-icing system for aeronautical applications. Kellerman et al. [6] investigated the potential of using existing

**Citation:** Pantelakis, S.; Strohmayer, A. Special Issue "9th EASN International Conference on Innovation in Aviation & Space". *Aerospace* **2021**, *8*, 110. https:// doi.org/10.3390/aerospace8040110

Received: 12 April 2021 Accepted: 13 April 2021 Published: 14 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

aircraft surfaces as heat sinks for the waste heat of a (hybrid-) electric drive train. The results pointed out that surface heat exchangers can provide cooling power in the same order of magnitude as the waste heat expected from (hybrid-) electric drive trains for all sizes of considered aircraft. Dvirnyk et al. [7] assessed the serviceability limit of the blades of the axial compressor of a helicopter engine operating in a dusty environment, and showed that the gradual loss of the stall margin over time determines the serviceability limits of compressor blades. The serviceability limits defined by the authors could enable helicopter users to significantly reduce operating costs by extending the remaining useful life (RUL) of the engines operated in a desert environment. Finally, in the study of Biser et al. [8], the possibilities of coupled, analytical models for sizing electric propulsion systems were demonstrated by considering the example of a propulsion system developed in the frame of a relevant EU-funded project. The approach proposed by the authors was found to favor model accuracy and low computational effort. Potential technical solutions to decrease the influence of the DC power transmission system were also given.

The editors of this Special Issue would like to thank the authors for their high-quality contributions and for making this Special Issue manageable. Additionally, the editors would like to thank Ms. Linghua Ding and the Aerospace editorial team.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Propeller Selection by Means of Pareto-Optimal Sets Applied to Flight Performance**

#### **Svatomir Slavik, Jan Klesa \* and Jiri Brabec**

Department of Aerospace Engineering, Czech Technical University in Prague, Prague 121 35, Czech Republic; svatomir.slavik@fs.cvut.cz (S.S.); j.brabec@fs.cvut.cz (J.B.)

**\*** Correspondence: jan.klesa@fs.cvut.cz

Received: 30 November 2019; Accepted: 28 February 2020; Published: 5 March 2020

**Abstract:** Selection process of the propeller for short take-off and landing (STOL) category aircraft is described. The aim is to achieve the highest possible performance with fixed propeller, i.e., high maximal horizontal and cruise speed, short take-off and high rate of climb. These requirements are contradictory and so Pareto sets were used in order to find the optimal propeller. The method is applied to a family of geometrically similar propellers that are suitable for 73.5 kW (100 hp) piston engine designed for ultralight category aircraft with maximal take-off weight of 472.5 kg. The propellers have from two to eight blades, blade angle settings from 15◦ to 40◦ and diameter from 1.1 m to 2.65 m. Pareto frontier is designed for each pair of flight conditions, and the optimal propeller is selected according to these results. For comparison, the optimal propeller selection from the propeller family by means of a standard single-optimal process based on the speed power coefficient *cs* is also used. Use of Pareto sets leads to considerable performance increase for the set of contradictory requirements. Therefore, high performance for a low price for the given aircraft can be achieved. The described method can be used for propeller optimization in similar cases.

**Keywords:** STOL aircraft; propeller; Pareto sets; propeller optimization

#### **1. Introduction**

Optimal propeller performance has been investigated from the beginning of aviation. The first scientific work on propeller aerodynamic optimization was performed by Betz and Prandtl [1]. A more precise method was described by Goldstein [2]. The theory of Betz and Prandtl was used by Larrabee for the aerodynamic design of the propellers with low loading [3]. This method is quite popular until today due to its simplicity and good results. It was later developed, e.g., by Adkins and Liebeck [4] or Hepperle [5]. All these methods can be used for the aerodynamic design of optimal propeller for given flight condition.

This paper is focused on the choice of optimal propeller from given family, i.e., group of propellers having identical blade shape and various diameter and number of blades for various flight conditions, i.e., requirements for the propeller are usually contradictory. Standard method for the solution of this problem by means of power speed coefficient *cs* is described in [6–8]. This procedure is relatively simple and determines optimal diameter and blade angle for the group of propellers having the same blade shape and number of blades. This works only for a single operating point.

Therefore, if propeller performance should be maximal in several flight conditions, some multiple objective optimization method should be used. This brings difficulties in identifying the propeller with the best performance and thus the optimal solution. It is possible to use Pareto-optimal sets for this task. Its application is quite frequent. It can be used for the multi-objective optimization of ship propellers [9,10], optimization of UAS powerplant [11] or propeller optimization with manufacturing constraints [12]. Other possibilities are, e.g., gradient-based methods [13], genetic algorithms [14] or calculus of variations [15].

Multi-objective propeller optimization described in [9–15] is focused on the propeller design. This paper brings a new point of view. It combines the approach from [6–8] with multi-disciplinary optimization using Pareto-optimal sets. The optimal propeller is selected from the family of propellers, i.e., propellers having the same blade shape but different diameter and number of blades. The main advantage of this problem definition is the prevention of possible structural problems. Increase in the number of blades together with lower blade cord leads to higher efficiency, but it leads also to lower bending stiffness and strength (e.g., for the blade chord decreased to 50% of its baseline value, bending stiffness decreases to 12.5% and bending strength to 25% of its baseline values). To avoid this, optimal propeller is selected from the set similar to the approach in [7,8].

#### **2. Problem Formulation**

The task is to choose the optimal fixed-pitch propeller for a given aircraft and engine. The objectives are to maximize both cruise and maximal horizontal flight speed and at the same time minimize the take-off distance of the aircraft. These results are contradictory, so some sort of trade-off must be accepted. The method itself can be used for any propeller driven aircraft. Further, it is presented in the example of the light short take-off and landing (STOL) category aircraft. Development of this aircraft also was main motivation for the research in this branch. The aircraft was originally equipped with constant-speed propeller. The usage of the fixed-pitch propeller was motivated by the decrease in both system complexity and cost.

#### **3. Methods**

#### *3.1. Pareto Sets of Flight Performance for Selection of Optimal Propellers*

Multi-objective optimization by means of Pareto sets is based on search for boundary values of objective functions dependent on the same variable *X* of these functions. Values of objective functions are evaluated for every value of the independent variable *X* and plotted in the graph (see Figure 1). An element of Pareto set represents mutual relation of all considered objective functions for one given value of the independent variable *X*.

**Figure 1.** Scheme of Pareto set (redrawn according to [16]).

In the case of multi-objective optimization with two objective functions, Pareto sets have the form of plane graphs. One axis represents values of the first objective function and the second one values of the second objective function, as shown in Figure 1. Corresponding pairs of both functions (for the same variable values *X*) create Pareto set as the area of the graph. The set area limits are bound to a defined condition range of an independent variable for each function. The bound of Pareto sets makes outer limits of achievable values of the objective functions.

Peak points of the bound (Points 1, 2, 3, and 4 in Figure 1) depict extremes (minimum, maximum) of one and the other objective functions. Parts of the boundary between extremes match optimal values of variable *X*, because if the extremes of both these functions are simultaneously required

(e.g., maximum 1 for *u*2(*X*) and maximum 2 for *u*1(*X*)), the extreme is achievable only for the first or second function. The change on Pareto boundary when dropping under the extreme of the first function approximates in the optimal way the solution to the extreme of the second function (for each variable *X*, between 1 and 2, both functions reach their maximum values). All variables *X* on Pareto boundary between these two extremes (Pareto-optimal front) are optimal (the most appropriate) because it is not possible to decide which combination of both objective functions is better. None of these solutions is worse or better—the solutions are mutually non-dominant.

In the case that the extreme of one function is also the extreme of the second function (for one value of the variable *X*, both extremes are reached) then Pareto-optimal front passes at one point (e.g., extremes 1 and 2 are identified) and the multi-objective optimization is one optimal solution only.

The goal of optimization is to find a non-dominant solution that requires the functions to be in contradiction (conflict). A change of variable *X* for one function towards its required extreme value delays the second function from its required minimum/maximum. A more detailed description of multi-objective optimization by means of Pareto sets can be found in literature [16] or [17].

In the optimization selection of propellers from the family of propeller, the discrete points represent individual propellers. The propeller geometry comprises both continuously geometric parameters (mainly distribution of the chord length, twist and thickness of the airfoils used along the blade, the propeller diameter) and discrete parameter—the number of blades. The propeller aerodynamic characteristics represent the dependence of the thrust and power coefficients on the advance ratio (*cT*(λ), *cP*(λ)) that together with an engine power curve enables to set the available isolated thrust curve of the power unit.

Pareto sets are composed from final number of points (equal to the number of propellers in the family of propeller) of appropriate pairs of the contradictory flight conditions objective functions: [maximum horizontal flight speed–take-off distance], [maximum horizontal flight speed–maximum rate of climb], [take-off distance–maximum rate of climb]. The pair of [maximum horizontal flight speed (continuous engine regime)–maximum horizontal flight speed (cruise regime)] is also included.

The optimization selection of the propeller by means of Pareto-optimal sets consists of the following sequential actions:


#### *3.2. Aerodynamic Characteristics of the Aeroplane*

A model study of a small two-seat sport airplane with a requirement for a short take-off and landing (STOL) was chosen. The airplane is designed as the high-wing arrangement in which 13 m2 trapezoidal wing of the aspect ratio equals to 7.2. The wing is equipped with a combined high-lift device on the leading and trailing edge—a slotted leading edge (slat) and Fowler flap. The whole wing leading edge is equipped with slat; Fowler flap is installed on 60% of the wing trailing edge. Wing airfoil was developed from GA(W)-1 in order to increase lift coefficient. The airplane is drawn up with the standard side-by-side seat arrangement and a fixed taildragger landing gear with fairing.

Lift and drag aerodynamic characteristics were determined according to the methodology described in Appendix B for the following flight lift device configuration:


The primary characteristics are determined without the propeller influence on the aerodynamic drag of the airplane and without the ground effect (see Figure 2). The characteristics correspond to moment-balanced states in the typical flight configuration. The primary lift and drag characteristics of take-off and landing configurations are corrected for the ground effect depending on the actual height above the ground. An additional airplane aerodynamic drag due to the propeller stream is included in the true propeller thrust.

**Figure 2.** Polar graph of the aircraft for cruise (i.e., flaps retracted) and take-off configuration (i.e., flaps 15◦).

#### *3.3. Powerplant*

The ROTAX 912 Aircraft Engine has been selected. The engine, primarily designed for this category of airplane, has take-off power of 73.5 kW and is equipped with a propeller speed reduction unit with gear ratio *i* = 2.43. The take-off regime corresponds to a maximum power at maximum permissible revolution per minute (RPM) for a short-term use, maximum 5 minutes (i.e., 5800 RPM at engine shaft); the continuous regime is 90% of the take-off power at reduced engine RPM (i.e., 5500 RPM at engine shaft) without any time limit; and the cruise regime (economic engine operation) means 75% of the continuous regime.

#### *3.4. Propellers*

#### 3.4.1. Propeller Family

The propeller family is created as a group of geometrically similar propellers. Propeller blade geometry is derived from a three-blade propeller suitable for a 73.5 kW (100 hp) piston engine. The three-blade propeller was extended to other numbers of blades: two-blade, four-blade, five-blade, six-blade and eight-blade propeller with eleven pitch blade angles from 15◦ to 40◦ at 75% of radial distance. Geometric characteristics of the propeller are presented in Figure 3 (blade-chord distribution), Figure 4 (blade twist distribution) and Figure 5 (relative airfoil thickness). Clark Y airfoil is used. Range of propeller diameter of the total propeller family was from 1.1 m to 2.65 m, but each number of blades for a given pitch blade angle is only a part of this range. An extremely overloaded propeller that reaches the permissible engine RPM practically at the beginning of the take-off limits the smallest diameter. The maximum diameter represents a very lightly loaded propeller, and thus, the take-off is unacceptably extended, and the climbing rate is reduced.

**Figure 3.** Propeller blade chord distribution.

**Figure 4.** Propeller blade geometric twist distribution.

**Figure 5.** Propeller blade relative thickness distribution.

#### 3.4.2. Aerodynamic Characteristics of Propellers

Aerodynamic characteristics of propellers needed to determine thrust curves are dependencies of the thrust and power coefficients on the advance ratio, i.e., *cT*(λ) and *cP*(λ). The vortex blade theory of an isolated propeller is used. Free helix vortex surfaces with a constant pitch leaving the boundary vortex of each blade generate a field of the induced velocities that adjust the magnitude and direction of free flow along the propeller blade. The aerodynamic forces acting along the blade can be considered as two-dimensional airfoil characteristics with the incoming free flow corrected to the induced angle of attack. The calculations of the propeller aerodynamic characteristics were performed using the numerical model [18,19]. Input geometric data involve, except diameter *D* and number of blades *N*, also the distribution of the chord length, twist and thickness of the airfoils used along the blade. The calculation was done for every blade pitch setting and number of blades. More detailed description of the method can be found in Appendix A.

#### *3.5. Thrust Curves of the Propeller*

The thrust curve of the power unit has the meaning of the thrust available in dependence on flight speed. First, the thrust curves of the isolated propellers were calculated, and then, these isolated thrusts were corrected for influence of the installation. The corrected thrust represents the true thrust. Because the additional drag of the airplane (both friction and pressure part) due to the increased speed of the propeller stream is thus dependent on the propeller thrust, it is preferable for the calculations of the flight performance to introduce so-called effective thrust. The effective thrust is the true thrust reduced by the additional propeller drag. The drag curve of the airplane without the additional propeller drag thus remains the same for all alternatives of the propeller propulsion units. The friction component depending on the thrust and wetted area influenced by the propeller flow was evaluated according to literature [20].

#### *3.6. Flight Performance for Pareto Sets*

Aircraft performance is computed by methods described in [6] and [21]. More detailed explanation can be found in Appendix B. Computation is performed for 0 meters international standard atmosphere (ISA), i.e., air pressure 101,325 Pa, air density 1.225 kg·m−<sup>3</sup> and air temperature 288.15 K (15 ◦C). Pareto sets require the pairing of such requirements, which are contradictory. The respective pairs of flight conditions required to achieve their extreme values (minimum, maximum) act on each other so that increasing one flight condition towards the extreme decreases the second flight condition from its desired extreme.

In general, it can be expected that the propellers with good take-off performance are not good for reaching maximum flight speeds and vice versa. Therefore, for pairs of the flight conditions for Pareto set, the following are used:


#### **4. Results**

Pareto sets of six-selected flight condition pairs discussed in the Section 3.6 are plotted in Figures 6–11. Every point in the graph represents aircraft performance for one propeller defined by its number of blades, diameter and blade angle setting. Flight performance is computed according to the methodology described in Section 3. Marker types depend on the number of blades. For clarity, only limited areas near Pareto-optimal fronts from the all propeller family combinations are depicted. The broken lines connect points from Pareto-optimal fronts. Selection by means of speed power coefficient *cs* is used for comparison (see *cs* opt points in Figures 6–11). Performance with original constant-speed propeller is also shown for comparison. Results are further discussed in Section 5. Tables 1–6 contain coordinates of the points on the Pareto fronts for corresponding Pareto graphs. Pareto fronts are formed only by two-blade propellers; propeller diameter and blade angle are mentioned for every point. Performance with selected optimal propeller is also plotted in Figures 6–11).

**Table 1.** Points on the Pareto front for maximum horizontal speed (continuous engine regime) and take-off distance (see Figure 6).


0D[LPXPOHYHOVSHHG9+PD[FRQWLQXRXVUHJLPH>NPK@

**Figure 6.** Pareto frontier for maximum horizontal speed (continuous engine regime) and take-off distance.

0D[LPXPOHYHOVSHHG9+&FUXLVHUHJLPH>NPK@

**Figure 7.** Pareto frontier for maximum horizontal flight speed (cruise engine regime)–take-off distance.

0D[LPXPOHYHOVSHHG9+PD[FRQWLQXRXVUHJLPH>NPK@

**Figure 8.** Pareto frontier for maximum horizontal flight speed (continuous engine regime) and maximum rate of climb (take-off engine regime).

0D[LPXPOHYHOVSHHG9+&FUXLVHUHJLPH>NPK@

**Figure 9.** Pareto front for maximum horizontal flight speed, (cruise engine regime) and maximum rate of climb (take-off engine regime).

**Figure 10.** Pareto frontier for take-off distance (take-off engine regime) and maximum rate of climb

(take-off engine regime).

0D[LPXPOHYHOVSHHG9+PD[FRQWLQXRXVUHJLPH>NPK@

**Figure 11.** Pareto frontier for maximum horizontal flight speed (continuous engine regime) and maximum horizontal flight speed (cruise regime).

**Table 2.** Points on the Pareto front for maximum horizontal speed (continuous engine regime) and take-off length (see Figure 7).


**Table 3.** Points on the Pareto front for maximum horizontal speed (continuous engine regime) and maximum rate of climb (take-off engine regime) (see Figure 8).



**Table 4.** Points on the Pareto front for maximum horizontal speed (cruise engine regime) and maximum rate of climb (see Figure 9).

**Table 5.** Point on the Pareto front for take-off distance (take-off engine regime) and maximum rate of climb (take-off engine regime) (see Figure 10).


**Table 6.** Points on the Pareto front for maximum horizontal flight speed (continuous engine regime) and maximum horizontal flight speed (cruise regime) (see Figure 11).


Figure 6 shows Pareto set for the Maximum level speed (continuous engine regime) vs. take-off length. Pareto front is located in the lower right-hand corner (short take-off and maximum speed). Pareto front is created uniquely by two-blade propellers. Two-blade propeller optimized by the means of the power speed coefficient *cs* is on the Pareto-optimal set. Figure 7 shows similar situation for the cruise speed. In this case, Pareto-optimal set contains also uniquely two-blade propellers, but the solution obtained by the means of the speed power coefficient is no more on the Pareto-optimal set.

Figures 8 and 9 present Pareto-optimal front for the cases maximum level speed (continuous engine regime) vs. rate of climb and cruise speed vs. rate of climb. Both parameters should be maximal, i.e., optimal solutions are in the top right-hand part of the graph. Like before, optimal solution obtained by means of speed power coefficient is part of optimal set only for the case of the maximal continuous engine power.

Figure 10 presents Pareto frontier for the combination of take-off length vs. rate of climb. Optimal set is created only by single point in this case. Figure 11 presents Pareto frontier for the combination of maximum level speed (continuous regime) vs. cruise speed. Both parameters should be maximal, so the optimal set is located in the top right-hand part of the graph. As in previous cases, optimal set is represented only by two-blade propellers. Optimal solution obtained by means of the power speed coefficient is the part of this set.

#### **5. Selection of Optimal Propeller**

Optimal propeller is selected from the group of near to the Pareto fronts in the Figures 6–9 and Figure 11. The choice itself depends on the weight factors given to the different parameters, i.e., trade-off is always necessary. Take-off length was selected as the most important parameter due

to the STOL category aircraft. Weight factors can be of course modified and tuned according to the designer requirements. The selected optimal propeller has two blades, diameter *D* = 2.05 m and pitch blade angle φ0.75 = 20◦. This point is part of the Pareto fronts in Figures 6, 8 and 11. It is also close to the Pareto front in other cases. Aircraft performance with the optimal propeller is described in Table 7.

**Table 7.** Aircraft performance with chosen optimal propeller (diameter *D* = 2.05 m and pitch blade angle φ0.75 = 20◦).


#### **6. Conclusions**

Multi-objective optimization using Pareto sets of flight performance of 472.5 kg small STOL sport airplane powered by 73.5 kW (100 hp) engine to select a fixed propeller from a family of geometrical similar propellers is presented. The propeller family included propellers from two-blade to eight-blade with pitch blade angles from 15◦ to 40◦ and diameters range from 1.1 m to 2.65 m.

The optimization criteria required maximum level speed, maximum rate of climb and minimal take-off. The optimal Pareto fronts were investigated for four pairs of opposing flight conditions:


Length of take-off path and maximum rate of climb relate to the take-off engine power regime; maximum horizontal flight speed was considered for both continuous and cruise power regime. Only four pairs of parameters are used for the optimization due to the fact that the last two pairs (take-off distance vs. rate of climb and cruise speed vs. maximal horizontal speed) in fact do not affect the choice of the optimal propeller.

Pareto-optimal study leads to the following conclusions:


Acceptable aircraft performance is reached for the propellers with two blades, diameter in the range from *D* = 1.95 to 2.2 m with a corresponding decrease of blade angle from φ0.75 = 22.5◦ to 12.5◦. Optimal propeller is then selected from this group.

To compare Pareto multi-objective optimization with standard selection propeller from a family of geometrically similar propellers by means of the speed power coefficient, the design was performed for design speed from the middle of Pareto-optimal front of maximal horizontal speed for continuous engine regime. The single-mode optimization selection for two-blade propeller confirms Pareto optimization selection: diameter *D* = 2.05 m and pitch blade angle φ0.75 = 20◦. The speed power coefficient methods are limited by single-regime design defined by given flight speed and engine regime, and the design is not clearly connected to flight performance. Two-blade propeller of this family considered as a constant speed propeller with diameter obtained by the speed power coefficient presents an increase in flight performance for the fixed pitch propeller.

**Author Contributions:** Author contributions are divided in the following way: conceptualization, S.S. and J.B.; methodology, J.K.; software, S.S. and J.K.; validation, S.S. and J.B.; formal analysis, S.S.; writing—original draft preparation, S.S. and J.K.; writing—review and editing, J.K.; visualization, S.S.; supervision, J.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the EU Operational Programme Research, Development and Education, under the Centre of Advanced Aerospace Technologies, project No. CZ.02.1.01/0.0/0.0/16\_019/0000826, Faculty of Mechanical Engineering, Czech Technical University in Prague.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**



#### **Appendix A. Computation of Propeller Characteristics**

Method for the computation of propeller characteristics is based on the combination of blade element theory and vortex theory. Dimensionless variables are used for the computation of isolated propeller performance. Standard formulations are used (see [19]). Flight velocity is expressed by advance ratio λ

$$
\lambda = \frac{V}{n\_{\text{\\$}}D}.\tag{A1}
$$

Thrust is expressed by thrust coefficient *cT*

$$cr = \frac{T}{\varrho n\_s^2 D^4}.\tag{A2}$$

Power coefficient *cP* is defined in a following way

$$\mathcal{L}\_P = \frac{N}{\varrho n\_s^3 D^5}.\tag{A3}$$

Radial coordinate r is replaced by

$$
\vec{r} = \frac{r}{R'} \tag{A4}
$$

Dimensionless velocities are defined in a following way

$$
\overline{\mathcal{U}} = \frac{\mathcal{U}}{\overline{\Omega}\mathcal{R}'} \tag{A5}
$$

and dimensionless circulation Γ is used in following form

$$
\overline{\Gamma} = \frac{\Gamma}{4\pi\Omega R^2}.\tag{A6}
$$

Thrust and power coefficients can be determined from [19]

$$
\kappa\_T = \pi^3 \int\_{\overline{\mathcal{T}}\_0}^1 \overline{\Gamma} \left( \overline{\mathcal{U}}\_1 - \frac{c\_D}{c\_L} \overline{\mathcal{V}}\_1 \right) d\overline{\mathcal{T}}\_r \tag{A7}
$$

$$
\omega\_P = \pi^4 \int\_{\mathcal{I}\_0}^1 \overline{\Gamma} \left( \overline{V}\_1 + \frac{c\_D}{c\_L} \overline{\mathcal{U}}\_1 \right) \overline{r} d\overline{r}. \tag{A8}
$$

Velocities on the propeller blade is shown in the Figure A1. *U*<sup>1</sup> and *V*<sup>1</sup> are defined by

$$
\overline{\mathcal{U}}\_1 = \overline{\mathcal{U}}\_0 + \overline{\mathfrak{u}}\_\prime \tag{A9}
$$

$$
\overline{V}\_1 = \overline{V}\_0 + \overline{v}\_\prime \tag{A10}
$$

where *u* and *v* represent induced velocities (during normal propeller operation, *u* is negative, and *v* is positive). Vortex theory is used for the determination of induced velocities. System of horseshoe vortices is used (see Figure A2). Induced velocities are computed by the way described by Okulov [22].

**Figure A1.** Flowfield and forces on the propeller blade.

**Figure A2.** Illustration of the vortex system on the propeller blade.

Lift coefficient *cL* is coupled with circulation by

$$
\overline{\Gamma} = \frac{1}{2} \underline{c}\underline{\iota} N \overline{\tilde{c}} \overline{\mathcal{W}}\_1. \tag{A11}
$$

Airfoil section lift and drag coefficients are computed by means of the mathematical model described in [23]. Iterative procedure in MATLAB environment is used for the solution of the system of equations. The procedure was tested on the case of the three-blade propeller 5868-R6 described in [8]. Comparison of the computed results with experimental data is presented in Figures A3 and A4.

**Figure A3.** Comparison of computed values of thrust coefficient *cT* (solid lines) with experimental results from [8] (discrete points) for the three-blade propeller 5868-R6.

**Figure A4.** Comparison of computed values of power coefficient *cP* (solid lines) with experimental results from [8] (discrete points) for the three-blade propeller 5868-R6.

#### **Appendix B. Computation of Flight Performance**

#### *Appendix B.1. Thrust Curve of the Propeller Power Unit*

The thrust curve of the power unit has the meaning of the available thrust in dependence on flight speed. First, the thrust curves of the isolated propellers are calculated, and then, the isolated thrust is corrected for influence of the installation. The corrected thrust represents the true thrust.

Because the additional drag of the airplane (both friction and pressure part) due to the increased speed of the propeller stream depends on the propeller thrust, it is preferable for the calculations of the flight performance to introduce so-called effective thrust. The effective thrust is the true thrust reduced by the additional propeller drag. The drag curve of the airplane without the additional propeller drag thus remains the same for all alternatives of the propeller propulsion units.

#### Appendix B.1.1. Thrust of Isolated Fixed-Pitch and Ground-Adjustable Propellers

Determination of the thrust curve of an isolated propeller requires, as the first step, the solution of the equilibrium propeller–engine rotational speed:

$$P(n) = \rho(H)n\_s^3 D^5 c\_P(\lambda). \tag{A12}$$

The solution determinates the root of the power equation expressed by the engine power curve *P*(*n*) on one side and the propeller power described by the power coefficient *cP*(λ) for a given flight speed *V* on the other side. The equilibrium speed is solved for a number of speeds from the full range of the flight speed. The pair of the flight speed and the equilibrium rotational speed determinates the advance ratio λ and by help of the thrust coefficient *cT*(λ) the isolated thrust can be evaluated:

$$T\_{\rm is} = \rho(H) n\_s^2 D^4 c\_T(\lambda). \tag{A13}$$

The power engine curve *P*(*n*) pertains to one of the three engine regimes: Take-off, Continuous and Cruise. In the case that the equilibrium rotational speed with an increasing flight speed reaches the maximum allowed value *nmax*, it is assumed that the pilot will maintain this limit rotational speed by the throttle. The equilibrium rotational speed in such overload regime is thus the maximum speed *nmax*.

#### Appendix B.1.2. Thrust of Isolated Constant-Speed Propellers

As the constant-speed propeller keeps the rotational speed constant independently of the flight speed, then, the regulated rotational speed is actually the equilibrium speed. It is therefore possible to calculate the thrust directly by means of the thrust coefficient curve *cT*(λ) determined for the constant-speed propeller working in the required power regime.

#### Appendix B.1.3. Thrust Curve with Influence of Height

For calculating of thrust curves at different heights, the engine power of the piston engine at zero altitude can be converted to the given height *H* by a multiple factor *kN* according to [24]

$$k\_N(H) = \ (1 + 0.132) \frac{\rho(H)}{\rho\_0} - 0.132 \ . \tag{A14}$$

The air density ρ at height *H* [m] corresponds to the standard dependence according to ISA:

$$
\rho(H) = \rho\_0 \left( 1 - \frac{H}{44308} \right)^{4.256},\tag{A15}
$$

where

$$
\rho\_0 = 1.225 \,\mathrm{kg/m^3},
\tag{A16}
$$

and the constant 44,308 m represents the ratio of two individual constants: the constant temperature decreases with high 0.0065 K/m and temperature of 288 K for *H* = 0 m.

The influence of the height on aerodynamic characteristics of propellers can be neglected.

#### Appendix B.1.4. True Thrust of Installed Propellers

The ratio of the true propeller thrust *T* and the thrust of isolated propeller *Ti* is estimated as the ratio of the actuator disc thrust of the propeller with a mean flow speed through the disc and thrust of the isolated actuator disc. The mean flow speed corresponds to the deceleration downstream due to a body behind the propeller. The mean value of the decelerated flow through the disc is determined from the equality of the flow momentum by the actuator disc with a constant and variable velocity.

The graphical dependence published in [19] is used to describe the change of the speed in the plane of the propeller disc caused by an engine nacelle. The graph shows the relative velocity drop (Δ*V*1*P*/*V*1*P*) in the plane of the propeller disc in the area of the central part of the propeller. The central part is defined by the cross-section of the engine nacelle *Sn*. The velocity drop ratio depends on the cross-sectional area *Sn* and propeller disc area *SP* and a compensatory analytical form of the graph has the form:

$$\frac{\Delta V\_{1P}}{V\_{1P}} = 0.2 \frac{S\_{\text{fl}}}{Sp} - 0.08 \sqrt{\frac{S\_{\text{fl}}}{Sp}} + 0.028 \,\text{.}\tag{A17}$$

The ratio of the thrusts is determined from the momentum of the induced velocity of the ideal actuator disc in the steady-state distance behind the disc and momentum of the ideal actuator disc with the mean speed. The ratio determinates the final correction factor to convert the isolated thrust to the true thrust:

$$k\_{\rm is} = \left(1 - \frac{S\_n}{S\_P} \left(\frac{\Delta V\_{1P}}{V\_{1P}}\right)\right) \frac{2V\_{1P}\left(1 - \frac{S\_n}{S\_P}\left(\frac{\Delta V\_{1P}}{V\_{1P}}\right)\right) - V}{2V\_{1P} - V}.\tag{A18}$$

The flow velocity *V*1*<sup>P</sup>* through the propeller disc can be expressed depending on the thrust at flight speed *V* according to the theory of the ideal actuator disc:

$$V\_{1P} = \frac{V}{2} \left( 1 + \sqrt{1 + \frac{T\_{\rm is}}{\frac{1}{2}\rho(H)S\_{\rm P}V^2}} \right). \tag{A19}$$

In the case of the single-engine airplane, the nacelle cross-section *Sn* is replied by the cross-section area of the engine part of the fuselage.

#### Appendix B.1.5. Effective Thrust

The effective thrust means the true propeller thrust reduced by the additional airplane drag generated by the propeller stream consisting of a friction and pressure component:

$$T\_{eff} = \,^T - \left(\Delta D\_{fr} - \Delta D\_{pr}\right). \tag{A20}$$

The friction component Δ*Dfr* is depended on the thrust and a wetted area of an airplane influenced by the propeller flow. Its mean value is presented in References [20,24]:

$$
\Delta D\_{fr} = 0.004 \ T \frac{S\_{\text{net}}}{S\_P}. \tag{A21}
$$

The pressure component Δ*Dpr* is taken from Reference [19] for the engine nacelle case:

$$
\Delta D\_{\rm pr} = \frac{\Delta V\_{1\rm P}}{V\_{1\rm P}} T.\tag{A22}
$$

The conversion factor of the true trust to the effective thrust can thus be expressed in the form:

$$k\_{eff} = 1 - 0.004 \left(\frac{S\_{wet}}{S\_P}\right) - 0.2 \frac{S\_{\rm u}}{S\_P} - 0.08 \sqrt{\frac{S\_{\rm u}}{S\_P}} + 0.028 \,\text{.}\tag{A23}$$

*Appendix B.2. Curves of Available and Required Thrust and Power*

#### Appendix B.2.1. Available Thrust Curve

The available thrust represents the effective thrust *Te*ff*;* see (A20). The available thrust curve is the dependence of the effective thrust on the flight speed *V*.

#### Appendix B.2.2. Required Thrust Curve

The required thrust is equal to the thrust that corresponds to the aerodynamic drag force without additional drag caused by propeller in a balanced horizontal flight at a steady flight speed *V*.

Calculation procedure for determining the required thrust:

Step 1: To determine the needed lift coefficient from the balance of lift and weight for the chosen flight speed *V*:

$$c\_L = \frac{mg}{\frac{1}{2}\rho(H)V^2c\_LS} \,. \tag{A24}$$

Step 2: To state the corresponding drag coefficient from the aerodynamic polar of the respective flight configuration (without additional aerodynamic drag caused by propeller flow) for the corresponding lift coefficient. Drag coefficient *cD* is determined from aerodynamic polar (Figure 2). Step 3: The required thrust is equal to the aerodynamic drag:

$$T\_{req} = \frac{1}{2}\rho(H)V^2\mathfrak{c}\mathfrak{D}.\tag{A25}$$

Steps 1–3 shall be repeated from *Vmin* (minimum horizontal flight speed at the given configuration–stalling speed) to the maximum flight speed *VHC* (or *VHmax*) corresponding to the intersection of the required thrust curve *Treq*(*V*) with the available thrust curve *Te*ff (*V*) for a given altitude. An example of the curve of available and required thrust of a small sport aircraft at a cruise regime is depicted in Figure A5.

**Figure A5.** Curves of available and required thrust for a small sport aircraft *m* = 473 kg at a cruise configuration with three-blade propeller: *D* = 1.75 m, ϕ0.75 = 22◦, maximum cruise-regime power 47.9 kW.

Stalling speed:

$$V\_{\rm min} = \sqrt{\frac{mg}{\frac{1}{2}\rho(H)V^2c\_{L\rm max}S}}\,\,\,\,\,\tag{A26}$$

Maximum flight speed: see (A31).

#### Appendix B.2.3. Required Thrust Curve with Influence of Height

#### 1. Ground influence

Due to the ground proximity, aerodynamic drag decreases and lift increases. A simple correction according to Reference [21] is introduced for induced part of the drag and the lift-curve slope:

$$
\omega\_{\text{Dyson}} = \left. \mathcal{L}\_{\text{Drain}} + (1 - \sigma)(\mathcal{c}\_{\text{D}} - \mathcal{c}\_{\text{Dmin}}) \right|, \tag{A27}
$$

where the correction factor σ is dependent on the height *hW* of the wing above the ground, the aspect ratio of the wing λ*<sup>W</sup>* and its area *S*:

$$\sigma(h\_W, \lambda\_W, S) = \left[ \epsilon^{-4.22 \left( \frac{h\_W}{\sqrt{l\_W s}} \right)} \right]^{0.766} \,\tag{A28}$$

$$c\_{l,ground} = \frac{c\_{l,\pi\lambda\eta}}{\left[ (1-\sigma)c\_{l,\} } . \tag{A29}$$

2. Out of ground influence

Changes of aerodynamic forces are considered only by changing of the air density with height. The effect of Reynolds number changes on aerodynamic coefficients is neglected.

Appendix B.2.4. Available Power Curve

$$P\_{eff}(V) = T\_{eff}(V) \,\text{V}.\tag{A29}$$

Appendix B.2.5. Required Power Curve

$$P\_{req}(V) = \; ^r\_{req}(V) \; V. \tag{A30}$$

#### *Appendix B.3. Maximum Horizontal Speed, Maximum Rate of Climb*

Standard calculation models based on the curves of the required and available thrust and power are used to determine the maximum horizontal speed and maximum rate of climb; see, for example, Reference [20].

Maximum flight speed:

$$T\_{nq}(V) \ = T\_{eff}(V) \ \ \rightarrow \ V \ \ = V\_{HC} \tag{A31}$$

Flight speed of maximum rate of climb *Vvymax*:

$$(P\_{eff}(V) - P\_{req}(V)\_) = \max \to V = V\_{\mathcal{V}y\text{max}}\,. \tag{A32}$$

Maximum rate of climb:

$$w\_{\text{ymax}} = \frac{P\_{\it eff} \left(V\_{\text{voymax}}\right) - P\_{\it req} \left(V\_{\it voymax}\right)}{m\text{g}} \,. \tag{A33}$$

#### *Appendix B.4. Take-O*ff *Distance*

The take-off distance is understood to be the horizontal distance of the classic take-off procedure consisting of the ground round phase and three airborne phases: the ground flight, transition and climb out. The standard dynamic models of the take-off, discussed in, e.g., References [20,25], are adopted for the individual phases. The normal take-off considers the airplane with the take-off configuration of the lift device.

#### Appendix B.4.1. Take-Off Ground Run

The take–off ground run determines the distance to achieve the safe lift-off speed *V*1. Length of the ground run:

$$L\_1 = \int\_0^{V\_1} \frac{V}{R\_{\text{ground } run}(V)} dV \,\,,\tag{A34}$$

where safe lift-off speed *V*<sup>1</sup> is within the range (110–115)% of the stall speed with the ground effect. Stalling speed with the ground run effect corresponds to the relationship (A26):

$$\left| V\_{\min\\_\text{ground\\_run}} \right| = \sqrt{\frac{mg}{\frac{1}{2}\rho(H)V^2c\_{\text{Lmax\\_ground\\_run}}S}}\,\text{}\,\text{}\,\text{}\,\tag{A35}$$

$$\text{for} : \ c\_{L\text{max}\_{\text{ground}}\,\text{ran}} = \frac{c\_{L\text{max}}\pi\lambda\_W}{\left[\left(1 - \sigma\left(h\_{W\_{\text{ground}}\,\text{ran}}, \lambda\_W, S\right)\right)c\_L^a\right]}\,. \tag{A36}$$

*Rground run* represents an accelerating force during the ground run:

$$R\_{\text{ground\\_run}}(V) = T\_{eff}(V) - mgf - \left(c\_{D\_{\text{ground\\_no}}} - c\_{L\_{\text{ground\\_no}}}f\right) \frac{1}{2}\rho(H)V^2S \,. \tag{A37}$$

with such a take-off angle of attack that satisfies the optimal condition:

$$\frac{d\mathcal{L}\_{\mathcal{I}\_{\text{ground}}\,mw}}{d\mathcal{L}\_{\mathcal{I}\_{\text{ground}}\,mw}} = \frac{1}{f} \,. \tag{A38}$$

Tangent to the lift-drag at the ground run point is equal to the inverse value of the rolling friction coefficient 1/*f*.

#### Appendix B.4.2. Take-Off Ground Flight

The take–off ground flight corresponds to the acceleration phase of the ground level flight from the safe lift-off speed *V*<sup>1</sup> up to the safe take-off climb speed *V*2.

Length of the ground flight:

$$L\_2 = \int\_{V\_1}^{V\_2} \frac{V}{R\_{\text{ground } fllight}(V)} dV\_\prime \tag{A39}$$

where safe take-off climb speed *V*<sup>2</sup> is required at least 120% of the stall speed in the ground flight regime. Stalling speed with the ground flight effect corresponds to the relationship (A26):

$$W\_{\text{min\\_ground\\_flkght}} = \sqrt{\frac{mg}{\frac{1}{2}\rho(H)V^2\mathfrak{c}\_{\text{Lmax\\_ground\\_flkght}}\mathfrak{S}}}\,\text{}\,\text{}\,\tag{A40}$$

$$\text{for} : \quad \mathbf{c}\_{L\text{max}\_{\text{ground}}\gamma\lambda\eta\mathbf{a}} = \frac{\mathbf{c}\_{L\text{max}}\pi\lambda\_{\text{IV}}}{\left[\left(1 - \sigma\{\mathbf{h}\_{\text{IV}\_{\text{ground}}}, \lambda\_{\text{IV}}, \mathbf{S}\right)\right] \mathbf{c}\_{L}^{\mathbf{a}}\right]}.\tag{A41}$$

*Rground flight* represents an accelerating force during the ground flight:

$$\mathcal{R}\_{\text{ground}\,f\,\text{light}}(V) = \,^\*T\_{eff}(V) - \mathbf{c}\_{\text{D}\_{\text{ground}\,f\,\text{light}}}(V)\frac{1}{2}\rho(H)V^2\mathbf{S}.\tag{A42}$$

The aerodynamic coefficients correspond to the point on the lift-drag polar of the ground-run configuration (corrected to the ground effect) at which the lift coefficient is equal:

$$\mathcal{L}\_{I\_{\text{ground}}/I\_{\text{light}}}(V) = \frac{m\mathcal{g}}{\frac{1}{2}\rho(H)V^2\,^\circ S} \,. \tag{A43}$$

Both of these coefficients thus depend on the ground flight speed *V*.

#### Appendix B.4.3. Take-Off Transition

The take–off transition phase means a curved flight from the safe take-off climb speed to the transition straight climb. The transition is taken as a flight close to a vertical circular arc with a constant speed equal to the safe take-off climb speed *V*2.

Horizontal length of transition arc:

$$L\_3 = \begin{array}{c} r \sin(\theta) \end{array} \tag{A44}$$

where the radius of transition arc r deepens on the load factor *ny* caused by the curvilinear flight:

$$r = \frac{V\_2^2}{\left[\mathcal{g}\left(n\_\mathcal{Y} - \cos(\theta)\right)\right]}\tag{A45}$$

The maximum achievable value of the load factor at the constant speed *V*<sup>2</sup> is

$$m\_{\text{ymax}} = \left(\frac{V\_2}{V\_{\text{min\\_ground\\_full}}}\right)^2. \tag{A46}$$

This maximum value corresponds to the *cLmax ground flight* and leads to the shortest length of the take–off transition phase, but it is on the edge of the sharp flow separation. A smaller value *ny* should be considered, e.g., 80% *nymax*.

The output angle from the transition arc is given by the difference of the available and required thrust at the transition speed:

$$\Theta\_{\parallel} = \arcsin\left[\frac{T\_{eff}(V\_2) - T\_{\text{req}}(V\_2)}{mg}\right],\tag{A47}$$

where the required thrust equals:

$$T\_{pq}(V\_2) = \ c\_{D\_{ground}/l\eta\otimes t}(V\_2)\frac{1}{2}\rho(H)V\_2^2S. \tag{A48}$$

As the ground effect during the transition flight is gradually decreasing, an approximate procedure based on a mean value height *hW* of the wing above the ground is applied. The mean value can be estimated, e.g., as a mean value between the height of arc *harc* calculated with the ground effect of the take-off ground flight and without the ground effect (free flight).

$$h\_{W\_{\text{reaction}}} = \frac{1}{2} \left[ h\_{\text{dref}} \left( \sigma \left( h\_{W\_{\text{ground}} \parallel \text{light}} \right) + h\_{\text{dref}} \sigma \left( h\_W \to \infty \right) \right) \right] \tag{A49}$$

where the arc height is equal:

$$h\_{\rm attr} = \left[ r \left[ 1 - \cos \left( T\_{eff} \left( V\_2 \right) - T\_{\rm req} \left( V\_2 \right) \right) \right] \right] \tag{A50}$$

#### Appendix B.4.4. Climb Out

The climb out is the phase of the steady climbing flight at the speed of the transition phase from the transition arc until the obstacle height *hobs* (usually 15.25 m = 50 ft)) is reached:

$$L\_4 = \frac{h\_{\text{obs}} - \left(h\_{\text{W}\_{\text{gmax}}f\text{light}} + h\_{\text{avc}}\right)}{t\chi(\theta)},\tag{A51}$$

where the arc height *harc* and take-off output angle θ are determined for *hW\_transition*.

#### Appendix B.4.5. Total Take-Off Distance

Total take-off distance is equal to the sum:

$$L\_{\text{table}\tooff} = L\_1 + L\_2 + L\_3 + L\_4. \tag{A52}$$

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **E**ff**ect of Hot-Wet Storage Aging on Mechanical Response of a Woven Thermoplastic Composite**

**Theofanis S. Plagianakos 1,**†**, \*, Kirsa Muñoz 2,**†**, Diego Saenz-Castillo 3,**†**, Maria Mora Mendias 3,**†**, Miguel Jiménez 2,**† **and Vasileios Prentzias 1,**†


Received: 27 November 2019; Accepted: 20 February 2020; Published: 24 February 2020

**Abstract:** The effect of hot-wet storage aging on the mechanical response of a carbon fiber polyether ether ketone (PEEK)-matrix woven composite has been studied. A wide range of static loads and selected cyclic load tests on the interlaminar fatigue strength were performed. Static tests were conducted in batch mode, including on- and off-axis tension, compression, flexure, interlaminar shear strength (ILSS) and fracture tests in Modes I, II and I/II. Respective mechanical properties have been determined, indicating a degrading effect of aging on strength-related properties. The measured response in general, as well as the variance quantified by batch-mode test execution, indicated the appropriateness of the applied standards on the material under consideration, especially in the case of fracture tests. The material properties presented in the current work may provide a useful basis towards preliminary design with PEEK-based woven thermoplastic composites during service in aerospace applications.

**Keywords:** carbon fibre thermoplastic composite; PEEK matrix; woven; aging; mechanical testing; static and fatigue

#### **1. Introduction**

Continuous-fiber reinforced thermoplastic composites are gaining attention in the aerospace industry for exhibiting advantages compared to thermoset composites, such as design and manufacturing flexibility, including multiple post-forming processes, and capability of being processed by a large range of traditional machining methods, fast cycle time and recyclability. As far as their mechanical performance is concerned, their enhanced impact resistance is a very attractive feature for selecting them in demanding lightweight applications. Moreover, in woven ply configurations they yield less anisotropic mechanical properties, which could be desirable in the context of conceptual design. In this context, the experimental determination of their mechanical properties over a wide range of static and fatigue mechanical tests, as well as quantification of the effect of aging on these properties, are extremely important for design and in-service monitoring purposes.

The main advantages of thermoplastic composites in comparison with thermoset composites have been very well described in the scientific literature to date. Excellent mechanical properties, good behaviour against impact, no need of cold storage owing to their long shelf-life and no chemical reaction during consolidation, which opens the possibility of short-time processing are among their

most appealing features [1–4]. There is a significant amount of literature focusing on mechanical characterization of continuous-fiber thermoplastic composites. Chu et al. [5] characterized 3-D braided Graphite/polyether ether ketone (PEEK) composites by static tensile tests and experimentally determined time-dependent mechanical properties at various temperatures by performing relaxation experiments and dynamic mechanical analysis (DMA) tests. Hufenbach et al. [6] studied the strain-rate dependency of the mechanical properties of three thermoplastic composite materials, including impact energy absorption, at different temperatures, and additionally studied the effect of fiber–matrix interphase modification on transverse tensile strength. Liu et al. [7] developed a damage model accounting for fiber failure and matrix cracking and validated it with an open-hole compression test on woven PEEK specimens. Kuo et al. [8] performed 4-point bending tests on thermoplastic composites and studied the effect of molding temperature on flexure properties and failure modes. Boccardi et al. [9] used infrared tomography to study the effect of frequency on the temperature of cantilever glass- and jute-based woven thermoplastic composite specimens subjected to cyclic bending. Sorentino et al. [10] fabricated and characterized polyethylene naphthalate (PEN) thermoplastic composites at 100 ◦C by 3-point bending, DMA and impact tests. As far as shear loading characterization is concerned, Hufenbach et al. [11] performed Iosipescu tests on woven thermoplastic composites and used high-speed camera and digital image correlation for studying the deformation and failure in order to develop material modeling strategies in commercial finite element software. Zenasni and coworkers [12,13] experimentally studied the effect of fiber material and weave pattern on the fracture response of polyetherimide (PEI)-matrix-based woven thermoplastic composite specimens in Mode I and Mode II.

The effect of aging on static and dynamic response of engineering thermoplastics has been studied since the early 1990s. As a general conclusion, PEEK polymers (and, therefore, PEEK-reinforced composites) are not the thermoplastic polymers which may absorb the highest level of water. It has been reported by Buchman and Isayev [14] that PEEK composites may absorb approximately 0.2% of humidity, in comparison with other thermoplastic polymers, which may absorb above 1%. This is mainly caused by the semi-crystalline structure of PEEK, in comparison with other thermoplastic polymers which have a larger amorphous part. Béland [15] concluded that semicrystalline thermoplastics/carbon composites may absorb also around 0.2% of humidity while carbon/amorphous thermoplastic may absorb beyond 0.8%. The effect of thermal aging on carbon-fibre reinforced PEEK has been studied by Buggy and Carew [16], who performed static and fatigue flexure tests. Aging has been applied by two different methods: by storage for up to 76 weeks at high temperatures and by thermal cycling from room temperature (RT) to 250◦C. They indicated a low sensitivity of the laminate flexural modulus and strength to ageing at 120 ◦C. Kim et al. [17] developed an analytical model for predicting the flexural properties degradation at high temperatures (540–640 ◦C) and performed 4-point bending tests for validation purposes. In the context of an implant application study, Schambron et al. [18] experimentally determined the effect of environmental ageing on static and cyclic flexure response of carbon-fibre/PEEK coupons and reported superior fatigue resistance compared to stainless steel. An experimental comparison of the bearing strength of woven-ply-reinforced thermoplastic or thermoset laminates at 120 ◦C after hygrothermal aging has been performed by Vieille et al. [19].

Thus, it can be concluded that the effect of aging on the mechanical response of woven thermoplastic composites has yet not been quantified over a wide range of mechanical tests. The present work aims to close this void by presenting an extensive test campaign performed to assess the mechanical properties of a high-performance woven carbon-fiber reinforced thermoplastic material in non-aged and aged conditions. Material characterization has been achieved by conducting mechanical tests according to (static tests) or based on (fatigue tests) ASTM and ISO standards, such as tension, compression, in-plane and interlaminar shear, flexure, Mode I, II and I/II fracture. Properties derived from static tests are reported in batch mode, to provide a measure of the test-type non-linearity and testing procedure repeatability. Moreover, the effect of aging is assessed by measuring mechanical properties after specimen environmental conditioning in hot-wet storage. Finally, albeit the main objective of

the current work is to provide experimental data, in order to check the applicability of linear elastic fracture mechanics based on effective material properties, interlaminar fracture test cases in Mode I and II have been modelled in a commercial finite element software and the predictions were compared with measured values. The main conclusion of the present work is that aging leads to significant degradation of the strength in engineering woven thermoplastic composites, while stiffness-related properties seem to be rather insensitive to aging.

#### **2. Materials and Methods**

#### *2.1. Specimen Manufacturing*

Thermoplastic composite coupons for mechanical properties' characterization were cut from a set of carbon fiber/PEEK plates manufactured by FIDAMC (Getafe, Spain). The engineering polymer selected for this study consisted of Tenax®TPCL PEEK-4-40-HTA40 3K supplied by Toho-Tenax (Tokyo, Japan). The material was carbon-fiber-reinforced PEEK (CF/PEEK) fabric 5HS (5 harness satin), 0.31 mm nominal thickness per ply, 40% of resin weight fraction and a fiber areal weight (FAW) of 285 gsm. The stacking sequence of the thermoplastic laminates varied depending on the mechanical testing and defined by the standards, explained in detail in the next sub-sections. It should be noted that, due to the fact that the used CF/PEEK fabric was a 5HS (non-symmetrical), the 0◦ direction of each ply was defined as the warp direction of the roll material, as suggested by the material supplier.

On the basis of previous experience with other candidate manufacturing methods [20], the selected manufacturing process for the consolidation of the laminates was compression molding by hot-platen press. The utilized equipment is located at the FIDAMC facilities (Getafe, Spain). Each laminate was first-hand laid-up with the help of a manual welder and then located in a metallic frame which acted as material retainer. Two polyimide sheets with a release agent were placed at both sides of the laminate, acting as release films. Two metallic caul plates were also used in order to obtain a proper flat-surface finish. The consolidation cycle consisted of a heating ramp at approximately 2 ◦C/min up to a consolidation dwell of 30 min at 400 ± 10 ◦C with an applied pressure of 1 MPa.

#### *2.2. Hot-Wet Storage Aging*

During its service life, an aircraft is exposed to high temperatures and high levels of humidity. The properties of composite materials may be affected as a consequence of moisture absorption and high temperatures. A faithful replication of the environmental exposure during aircraft operation would require a cyclic conditioning procedure between hot/humid and cold/dry conditions, as dictated by relevant standards (MIL-STD-810 or other). In the proposed work, we have followed an accelerated conditioning procedure on the basis of common practice [21], which would form a small part of an extended experimental aging campaign towards material airworthiness certification.

In order to evaluate the degradation of mechanical properties, an environmentally conditioned testing scenario was considered according to ASTM D5229/D5229M. The selected method is a recommended pre-test conditioning method, consistent with the recommendations of the Composite Materials Handbook-17 [22]. Specifically, the procedure that has been followed in these test series is a conditioning procedure, BHEP, which covers non-ambient moisture conditioning of material coupons from the same batch as the ones tested, widely known as traveler specimens, in a humidity chamber (-H-) at a prescribed constant, conditioning the environment to equilibrium (-E-), periodic (-P-) coupon weighing being required. Two different types of test series were covered. On the one hand, room temperature (RT) tests were conducted on specimens without any previous conditioning. Room Temperature conditions are controlled to meet standard laboratory atmosphere conditions, which according to ASTM D5229/D5229M, are 23 ± 2 ◦C and 50% ± 10% relative humidity. On the other hand, parallel test series were performed on specimens that have been previously conditioned. These are known as Hot Wet (HW) tests.

The conditioning parameters are usually fixed according to the conditions to which aircraft structure may be subjected during its service life. This commonly means an equilibrium moisture weight in an 85% relative humidity environment and a temperature of 70 ◦C. Nevertheless, in order to achieve a lower completion time for the testing campaign, accelerated conditioning was carried out. In order to achieve a faster aging process, conditioning was performed at a higher temperature (80 ± 3 ◦C), under the same relative humidity (85% ± 5%). It was checked that glass transition temperatures for the assessed material, as provided by the manufacturer, were significantly higher than the accelerated conditioning temperature. By proceeding this way, a decrease in conditioning periods was achieved, speeding the rate of testing.

#### *2.3. Experimental Methods*

All mechanical tests presented in this work have been performed according to (static tests) or based on (fatigue tests) relevant ASTM or ISO standards. Each standard has been reported in the title of relevant the subsection of the following section. Static tests have been performed in a five specimen batch mode to provide a statistically valid representation of a material sample response under static loads. For the sake of cost savings in the case of such an extended test campaign, three specimens per test type have been tested in fatigue, each one at a load level representing a fatigue cycle regime: low cycle-(1.0 <sup>×</sup> 104 cycles), medium cycle-(1.0 <sup>×</sup> 10<sup>5</sup> cycles) and high cycle-fatigue (1.0 <sup>×</sup> 10<sup>6</sup> cycles). Each test was assumed to have concluded with a failure if either the specimen failed due to rupture or if a load-displacement slope decrease of at least 10% was detected compared to the slope at the 100th loading cycle. A stress ratio of *R* = σmin/σmax = 0.1 has been considered, which is typical for the characterization of carbon-fiber reinforced plastics (CFRP) under dynamic loading [23,24]. Regarding testing frequency, high frequencies may produce an increase in the specimen temperature, which leads to mechanical properties degradation [25]. To avoid this, a frequency of 5 Hz was applied in all fatigue tests and the temperature of the specimens was monitored by means of a thermocouple placed on each. The laboratories involved in the campaign considered a range of specimen temperatures from RT up to 35 ◦C to be acceptable. For the selected frequency of testing, no specimen temperature reached 35 ◦C at any of the tests performed, and thus no effect of frequency on the results was assumed. According to the testing strategy, the fatigue endurance limit was set to 1.0 <sup>×</sup> 106 cycles.

All tests of the non-aged composite system (RT) have been conducted in Hellenic Aerospace Industry (HAI) facilities on an INSTRON 8801 (Norwood, MA, USA) hydraulic testing machine at room temperature (*T* = 25 ± 1 ◦C, 45% ± 5% humidity), equipped by an INSTRON load cell with a range up to 100 kN. In the case of fracture tests, a 5 kN Omega load cell has been used. In the case of double-cantilever beam (DCB) and mixed-mode bending (MMB) tests, a Philips SPC2050NC digital camera and Debut v4.08 by NCH software (Greenwood Village, CO, USA—non-commercial use edition) have been used to record crack propagation, whereas, in end-notch flexure (ENF) tests, crack propagation has been visually monitored. The camera-to-specimen distance was ca. 130 mm and focus has automatically been applied at a 640 × 80 pixel resolution and a rate of 30 frames per second.

HW tests have been performed on different Universal Testing Machines in Element facilities. Static tension (0◦, ±45◦, 90◦) and compression (0◦, 90◦) tests were performed on a Zwick Z100 BS1, (Kennesaw, GA, USA) equipped with a Zwick load cell up to 100 kN. Static flexure, interlaminar shear strength (ILSS) and interlaminar fracture toughness tests were performed on a Zwick Retroline (INSTRON 5866) Universal Testing Machine, equipped with an INSTRON load cell having a range up to 10 kN. All static HW tests have been conducted under controlled temperature (70 ◦C) using thermocouples inside a temperature chamber (Thermcraft). Regarding fatigue, HW tests, tension (0◦, ±45◦ laminates) and ILSS tests were performed on a Universal Dynamic Testing Machine INSTRON 8872, equipped with an INSTRON load cell up to 25 kN. For fatigue tension (90◦ laminate), a Universal Dynamic Testing Machine INSTRON 8801, equipped with an INSTRON load cell up to 100 kN, was used. In all fracture tests crack propagation has been visually monitored.

HW fatigue tests were performed at RT right after being environmentally conditioned. That choice was based on the assumption that there is an irreversible degradation of the material due to specimen conditioning. As reported in the literature [26], long-term environmental conditioning leads to irreversible changes that cause permanent property alterations within the matrix, the fiber surfaces and the fiber/matrix interface. Therefore, by testing unaged and aged specimens under fatigue loading at RT conditions, it is possible to evaluate the effect of permanent degradation caused by hygrothermal conditioning.

#### *2.4. Finite Element Models*

The numerical simulations for pure fracture Modes I and II have been performed using MSC MARC finite element (FE) software [27]. In both cases, cohesive elements have been integrated into the models in order to evaluate the prediction of delamination propagation with experimental results. The composite thermoplastic material is modelled using property values at ply level extracted from mechanical tests. The constitutive relation of the cohesive elements is based on Linear Elastic Fracture Mechanics (LEFM) and an exponential damage law is used [28]. The cohesive properties of the interface have been derived by experimentally determined values for normal traction (strength—Section 3.3.1), interlaminar shear traction (strength—Section 3.5.1) and critical energy release rates in Mode I (Section 3.5.3) and Mode II (Section 3.5.4). From these values, normal traction and fracture toughness in Mode I have been adapted according to the strategy formulated in [29] in order to yield a critical opening displacement allowing for solution convergence. As discussed in Sections 3.5.3 and 3.5.4, the modeling approach adopted herein, which has been successfully applied in unidirectional thermoset materials [28], failed to accurately estimate the maximum load in both fracture modes of the woven thermoplastic material.

In order to simplify the analysis of both DCB and ENF tests, 2D models have been created using plane strain full-integration elements (Type 11) for the bulk material and cohesive elements (Type 186) for the interface. There are four elements through the thickness for the DCB model and five elements for the ENF model because of the difference in specimen thickness. An element length of 0.75 mm has been selected for composite and interface in order to achieve an acceptable convergence rate.

#### **3. Results and Discussion**

The test campaign for material characterization in terms of mechanical properties includes static and fatigue tests. Static tests have been performed for the determination of basic mechanical properties to feed the structural design phase of thermoplastic parts, as well as for quantification of the degradation of these properties due to aging during service. Fatigue tests provide an overview of fatigue endurance of the material via measured S–N curves and validate the trend observed in the static tests concerning the effect of aging. Specifically, static tests include tension, compression, in-plane shear, flexure, interlaminar shear, and interlaminar fracture toughness in Mode I, II and I/II. Fatigue tests include interlaminar shear strength. The specimen nominal dimensions in fatigue were identical to the ones used in the corresponding static case. Wherever available, the measured data are compared with values provided by the manufacturer [30].

#### *3.1. Specimen Dimensions and Lamination*

The specimen dimensions have been selected in order to comply with the relevant standard. The lamination has been determined on the basis of relevant standards, while available experimental capabilities have also been taken into account, in order to remain within the measurable load range of available equipment. In Table 1, the geometry and lamination of the specimens are presented. The test type description is provided as T-0 (tension warp), T-90 (tension weft), C-0 (compression warp), C-90 (compression weft), FLEX (flexure), T-45 (tension ±45), ILSS (interlaminar shear strength), ILFT I (interlaminar fracture toughness Mode I), ILFT II (interlaminar fracture toughness Mode II), ILFT I/II (interlaminar fracture toughness Mixed-Mode).



#### *3.2. Specimen Aging*

For each test type, the respective moisture absorption curve, measured from corresponding traveler specimens (in-plane dimensions 25 × 25 mm), is reported. Figure 1 provides an overview of the conditioning process. On the basis of these curves, the moisture content uptake rate and water content at saturation during conditioning have been extracted in Table 2. For the determination of the slope, a measurement at saturation and the initial measurement have been used.

**Figure 1.** Conditioning curves for test types considered; the exposure time is expressed in square root of days.

**Table 2.** Moisture uptake parameters extracted from conditioning curves.


In order to comment on the trends observed in the conditioning process, the thickness of the relevant traveler specimen should be taken into account, since in-plane dimensions have been identical for all traveler specimens. It may be concluded that thicker traveler specimens, such as those considered for in-plane and interlaminar shear tests (including fracture toughness) yield generally higher moisture uptake rates and maximum values. This trend is attributed to the fact that moisture penetrates more easily across the sides of the specimens along thickness, which are uncoated compared to the upper and lower faces. Moreover, it is highly related to the effect of aging on the response in each test, as will be explained in the following sections.

Concerning comparison with data reported in the open literature, the maximum value of moisture uptake obtained for T-0 is very close to data reported by other researchers [31] (0.13%) and within the range set by previous studies on the moisture absorption of engineering thermoplastics [14,15].

#### *3.3. Determination of Fiber-Dominated Properties*

#### 3.3.1. Static Tension of Woven Laminate Warp Direction (ASTM D3039)

The head velocity was 2 mm/min, as dictated by the standard. Strains were measured by strain gages: biaxial rosettes have been attached on each specimen, whereas an additional longitudinal strain gage has been attached on the opposite side of one specimen to ensure that bending remained below 10% throughout the test. Data for stress vs. longitudinal strain for RT and HW specimens are presented in Figure 2. It should be noted that the maximum load may differ from the value shown, as strain gages occasionally got detached at a lower load level. Moreover, in the legend of Figure 2 the custom naming convention pattern followed in all subsequent data figures is provided.

**Figure 2.** Static tension of warp-directional laminate: stress vs strain data and characteristic failure mode at RT. The description of tests follows a custom naming convention: S—static, T—tension, 0—warp direction, RT—room temperature test, HW—hot wet test, last digit—specimen number.

Elaboration of the above static tension test data yields the mechanical properties shown in Table 3. Properties are provided in terms of mean value (MV) and coefficient of variance (CV). CV is directly related to the standard deviation (STDEV) as: CV=STDEV/MV, and thus STDEV could be alternatively used for reporting the deviation intervals in the experimentally determined values.

The strain range between 1000 and 3000 με has been used for extraction of the tensile elastic modulus. The respective values provided in the manufacturer datasheet according to ISO 527-4 Type 3 are also listed for comparison.

**Table 3.** Effect of hot-wet storage aging on mechanical properties derived from static tension tests on warp-directional laminates.


It may be observed that the elastic modulus is insensitive to aging (2% difference), whereas strength is strongly affected in a negative sense (21% reduction). This trend agrees well with results obtained by Solvay for PEEK thermoplastic polymer prepregs [32], and may be at least partially explained by the fact that the modulus is evaluated in the linear strain range (1000–3000 μm/m), while strength involves larger strains. It is reasonable to assume that in such a high demanding situation in terms of strain, the load-bearing capacity of the material appears to be sensitive to aging. As far as the Poisson ratio is concerned, it also seems to be affected, however, its measure is small and thus not expected to be of major importance as a design variable.

A typical failure mode is also illustrated in Figure 2. Lateral failure modes at multiple areas and various locations (LMV according to ASTM D3039) have been primarily observed. This trend has been expected, since the woven pattern of the fabric in the laminate prevents the occurrence of exploding modes. Similar failure patterns have been observed in RT and HW specimens, indicating that failure mode is rather insensitive to aging. Most failure patterns (8 out of 10) included failure near the tabs, indicating the redistribution of stresses along the specimen near/at failure load. The subsequent failure mode is related, partially at least, to inertia forces following the first failure, to a 3-D stress state involving local stress fields near the tabs, and to free-edge effects. The present approach focuses on the determination of the effect of hot-wet storage aging on the global mechanical properties of the woven thermoplastic material, while a detailed study of the failure mechanisms involved [11] would require special equipment, which has not been available in the context of this work.

#### 3.3.2. Static Tension of Woven Laminate Weft Direction (ASTM D3039)

The selection of a different lamination compared to the warp-directional specimens, as listed in Table 1, is explained by the fact that the tests shown in this work are part of a larger test campaign where unidirectional laminates of various material systems are mechanically evaluated. Wherever possible, identical specimen dimensions for all material systems have been selected with the purpose of maximizing comparability for the same mechanical property evaluation. In this context, the tests on the weft-directional specimens have been similar to the warp-directional ones in terms of head velocity and strain gauge placement.

**Figure 3.** Static tension for weft-directional laminate: stress–strain data and characteristic failure mode.

Elaboration of the above static tension test data yields the mechanical properties shown in Table 4. The same analysis was applied as in the case of the warp-directional static tension tests. The respective values provided in the manufacturer datasheet according to ISO 527-4 type 3 are also listed for comparison.

**Table 4.** Effect of hot-wet storage aging on mechanical properties derived from static tension tests on weft-directional laminates.


Concerning the effect of aging, a similar trend is observed, as in the case of the warp-directional specimens. The elastic modulus appears to be slightly (9%) increased due to aging, while strength is slightly (7%) reduced. Such a trend for the modulus has been expected, as both laminations are similar, whereas the reduction in strength is less than expected.

A typical failure mode of an RT specimen is also illustrated in Figure 3. Lateral failure modes have been observed in all specimens, being less distributed along the specimen compared to the case of warp-directional laminates. Most specimens (7/10) failed in the vicinity of the tabs. These failures are attributed to pre-stress, due to excessive pressure on the tabs at the grips of the machine in order to avoid slippage during loading.

3.3.3. Static Compression of Woven Laminate Warp Direction (ASTM D695)

Two batches in total, each consisting of 5 RT and 4 HW specimens, were tested. The head velocity was 1.3 mm/min. Strains were measured by two unidirectional strain gages, each attached on the narrow part of each side of the dogbone specimen to ensure that bending remained below 10% up to buckling.

A typical load vs. longitudinal strain measurement is presented in Figure 4. A typical failure mode is also shown. All specimens failed under intralaminar shear.

**Figure 4.** Compression test for warp-directional laminate: typical measured load vs longitudinal strain data and typical failure mode (SG: strain gauge).

Elaboration of the above static compression test data in the strain range between 1000 and 3000 με yields the elastic modulus shown in Table 5. The respective values provided in the manufacturer datasheet according to EN 2850 are also listed for comparison.


**Table 5.** Effect of hot-wet storage aging on compressive properties of warp-directional laminates.

It may be observed that compressive modulus and strength exhibit low sensitivity (less than 10%) to aging. This trend agrees with the results reported in [32] for a similar material. The values measured for the modulus are close to those provided by the manufacturer, as they do not depend on the specimen geometry accounted for in the two standards. On the contrary, the higher thickness in most of the specimens due to the consideration of tabs in EN 2850 leads to higher strength values compared to the ones measured in the current study.

#### 3.3.4. Static Compression of Woven Laminate Weft Direction (ASTM D695)

Two batches in total, each consisting of 5 RT and 4 HW specimens, were tested. The head velocity was the same as in the warp-directional case and strains were also measured accordingly. A typical load vs. longitudinal strain data curve is presented in Figure 5 along with a typical failure mode. All specimens failed under intralaminar shear.

**Figure 5.** Compression test for weft-directional laminate: typical measured load vs. longitudinal strain data and typical failure mode.

Elaboration of the above static compression test data in the strain range between 1000 and 3000 με yields the elastic modulus shown in Table 6.

**Table 6.** Effect of hot-wet storage aging on compressive properties of [90] laminates.


As expected, the trends observed in the weft-directional lamination are similar to those reported in the warp-directional one (Table 5). Compressive modulus and strength exhibit low sensitivity (less than 10%) to hot-wet aging. The moduli measured are close to the value provided by the manufacturer, whereas there is a significant deviation in strength due to the application of different standards encompassing different specimen geometries.

#### 3.3.5. Static Flexure (ASTM D7264)

Four-point bending tests have been performed in order to determine static flexural properties. Support span was 128 mm in order to comply with the span to thickness ratio suggested by the ASTM standard. The deflection was measured using an extensometer in a configuration similar to an LVDT. The head velocity was 1.5 mm/min.

The load vs. mid-span deflection for RT and HW specimens is presented in Figure 6, along with a characteristic failure mode.

**Figure 6.** Static flexure: load-deflection data and typical failure mode.

Elaboration of the above test data yields the mechanical properties shown in Table 7. Respective values provided for a thinner lamination in the manufacturer datasheet according to EN 2562-A are also listed for comparison.

**Table 7.** Effect of hot-wet storage aging on the mechanical properties derived from static flexure tests.


The elastic modulus appears to be slightly increased due to aging (10% increase), as in the case of tensile tests. On the other hand, aging leads to a drop of 16% in flexural strength and 20% in failure strain. Experimental studies reported in the literature [16] showed that thermal aging at 120 ◦C did not affect flexural modulus and strength. Therefore, the results presented herein indicate that it is the wet storage that mainly leads to strength degradation in a PEEK-matrix based composite configuration.

As far as failure modes of the specimens are concerned, all specimens failed under tension at the bottom surface, either at (5/10, failure description TAB according to ASTM D7264) or between loading noses (5/10, failure description TBB). According to the standard, these failure modes indicate valid flexural strength. This point is further supported by the acceptable variance in experimentally determined values for strength.

#### *3.4. Determination of In-Plane Matrix-Dominated Properties*

#### Static Tension of ±45 Woven Laminate (ASTM D3518)

In order to measure in-plane shear properties of the woven composite material, static tension tests on ±45 laminates have been performed. The head velocity was 2 mm/min, whereas strains were measured in the same manner as the other tensile tests presented above.

The stress vs. shear strain data for RT and HW specimens are presented in Figure 7, respectively. Each strain gage failed at approximately 15,000 με. The achieved range allowed for the evaluation of the 0.2% offset (yield) strength. Nevertheless, the load continued to rise after the failure of the strain gages and eventually reached a plateau at about 23 kN. Thus, the ultimate strength has been evaluated from the maximum load. At that stage, the plastic deformation of the specimens was visible with the naked eye, albeit no rupture of the specimen has occurred. That trend indicates a ductile response of the material, similar to that of metals, excluding failure in terms of breakage.

**Figure 7.** Effect of aging on stress–strain curves of in-plane shear strength (IPSS) test.

Elaboration of the above static tension test data yields the mechanical properties shown in Table 8. Respective values provided in the manufacturer datasheet (measured by using a proprietary testing method) are also listed for comparison.

**Table 8.** Effect of hot-wet storage aging on the mechanical properties derived from static tension tests on [±45] laminates.


As was also observed in the behavior of the tensile modulus in the previous tensile cases, the shear modulus, which is evaluated in a range of 2000 to 6000 μm/m, appears to be insignificantly affected by aging (3% increase). On the other hand, the tangential modulus beyond linear regime, as well as offset strength and maximum stress, are significantly lower in the case of the aged material (10% and 16%, respectively). Thus, the experimental results indicate that aging is more effective in larger strains occuring in the non-linear part of the response, which is prominent in that type of test. Considering that the maximum moisture absorption of in-plane shear specimens is slightly higher than in the tensile specimens, although not much higher, a reasonable phenomenological explanation of this trend is attributed to the thickness to width ratio of the specimens. In the in-plane shear case, the thickness to width ratio is approximately three times higher than in the tensile cases, and thus there is a wider distribution of moisture along the laminate plane, as moisture penetrates across the sides along thickness, which are not protected by a coating layer. As the strains increase, this trend gets more effective on the response and is enhanced by edge effects, which grow dominant. The experimental verification of that phenomenological explanation would be interesting, however, no appropriate equipment was available in the context of the present work to study the micromechanics involved during loading.

#### *3.5. Determination of Out of Plane Matrix-Dominated Properties*

#### 3.5.1. Static Interlaminar Shear Strength Tests (ILSS) (ASTM D2344)

In the case of ILSS tests the support span was 24 mm, thus yielding a thickness aspect ratio of 4. The deflection was measured using the head displacement of the testing machine. The head velocity was 1 mm/min.

The load vs. mid-span deflection for RT and HW specimens is presented in Figure 8. A typical failure mode is also shown.

**Figure 8.** Effect of aging on load-deflection curves measured in ILSS test.

Elaboration of the above test data yields the mechanical properties shown in Table 9. For this test type, no data are included in the manufacturer datasheet.

**Table 9.** Effect of hot-wet storage aging on the mechanical properties derived from static ILSS tests.


The interlaminar shear strength appears to be the most sensitive to aging among the properties studied so far, as a reduction of 25% has been experimentally determined. Moreover, the load-displacement slope is also significantly affected within both the linear and non-linear part of the response. This trend has been expected, since the ILSS test is a matrix-dominated test type. In addition, taking into account the dominance of the free-edge effects in such a low thickness aspect ratio (length/thickness) specimens, the ILSS experimental data led credibility to the phenomenological explanation provided in Section 3.4 concerning the effect of thickness to width ratio on the response. The thickness to width ratio of the ILSS specimens is 2.5 and 30 times higher than that of IPSS and tensile specimens, respectively, and thus they are more prone to a wider distribution of moisture in the specimen.

As far as the failure mode is concerned, the specimens failed under inelastic deformation and interlaminar shear. Each of these failure modes are considered acceptable in the standard, however, in the current case the failure appeared to include both modes as a result of the ductility of the matrix compared to thermoset composites.

#### 3.5.2. Fatigue Interlaminar Shear Strength (ILSS) Tests

The effect of aging on interlaminar shear fatigue strength is presented in Figure 9. Aging has a clear degrading effect on the fatigue strength of the thermoplastic material, which is obvious in low-, medium- and high-cycle regimes. This trend is in-line with the experimental results of corresponding static tests and supports the argumentation based on the effectiveness of the thickness to width ratio on the enhanced sensitivity of ILSS specimens to hot-wet storage aging. It should be noted that all failures in fatigue have been attributed to a reduction in slope in the load-displacement loop beyond 10%. It should be also indicated that, as these tests are limited in quantity, statistical means cannot be drawn.

**Figure 9.** Effect of aging on fatigue interlaminar shear strength. Arrows to the right indicate that no failure has occurred.

#### 3.5.3. Static Interlaminar Fracture Toughness Tests Mode I (ISO 15024)

According to the standard, an initial delamination of 60 mm has been built in by the intervention of a release film at the middle of the stacking sequence in the fabrication process. A head velocity of 1 mm/min has been applied. An initial loading has been considered for creating a crack propagation of 3–5 mm to eliminate fiber bridging effects, followed by a reloading stage for further crack propagation. In Figure 10a, measured load vs. opening displacement data for RT and HW specimens at the reloading stage are presented. A predicted instance of crack propagation is shown in Figure 10b. An interesting observation during the test was the "stick-slip" response during crack propagation, which resulted

in the usability of few data points for the evaluation of fracture toughness. While in the case of unidirectional (UD) thermoset materials, a stable propagation has been measured [28], it being common practice and, as such, considered in the related standards, in the case of the woven thermoplastic material propagation occurred in maximum of four steps. As shown in Figure 10c, both the upper and lower face have been excessively deformed at each propagation step. Thus, strain energy was stored in the face up to a sudden release. At that point, the crack propagated to an arbitrary distance, indicating a highly non-linear response and leading to a relatively large variance in fracture toughness. This stick-slip response has not been observed in similar tests on PEI matrix-based woven thermoplastics reported in the open literature [13], whereas the values obtained for fracture toughness are in a similar range. Thus, the stick–slip response could be attributed to the type of matrix, adhesion between fibers and matrix and weave type.

**Figure 10.** Mode I fracture test: (**a**) measured and predicted load vs. opening displacement data; (**b**) instance of crack propagation prediction; (**c**) typical deformed shape during test.

The FE model, which is based on LEFM, succeeds at predicting the slope of the load-displacement curves, whereas deviation is observed for maximum load. The deviation observed indicates that detailed unit-cell models and the application of appropriate failure criteria for woven thermoplastic materials should be applied in order to improve this prediction.

Elaboration of the measured data on the basis of Corrected Beam Theory (accounting for large-displacement correction) yields the interlaminar fracture toughness values presented as mean values in Table 10. For this test type no data are included in the manufacturer datasheet.


**Table 10.** Effect of hot-wet storage aging on interlaminar fracture toughness in Mode I.

The aged specimens exhibited higher fracture toughness than the RT ones, which is not compatible with the trends observed in all other matrix-dominated test types, including the subsequent Mode II fracture tests. It could be partially related to plasticization effects at the crack front and the deformation mechanism of the faces. However, since the propagation evolves in limited steps, the contribution of each mechanism (strain stored in the faces within a geometrically non-linear state and purely fracture energy) can hardly be identified. For instance, the HW specimens may endure a larger deformation of the faces due to being more ductile as a result of conditioning process, which leads to the determination of higher GIc values.

#### 3.5.4. Static Interlaminar Fracture Toughness Tests Mode II (ASTM D7905)

According to the standard, an initial delamination of 45 mm has been considered by applying a 0.012 mm release film at the middle of the lamination. A head velocity of 0.5 mm/min was applied. A non-precrack (NPC) and a subsequent precrack (PC) loading stage have been applied, each consisting of two tests for compliance calibration and an additional one for crack propagation. In Figure 11, measured load vs. transverse displacement data for RT and HW specimens in the crack propagation test of the reloading stage are presented, whereas measured and predicted instances of crack propagation are also shown. In contrast to the observations made in the DCB test, the woven thermoplastic material exhibits a stable response, similar to the one observed in the case of thermoset UD materials [28]. The FE model succeeds in accurately predicting the slope of the measured curve, whilst it fails to adequately estimate the measured maximum load for crack propagation. This shortcoming of the numerical model is largely attributed to the lack of implementation of appropriate failure criteria for woven thermoplastic materials encompassing interlaminar shear effects.

**(a) Figure 11.** *Cont*.

**Figure 11.** Mode II fracture test: (**a**) measured and predicted load vs. transverse displacement data; (**b**) typical deformed shape during test; (**c**) instance of crack propagation prediction.

Elaboration of the measured data yields the interlaminar fracture toughness values for RT and HW conditions, respectively, presented in Table 11. For this test type, no data are included in the manufacturer datasheet.

**Table 11.** Effect of hot-wet storage aging on interlaminar fracture toughness in Mode II.


The aged specimens exhibited lower fracture toughness than the RT ones. Since interlaminar fracture toughness is a matrix-dominated mechanical property, that trend is in agreement with the behavior observed in other matrix-dominated tests, such as the IPSS and ILSS tests presented above. Nevertheless, the thickness to width ratio of the specimens involved is in the range of the IPSS ones, which is an additional argument for justifying the higher sensitivity to aging compared to tensile and flexure tests.

3.5.5. Static Interlaminar Fracture Toughness Tests Mode I/II (ASTM D6671)

A mode mixture of *R* = GII/(GI + GII) = 0.2 has been selected for the RT specimens and a corresponding lever length of *c* = 93 mm. An initial delamination of 45 mm has been considered by applying a 0.012 mm release film at the middle of the lamination. A head velocity of 0.5 mm/min was applied. In Figure 12, measured load vs. opening displacement data for RT specimens are presented, along with instances of measured crack propagation. The "stick-slip" response during crack propagation has been also observed in this case, leading to relatively few data points and large variances. In the case of HW tests, the specimen failed at the root of the upper face at a maximum opening displacement between 3 and 10 mm. Thus, various mode mixtures (*R* = 0.25, 0.3 and 0.5) have been considered, albeit without success in eliminating that trend.

**Figure 12.** Mode II fracture test: (**a**) measured and predicted load vs opening displacement data RT specimens; (**b**) typical deformed shape during RT test; (**c**) failure mode at root of upper face during HW test.

Elaboration of the measured data yields the interlaminar fracture toughness values presented in Table 12. For this test type no data are included in the manufacturer datasheet.


**Table 12.** Effect of hot-wet storage aging on interlaminar fracture toughness in Mode I/II.

The large variance observed indicates the highly non-linear nature of this test on woven thermoplastic PEEK materials. The aged specimens exhibited significantly lower fracture toughness than the RT ones. However, premature failure of the upper face, as well as the fact that the crack evolves in limited steps, did not allow for a proper estimation of fracture toughness in the case of the HW specimens. The latter two points might indicate that the standard should be further elaborated in order to prescribe a stable and uniform crack growth in woven thermoplastic laminates, especially under hot-wet aged conditions.

#### **4. Conclusions**

On the basis of the presented study, the following major conclusions may be drawn:


**Author Contributions:** Contribution of each co-author to the current work may be distinguished as: Conceptualization, all; Investigation, K.M. and T.S.P.; Methodology, M.J., K.M., T.S.P., and V.P.; Resources, D.S.-C. and M.M.M.; Validation, M.J., K.M., T.S.P., and V.P.; Visualization, T.S.P.; Writing—original draft preparation, all; Writing—review and editing, T.S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** The current work has received funding from EU Horizon 2020 Clean Sky II project SHERLOC (Structural Health Monitoring, Manufacturing and Repair Technologies for Life Management of Composite Fuselage) under Grant Agreement No. CS2-AIR-GAM-2014-2015-01.

**Acknowledgments:** The authors from HAI would like to thank our ex-colleague Stavros Kalogeropoulos for his major assistance with the experimental work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
