**Estimation of Performance Parameters of Turbine Engine Components Using Experimental Data in Parametric Uncertainty Conditions**

#### **Olexandr Khustochka 1, Sergiy Yepifanov 2, Roman Zelenskyi <sup>2</sup> and Radoslaw Przysowa 3,4,\***


Received: 30 November 2019; Accepted: 10 January 2020; Published: 16 January 2020

**Abstract:** Zero-dimensional models based on the description of the thermo-gas-dynamic process are widely used in the design of engines and their control and diagnostic systems. The models are subjected to an identification procedure to bring their outputs as close as possible to experimental data and assess engine health. This paper aims to improve the stability of engine model identification when the number of measured parameters is small, and their measurement error is not negligible. The proposed method for the estimation of engine components' parameters, based on multi-criteria identification, provides stable estimations and their confidence intervals within known measurement errors. A priori information about the engine, its parameters and performance is used directly in the regularized identification procedure. The mathematical basis for this approach is the fuzzy sets theory. Synthesis of objective functions and subsequent scalar convolutions of these functions are used to estimate gas-path components' parameters. A comparison with traditional methods showed that the main advantage of the proposed approach is the high stability of estimation in the parametric uncertainty conditions. Regularization reduces scattering, excludes incorrect solutions that do not correspond to a priori assumptions and also helps to implement the gas path analysis with a limited number of measured parameters. The method can be used for matching thermodynamic models to experimental data, gas path analysis and adapting dynamic models to the needs of the engine control system.

**Keywords:** gas turbine engine; performance model; gas path analysis; robust estimation; identification; regularization; fuzzy set; membership function

#### **1. Introduction**

Aircraft gas turbine engines are the most expensive and vital components of the aircraft. Their maintenance and operating costs depend on the engine performance, which is mainly determined by a condition of the engine gas path components such as the compressor, turbine, combustor, ducts, etc. For this reason, engine maintenance is usually supported by a prognostics and health management system (PHM), which usually implements a diagnostic procedure known as gas path analysis (GPA).

Nowadays, different approaches are used for processing flight and test-cell data to estimate component performance and detect faults [1]. These methods can be divided into two groups [2]:


Model-based and hybrid systems [5] are effective only if a validated engine model is available.

The non-linear thermo-gas-dynamic models of turbine engines [6] and identification procedures have been applied in diagnostics for more than 40 years. Identification adjusts the model to make its output parameters as close as possible to the experimental data [7–14]. Besides the significant improvement in the gas path simulation, the estimated parameters contain information about the health of each component. Model identification involves the minimization of a functional which is a sum of squared deviations, calculated as differences of calculated and measured values. This is the Least Squares Method (LSM) which is usually implemented in non-linear GPA solvers using iterative Newton–Raphson algorithm.

Matching the engine model to experimental data is characterized by the presence of multiple parameters, which can be used for the model correction. These parameters can be strongly correlated. At the same time, a quantity of measured parameters is strongly limited. These reasons decrease the stability of the correction procedure, which is based on LSM, and forced researchers to look for methods to improve the stability. Hence, matching engine models needs regularization.

For more stable solutions, V. Borovick and Ye. Taran applied the least modulus method [15], which is more robust to outliers. S. Yepifanov implemented the Levenberg–Marquardt method [16] as well as Singular Value Decomposition and ε-structuration [17,18]. A. Volponi et al. applied the Kalman filter [19], and this approach is followed and developed in many papers [20–24], whose authors improved stability of the algorithm and its applicability to a non-linear engine model. X. Chang et al. applied an alternative method based on the non-linear filtration (sliding mode observer) [25–27].

Stability of identification procedures can be improved by regularization. This tool has been long known and used in statistical parameter estimation but was not common in system identification until recently [28]. T. Breikin et al. presented a regularization-based approach to the estimation of engine model coefficients [29]. It was based on the regularity of the frequency response pattern over the operation range of the engine, but that method is not universal. This paper is based on the Levenberg–Marquardt algorithm and Tikhonov regularization which were also used by S. Guseynov et al. [30,31]. These papers did not analyse the bias nor consider how the regularization coefficient influences the bias and the noise, which could impede practical applications.

It is well-known that regularization reduces the modification of the LSM functional by adding the regularizing component to it. This causes biased estimates of the model coefficients. B. Roth, D. Doel et al. demonstrated this clearly [32]. The biased estimates cause errors in GPA diagnosing solutions. The errors of the estimated parameters caused by the bias need to be monitored, but the traditional methods of regularization essentially prevent this. There is no universal recommendation on a choice of the weight coefficient for the added regularizing component.

This paper presents a new method for the regularized identification of gas-turbine models based on a priori information (Figure 1). First, the traditional Tikhonov regularisation is applied to simulated engine data to analyse the impact of regularization on a non-linear least squares solution and select the optimal value of the regularisation coefficient (Section 2.3). Next, membership function is used to implement a priori information into non-linear identification procedure, which is based on a generic algorithm (Section 2.4). Finally, the proposed method is tested with a turboshaft model and test-cell data and compared with traditional approaches.

**Figure 1.** Regularized multi-criteria identification of gas-turbine model.

The proposed approach is a significant modification of the standard model-based procedure for estimation of the performance parameters of engine components. The ability to implement prior knowledge in different forms to stabilize the solution is the main contribution of the paper. Regularized multi-criteria identification enables engine users to implement Non-Linear GPA [1,2] in ill-conditioned configurations with a limited number of measured parameters.

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

#### *2.1. Basic Identification Procedure*

The model calculates parameters of the engine gas path <sup>→</sup> Y at steady-state modes depending on an operating point, external conditions <sup>→</sup> <sup>U</sup> and component performance <sup>→</sup> Θ. Hence, in the general case, it is presented as:

$$
\overrightarrow{\mathbf{Y}} = \mathbf{F}(\overrightarrow{\mathbf{U}}, \overrightarrow{\Theta})\tag{1}
$$

The linear model may be formulated as

$$
\vec{\delta Y} = \vec{\text{H}} \vec{\delta \Theta} \tag{2}
$$

which relates small deviations of the gas path parameters δ → Y and parameters of components' performances δ → Θ at a single operating point. H is an influence coefficient matrix (ICM). The linear model (2) is a component of the identification algorithm for the non-linear model.

Identification of the non-linear model (1) by measuring results <sup>→</sup> Y ∗ and <sup>→</sup> U ∗ requires determining estimations <sup>→</sup><sup>ˆ</sup> Θ, which are solutions of the following optimization task:

$$\Phi(\stackrel{\rightarrow}{\Theta}) = \|\stackrel{\rightarrow}{Y}^\* - \stackrel{\rightarrow}{Y}(\stackrel{\rightarrow}{\bullet}^\*, \stackrel{\rightarrow}{\Theta})\|, \quad \Phi(\stackrel{\rightarrow}{\Theta}) = \min \Phi(\stackrel{\rightarrow}{\Theta}) \tag{3}$$

*Aerospace* **2020**, *7*, 6

where is the Euclidean norm of the vector, Φ–least squares functional to be minimized.

The precision of estimation will improve if the amount of a priori information is higher than the number of unknown parameters. Therefore, the identification can be provided by a set of measurements that are done at N different operating points. For this purpose, the generalized vector of residuals is minimized:

$$
\overrightarrow{Z}(\overrightarrow{\Theta}) = \begin{bmatrix}
\overrightarrow{\mathbf{Y}}\_1^\* - \mathbf{Y}(\overrightarrow{\mathbf{U}}\_1^\*, \overrightarrow{\Theta}) \\
\overset{\rightarrow}{\to} \overrightarrow{\mathbf{Y}}\_2 - \mathbf{Y}(\overrightarrow{\mathbf{U}}\_2, \overrightarrow{\Theta}) \\
\vdots \\
\overset{\rightarrow}{\to} \overleftarrow{\mathbf{Y}}\_N - \mathbf{Y}(\overrightarrow{\mathbf{U}}\_N, \overrightarrow{\Theta})
\end{bmatrix} \tag{4}
$$

Model (1) which is included in (4) is numerical. Therefore, the minimization of residuals (4) is a numerical iterative procedure, which in each iteration determines the solution as a sum of the previous solution and current correction:

$$
\stackrel{\rightarrow}{\Theta}^{n+1} = \stackrel{\rightarrow}{\Theta}^n + \stackrel{\rightarrow}{\Lambda}^{n+1} \tag{5}
$$

and iterations continue until corrections become negligible.

The correction Δ → Θ n+1 is determined as a solution of the overdetermined system of linear algebraic equations:

$$\stackrel{\rightarrow}{Z}(\stackrel{\rightarrow}{\Theta}^{\text{n}})\stackrel{\rightarrow}{\Delta\Theta}^{\text{n}+1} = Z(\stackrel{\rightarrow}{\Theta}^{\text{n}})\tag{6}$$

where

$$\mathbf{C}(\stackrel{\rightarrow}{\Theta}^{\mathbf{n}}) = \begin{bmatrix} \mathbf{H}\_1(\stackrel{\rightarrow}{\Theta}^{\mathbf{n}}) \\ \stackrel{\rightarrow}{\mathbf{H}\_2(\stackrel{\rightarrow}{\Theta}^{\mathbf{n}})} \\ \vdots \\ \stackrel{\rightarrow}{\mathbf{H}\_N(\stackrel{\rightarrow}{\Theta}^{\mathbf{n}})} \end{bmatrix} \tag{7}$$

is a generalized matrix of influence coefficients, which is composed of elementary matrixes Hi( → Θ n ) determined for each operating point. In accordance with Equation (6), for each iteration, such correction Δ → Θ n+1 →

of required parameters is found, which minimizes residuals <sup>→</sup> Z( Θ) (4).

The known solution of the linear system (6) by the least squares method takes the form:

$$
\stackrel{\rightharpoonup^{n+1}}{\Delta\Theta}^{n+1} = \text{A}^{-1}\text{C}^{\text{T}}\text{G}\stackrel{\rightharpoonup}{\text{Z}}\tag{8}
$$

where A = CTGC–Fisher information matrix, G–weight diagonal matrix, whose elements are inverse to variations of measuring errors σ<sup>2</sup> y i.

Unfortunately, the least squares method is sensitive to outliers in the right part of the system (6), which can be related to faults in experimental data. Therefore, practical applications of the classic LSM suffer from the instability of estimations Δ → Θ n+1 and poor convergence of the identification algorithm as a whole. In some cases, the estimations <sup>→</sup><sup>ˆ</sup> Θ are far from the expected values. The reasons for such results may be the lack of empirical information, an excessive number of estimated parameters or correlation between two or more state parameters. This causes ill-conditioning of the Fisher Matrix and excessive deviations of estimations. These estimations lose a physical sense (are out of range of possible components' performances parameters variation—for example, efficiency is more than one) and can cause the model calculation program to crash.

Hence, the LSM identification procedure needs modification to provide its stability and physically adequate estimations.

#### *2.2. Regularized Identification Procedure*

In the known regularization procedure named for Tikhonov, the generalized functional is minimized, which besides the residuals of measured parameters includes a norm of the finding parameter vector with a weighting factor α. The identification task takes the following form:

$$\Phi'(\stackrel{\textstyle \rightarrow}{\Theta}) = \sum\_{i=1}^{m} \lg\_i \sum\_{j=1}^{N} \left[ \mathbf{y}\_i(\stackrel{\textstyle \rightarrow}{\mathbf{U}}\_{\mathbf{\hat{\jmath}}} \stackrel{\textstyle \rightarrow}{\Theta}) - \mathbf{y}\_{\mathbf{\hat{\jmath}}}' \right]^2 + \alpha \sum\_{k=1}^{r} \theta\_{\mathbf{k}'}^2 \cdot \Phi'(\stackrel{\textstyle \rightarrow}{\Theta}) = \min \Phi'(\stackrel{\textstyle \rightarrow}{\Theta}) \tag{9}$$

where gi is the weight coefficient that describes measurement error (**gi** = <sup>1</sup> σ<sup>2</sup> **i** ), yi is a component of vector <sup>→</sup> **<sup>Y</sup>** and <sup>θ</sup>**<sup>k</sup>** is a component of vector <sup>→</sup> Θ. T

The regularized identification is based on the least squares method described above. As it is shown in Equation (9), the minimized functional is extended by elements which are related to parameters → Θ to be found. Thus, the main changes in the identification procedure are the generalized vector of residuals <sup>→</sup> Z( → Θ) and generalized matrix of influence coefficients C( → Θ). The modified terms have the following structure:

$$
\stackrel{\rightarrow}{Z}'(\stackrel{\rightarrow}{\Theta}) = \left[ \stackrel{\rightarrow}{\underset{\alpha\Theta}{\rightarrow}} \stackrel{\rightarrow}{\Theta} \right] , \text{C}'(\stackrel{\rightarrow}{\Theta}) = \left[ \stackrel{\rightarrow}{\text{C}(\stackrel{\rightarrow}{\Theta})} \right] \tag{10}
$$

V

where I is an identity matrix.

This identification procedure will deliver the non-regularized solution if the regularization coefficient is zero. Otherwise, the influence of coefficient α depends on the proportion between components: m i=1 gi N j=1 yi ( → Uj, → Θ) − y<sup>∗</sup> ij2 and <sup>r</sup> k=1 θ2 <sup>k</sup> of the generalized functional Φ (Θ). Therefore, the effect of regularization is determined by the conditions of identification such as the number of measured parameters *m*, number of operating points *N*, measuring errors G and number of estimated parameters *r*.

#### *2.3. Numerical Simulation*

To assess the impact of the regularization coefficient α, the identification procedure was tested with the zero-dimensional thermodynamic model of a turboshaft engine [17] (Figure 2), based on experimental component maps. Performed numerical experiments involved calculation of N = 12 operating points for two scenarios:


**Figure 2.** Scheme and measured parameters of the engine.

Five gas path parameters Yi were used for model identification:


*Aerospace* **2020**, *7*, 6

(4) Rotational speeds of rotors N2 and N1.

Initial deviations of two compressor performance parameters were assumed as follows:


Four estimated parameters included:


Figure <sup>3</sup> presents estimations <sup>→</sup><sup>ˆ</sup> <sup>Θ</sup> and simulation results <sup>→</sup><sup>ˆ</sup> Y = F( → U, →ˆ Θ) for different values of α, which varied in a range from 0 to 1. ΔYav is the average deviation between initial and estimated models:

$$
\Delta \mathbf{Y\_{av}} = \sqrt{\frac{1}{\text{mN}} \sum\_{j=1}^{N} \sum\_{i=1}^{m} \left[ \mathbf{y\_i(\overrightarrow{\mathbf{U}\_{j\cdot}} \widehat{\Theta}) - \mathbf{y\_i(\overrightarrow{\mathbf{U}\_{j\cdot}} \overrightarrow{\Theta})} = 0 \right]^2}} \tag{11}
$$

and ΔY∗ av is the average deviation between measurements and estimated model:

$$
\Delta Y\_{\text{av}}^{\*} = \sqrt{\frac{1}{\text{mN}} \sum\_{j=1}^{N} \sum\_{i=1}^{m} \left[ \mathbf{y}\_{i} (\overrightarrow{\mathbf{U}}\_{\text{/}} \overset{\triangle}{\Theta}) - \mathbf{y}\_{\text{\textsuperscript{\text{\tiny}}}}^{\*} \right]^{2}} \tag{12}
$$

**Figure 3.** Influence of the regularization coefficient on estimations when noise is absent.

Figure 3 shows that despite the initial growth of Θ<sup>3</sup> and Θ<sup>4</sup> absolute values, the average value Θav decreases monotonically, whereas the residual between measurements and the corrected model increases. This corresponds to theoretical predictions of the influence of regularization on the identification process.

The simulation results help to choose the value of the regularization coefficient α. If the deviation of parameters <sup>→</sup> Y is to be lower than 5 % and the deviation of parameters Θ lower than 1%, then the value of α cannot exceed 0.03.

To obtain precise values of variance and to analyse the distributions of estimations <sup>→</sup><sup>ˆ</sup> Θ, a sequence of data sets with random measuring errors with variance σ<sup>2</sup> Yi was generated and identified 1000 times. This provided the average precision of about 1% of initial residuals for parameters Y and about 0.5 %

of simulated deviation 0.03 for the estimated parameters <sup>→</sup><sup>ˆ</sup> Θ.

Figure 4 presents histograms for calculations with two estimated parameters at α = 0 and α = 0.04. Both non-regularized and regularized estimations have distributions close to normal ones. Total scatter of estimations is also similar, but centres of regularized estimations moved from their true values by 0.008 and 0.005, respectively.

**Figure 4.** Distributions of estimations when two parameters are estimated.

When four parameters were estimated (Figure 5), estimations Θˆ <sup>1</sup> and Θˆ <sup>3</sup> preserved their distributions while estimations Θˆ <sup>2</sup> and Θˆ <sup>4</sup> had significant differences. Such abnormal behaviour of parameters Θ<sup>2</sup> and Θ<sup>4</sup> is explained by their correlation.

**Figure 5.** Distributions of estimations when four parameters are estimated.

At low values of the regularization coefficient α, centres of distributions are close to their true values, but scatter is high, and distributions are strongly asymmetrical. When α increases, then mean values move, but scatter decreases, and its distribution approaches the normal one.

As shown above, the value of the regularization coefficient α has to be set appropriately to provide meaningful results of regularized least squares. This procedure is useful in ill-conditioned engine model fitting but needs preliminary adjustment and sensitivity analysis.

#### *2.4. Regularized Identification Procedure Development Using a Priori Information*

This paper demonstrates an idea to implement engine heuristics involving component parameters and performances to improve the precision and stability of model identification. The main difficulty is in the diversity of this information, which is presented in one of the following forms:


The next difficulty is in the formalization of parameters that characterize the model quality. These parameters are set on the basis of subjective preferences of decision makers (DM). The same difficulties appear in the ranking of partial criteria and limitations according to their significance for estimation of the model quality. The analysis showed that the main problem is that theoretical-probabilistic methods are difficult to apply in the presence of uncertainties, which are related to subjective preferences whose nature is not statistical. Actually, the choice of the model's structure is the DM procedure, which in multi-criteria case is inevitably highly subjective. If the complexity of the task increases, the role of quality factors will grow. Therefore, it is possible to take into account all criteria using the proper mathematical tool.

We propose to use the fuzzy set theory [33] as this tool. This theory makes the uniform base to describe the information given in all the forms listed above, thus providing the correct mathematical definition of the identification task.

An example of applying fuzzy sets in the GPA and engine model matching is given by M. Zwingenberg et al. [34]. They used fuzzy logic for the evaluation of sensor failures. In contrast, we introduce a priori information about engine performance and experimental data directly into the stabilizing functional through the fuzzy logic approach.

Generalization of the functional (9) gives:

$$
\vec{\Phi(\Theta)} = \vec{\Phi\_{\mathfrak{e}}(\vec{\Theta})} + \vec{\Phi\_{\mathfrak{z}}(\vec{\Theta})} \tag{13}
$$

where Φe( → Θ) is the least squares functional (3), which minimizes measurement error, so it is called the empirical risk functional and Φa( → Θ) is the stabilizing functional, which is considered as the functional of a priori risk.

The determined a priori information is set as a limitation, which may take the form of an equality or inequality. The more general form of setting a priori information is its representation as a fuzzy set. The fuzzy set of parameter x is represented as a definition domain and a membership function μ(x) in this domain. For the considered task, the parameters x may be the model parameters <sup>→</sup> Θ or the output (calculated) variables <sup>→</sup> X.

Next, we will use limited normal fuzzy sets with limited definition domain (*xmin*, *xmax*) and 0 ≤ μ(x) ≤ 1. We will express each a priori information as a particular functional of a priori risk

Φa q( → θ) and will determine the general functional of a priori risk as a linear composition of particular functionals:

$$\Phi\_{\mathfrak{a}}(\overrightarrow{\theta}) = \sum\_{\mathfrak{q}=1}^{\mathbb{Q}} \alpha\_{\mathfrak{q}} \Phi\_{\mathfrak{a}\ \mathfrak{q}}(\overrightarrow{\theta}) \tag{14}$$

where α<sup>q</sup> are weight coefficients.

Let us consider some types of a priori information and its expression in view of fuzzy sets:

**Case 1.** Limitations of some parameters are known. For example, it is known that 0 < η < 1 and 0 < σ < 1, etc. Using experience and calculation results, these limits can be significantly reduced; for example, 0.5 < η < 0.9 and 0.9 < σ < 0.99 (Figure 6).

**Figure 6.** Examples of membership functions.

**Case 2.** The a priori mathematical model is known. This can be the model with design maps of components or the model of the average engine, which is matched with previous testing results.

This information may be expressed as a relationship between μ<sup>a</sup> or Φ<sup>a</sup> and the difference between parameters that correspond to the matched and a priori models. These parameters may be the model parameters <sup>→</sup> Θ as measured (for example fuel flow) or non-measured calculated parameters (for example thrust). The membership function, in this case, can be of symmetrical triangular shape and the functional of a priori risk Φa q( → θ) that characterizes the similarity between values of the parameter xq of matched and a priori models can be formed as

$$\Phi\_{\mathbf{a}\ \mathbf{q}}(\overrightarrow{\Theta}) = \sum\_{\mathbf{j}=1}^{N} \mu\_{\mathbf{q}} (\Delta \mathbf{x}\_{\mathbf{q}\ \mathbf{j}}) \cdot (\Delta \mathbf{x}\_{\mathbf{q}\ \mathbf{j}})^2 \,\Delta \mathbf{x}\_{\mathbf{q}\ \mathbf{j}} = \mathbf{x}\_{\mathbf{q}}(\overrightarrow{\mathbf{U}\_{\mathbf{j}}}, \overset{\frown}{\Theta}) - \mathbf{x}\_{\mathbf{q}\ 0}(\overrightarrow{\mathbf{U}\_{\mathbf{j}}}, \overset{\frown}{\Theta}) \tag{15}$$

where xq is a calculated parameter of the engine; xq( → Uj, →ˆ Θ) is the value calculated by the model to be matched; xq 0( → Uj, → <sup>Θ</sup>0)—the value calculated by the a priori model; <sup>→</sup> Θ<sup>0</sup> are the parameters of the a priori model.

**Case 3.** Information about the confidence in different sets of experimental data obtained in different conditions with different precision.

**Case 4.** Confidence in the available maps of the engine components. For example, we know that the compressor map used in the model corresponds to the old version of the engine and is far from the actual map. Thus, we can express this knowledge in view of the confidence functions that are in this case the membership functions.

The empirical risk functional Φe( → Θ) in Equation (3) may be considered as a square of Euclidian distance between two sets, one of which contains experimental data, and the other is composed of simulation results:

$$\Phi\_{\mathbf{e}}(\overrightarrow{\Theta}) = \frac{1}{\text{mN}} \sum\_{i=1}^{\text{m}} \text{g}\_{i} \sum\_{\mathbf{j}=1}^{\text{N}} \left[ \mathbf{y}\_{i}(\overrightarrow{\mathbf{U}}\_{\mathbf{j}}, \overrightarrow{\Theta}) - \mathbf{y}\_{i\mathbf{j}}^{\*} \right]^{2} \tag{16}$$

By analogy, the stabilizing functional can be formed as a second power of a distance from estimated parameter <sup>θ</sup><sup>k</sup> (or function x(<sup>→</sup> Θ) that is determined using the estimated parameters).

*Aerospace* **2020**, *7*, 6

The main problem in this analogy is that the fuzzy set that contains a priori information is infinite, so the sum in the above equation must be replaced with integral. For example, the second power of the distance from the engine map parameter θ = x to the fuzzy set with a given membership function μ(θ) equals:

$$\Phi\_{\mathbf{a}}(\stackrel{\textstyle \rightarrow}{\Theta}) = \frac{\int^{\infty} \mu(\theta)(\mathbf{x} - \theta)^{2} d\theta}{\int^{\infty} \mu(\theta) d\theta} \tag{17}$$

In the iteration process of the task solution, each integral must be calculated numerically. This requires a lot of time and high computational capacity. A numerical solution can be found for certain types of the membership function. We considered the trapezoidal function, the particular cases of which are the rectangle, triangle and isosceles trapezium (Figure 7).

**Figure 7.** Considered membership functions.

The functional derived for the trapezium:

$$\begin{split} \boldsymbol{\Phi}\_{\mathbf{a}} \cdot \mathbf{q} \stackrel{\leftrightarrow}{(\boldsymbol{\Phi}}) &= \mathbf{X}\_{\mathbf{q}}^{2} - \langle \mathbf{X}\_{\mathbf{q}} \Big[ \mathbf{x}\_{1} (\mathbf{a} + 2\mathbf{b} + \mathbf{c}) + \frac{1}{3} \Big( 2\mathbf{a}^{2} + 3\mathbf{b}^{2} + \mathbf{c}^{2} + 6\mathbf{ab} + 3\mathbf{ac} + 3\mathbf{b} \mathbf{c} \Big) \Big] + \\ &\frac{\mathbf{a}\mathbf{x}\_{1}^{2}}{2} + \frac{2}{3} \mathbf{a}^{2} \mathbf{x}\_{1} + \frac{\mathbf{a}^{3}}{4} + \mathsf{b} \Big( \mathbf{x}\_{1}^{2} + 2\mathbf{a} \mathbf{x}\_{1} + \mathsf{b} \mathbf{x}\_{1} + \mathsf{a}^{2} + \mathsf{ab} + \frac{\mathbf{b}^{2}}{3} \Big) + \\ &\frac{\mathbf{a}}{12} \Big[ 6\mathbf{x}\_{1}^{2} + (12\mathbf{a} + 12\mathbf{b} + 4\mathbf{c}) \mathbf{x}\_{1} + 6\mathbf{a}^{2} + 6\mathbf{b}^{2} + 12\mathbf{ab} + 4\mathbf{ac} + 4\mathbf{b} \mathbf{c} + \mathbf{c}^{2} \Big] / \left( \frac{\mathbf{a}}{2} + \mathsf{b} + \frac{\mathbf{c}}{2} \right). \end{split} \tag{18}$$

The proposed modifications of the objective functional cause the non-linearity of the estimation problem, so traditional solution methods are not applicable. This problem is overcome by adapting genetic optimization algorithm to the specifics of the engine model matching [35]. Genetic algorithms are increasingly used in gas turbines to solve complex optimization problems [4,10,36,37].

#### **3. Results**

The designed methods were implemented using experimental data obtained during the test-cell testing of the helicopter turboshaft engine. Table 1 contains the values of measured parameters. The last row presents mean squared measurement errors σi, which were used to form the diagonal weight matrix G (Equation (8)).


**Table 1.** Parameters acquired during bench testing.

The model parameters to be corrected:

πC\* (c)—Scaling factor of CPR (Compressor Pressure Ratio);

πC\* (a)—Factor of rotation in CPR-*W*Ccor plane of the compressor pressure map;

πC\* (b)—Factor of rotation in CPR-*N*cor plane;

*W*HPT(c)—HPT flow coefficient;

*W*PT(c)—Power turbine flow coefficient; ηC\* (a)—Factor of rotation in the ηC\* –*W*Ccor plane of the compressor efficiency map; ηC\* (c)—Scaling factor of the compressor efficiency map; ηCC—Scaling factor of combustion efficiency;

#### *3.1. Least Squares Identification*

The initial values of model parameters <sup>→</sup> Θ<sup>0</sup> that correspond to the average engine (a priori information) are shown in Table 2.



In the first step of the iteration procedure, the following corrections to component maps were determined:

$$
\delta\pi^\*\_{\mathsf{C}}(\mathsf{a}) = 342 \,\mathsf{\upmu}; \,\delta\mathsf{W}\_{\mathsf{H}\mathsf{PT}} = 0.93 \,\mathsf{\upmu}; \,\delta\mathsf{W}\_{\mathsf{PT}} = 0.85 \,\mathsf{\upmu}; \,\delta\mathfrak{n}^\*\_{\mathsf{C}}(\mathsf{a}) = 26 \,\mathsf{\upmu}; \delta\mathfrak{n}^\*\_{\mathsf{C}}(\mathsf{c}) = 1.4 \,\mathsf{\upmu} \tag{19}
$$

Such large corrections make it impossible to calculate part-load performance and influence coefficient matrices for the next iteration. Therefore, the calculation procedure was modified: Corrections to the component maps were limited. At corrections as high as 7%, the procedure slowly becomes convergent (by 7–9 iterations). The solution corresponds to small (less than 1%) corrections to the component maps.

The results of the standard least squares (the Newton–Raphson method) show that this method has low stability due to a weak conditionality of the engine model. Therefore, practical application requires improvement of the method's stability.

#### *3.2. Regularized Identification*

The same experimental data were processed by the regularized method based on Equation (9). Parameter α was set at 0.01. Table 3 shows corrections δ → Θ s for the first three iterations while Table 4 presents the parameters calculated by the corrected model.


**Table 3.** Corrections for three first iterations.

**Table 4.** Values of measured parameters determined by the corrected model.


Table 3 shows the high stability of the procedure: Estimations change slowly, without significant oscillations. The changes become smaller from iteration to iteration (this is a sign of stability). Comparison of Tables 1 and 4 shows that parameters determined by the corrected model are well correlated with the measured parameters. Thus, the developed procedure of component performance estimation and the engine model matching can be implemented in practice.

#### *3.3. Regularized Multi-Criteria Identification*

Finally, the above task of multi-mode diagnostics was considered with a priori information about the model of the average engine, the parameters of which were estimated by previous testing data (Table 2). It is also known that a scatter of measured parameters of different engines at *P*net = const, *N*<sup>1</sup> = const follows the normal distribution with the following mean squared error (Table 5).

**Table 5.** Mean squared error of performance parameters in the engine population.


This information was then transformed into the functional Φ<sup>a</sup> (Equation (13)), and taking into account that in this case, the functionals Φ<sup>e</sup> and Φ<sup>a</sup> are homogeneous as they contain residuals by parameters of the same names. This facilitated the setting of weight coefficients in Equation (14). They had to relate to a scatter of measurements and a scatter of the engine parameters by series. This provided the composition of the functionals:

$$\Phi(\stackrel{\textstyle \Theta}{\Theta}) = \Phi\_{\text{e}}(\stackrel{\textstyle \Theta}{\Theta}) + \frac{\sigma\_{\text{meas } \text{nv}}^{2}}{\sigma\_{0}^{2}} \Phi\_{\text{a}}(\stackrel{\textstyle \Theta}{\Theta}) \tag{20}$$

In the considered example <sup>σ</sup><sup>2</sup> meas av σ2 0 ≈ 8.

Table 6 shows the estimated parameters of the model in subsequent iterations. The final results of the corrected model are presented in Table 7. The comparison of Tables 1, 4 and 7 confirms that the proposed procedure is stable. For the corrected model, deviations of output parameters from experimental data are within the measurement error. The introduction of a priori information results in a smaller deviation of results compared to the conventional regularization method.


**Table 6.** Corrections for three first iterations.

**Table 7.** Values of measured parameters determined by the corrected model.


Figure 8 presents the experimental data and modelling results for part-load performance. As the initial model is far enough from the field data, the diagnostic algorithm needs a preliminary adjustment to make the model closer to the average engine and to improve the precision of GPA. The visible lines correspond to the initial model before the correction, the average engine model and corrected model of the analyzed engine. The deviations of performance parameters from the reference indicate the engine's health status.

**Figure 8.** Experimental data vs models. (**a**) fuel flow; (**b**) high pressure rotor rotation speed; (**c**) turbine inlet temperature; (**d**) compressor pressure ratio.

Figure 9 illustrates the correction of compressor maps and confirms that the procedure was effective even though the initial model was inaccurate.

**Figure 9.** Compressor maps with the operating lines.

#### **4. Discussion**

Matching turbine engine models to experimental data is an inverse problem of mathematical modelling which is often characterized by parametric uncertainty. This results from the fact that the number of measured parameters is significantly lower than the number of components' performance parameters needed to describe the real engine. It was shown that in such conditions, even small measurement errors resulted in a high variation of results, hence the obtained efficiency, loss factors, etc., exceeded physically feasible ranges.

It was confirmed that extending the least squares functional by regularizing term improved stability of estimation but caused significant bias, which could lead to incorrect diagnoses. Stable operation of the algorithm and tolerant bias of the estimates in wide range of engine operation conditions are difficult to achieve as only selected values of the regularization coefficient α can ensure that. Finding its optimal solution is complicated because the regularizing term is not related to the physical properties of the object.

The new regularization procedure of multi-criteria identification, introduced in Section 2.4, uses a priori information about the engine, its measurement system, mathematical model and expected performance (Figure 1). This information is heterogeneous: Partially, it is presented in the form of statistic parameters and partially in the form of heuristic expert representations. Some other methods, e.g., Extended Kalman Filter can also use a priori information [22]. Undoubtedly, the quality of implemented expert knowledge should be carefully checked in advance.

The framework that supported the formalization of the regularization problem used fuzzy sets theory [33,37]. Similarly to the solution of Tanaka et al. [38], membership functions were introduced directly into the minimized functional, which contained a priori information about an acceptable range of estimated parameters. However, fuzzy logic is not used here for the non-linear mapping of parameters like in AI-based diagnostics [39]. In this work, inputs and outputs are measured variables with Gaussian noise.

The integral form of the stabilization functional was proposed which operates with continuous membership functions. To accelerate numerical calculations, the trapezoidal membership function was recommended. The analytical equation was derived for the stabilization functional with this function, which was used in place of numerical integration, thus decreasing the computing capacity required for the estimation procedure.

In the example presented in Section 3.3, available component maps were significantly different from the actual maps of the tested engine. Nevertheless, the estimation algorithm worked stable and provided a perfect matching of the model to the experimental data. Thanks to regularizing terms implemented in the estimation procedure, intermediate estimations did not exceed their limits when the computational algorithm was unstable.

Further studies and more field data are needed to quantify the accuracy and the stability of the procedure in relation to the number of estimated parameters. However, when a large number of engine parameters is measured, measurement uncertainty is low, and the number of coefficients of the model to be found is significantly less than the number of measured parameters, the estimates will not exceed the acceptable range, stability is not a problem, and the regularization does not have a positive effect and is not required. Hence, the proposed method is efficient when the number of measured parameters is small, and errors of measurement are significant. This is often the case in real applications.

#### **5. Conclusions**

A new method for the stable estimation of engine performance parameters is proposed. It uses a priori information about the engine and data acquisition system, i.e., the expected performance and acceptable variations of parameters. The method improves the accuracy of gas-turbine performance models in off-design conditions. A priori information is introduced in the model identification procedure in the form of membership functions, which are transformed into additional stabilization terms of the functional to be minimized.

The proposed approach was used to process test-cell data from a helicopter engine. The comparison with the conventional least squares and Tikhonov regularization proved the following:


**Author Contributions:** Conceptualization, O.K. and R.Z.; methodology, O.K.; software, R.Z.; validation, R.P. and S.Y.; formal analysis, S.Y. and R.Z.; investigation, O.K.; writing—original draft preparation, S.Y.; writing—review and editing, R.P.; visualization, S.Y., R.Z. and R.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This publication was prepared within the framework of the AERO-UA project, which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 724034, and with the support of the Ministry of Education and Science of Ukraine (research project No. D203-3/2019-Π).

**Acknowledgments:** The authors would like to acknowledge the constructive comments, and suggestions provided by the anonymous reviewers that greatly improved the quality of the article. We are also grateful to Michail Ugryumov and Anatolyi Mazurkov for their help in the application of the genetic optimization procedure and Kacper Grz ˛edzi ´nski for his comments on an earlier version of the manuscript.

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

#### **Nomenclature**


#### **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* **Superhydrophobic Coatings as Anti-Icing Systems for Small Aircraft**

**Filomena Piscitelli 1,\*, Antonio Chiariello 1, Dariusz Dabkowski 2, Gianluca Corraro 1, Francesco Marra <sup>3</sup> and Luigi Di Palma <sup>1</sup>**


Received: 26 November 2019; Accepted: 25 December 2019; Published: 2 January 2020

**Abstract:** Traditional anti-icing/de-icing systems, i.e., thermal and pneumatic, in most cases require a power consumption not always allowable in small aircraft. Therefore, the use of passive systems, able to delay the ice formation, or reduce the ice adhesion strength once formed, with no additional energy consumption, can be considered as the most promising solution to solve the problem of the ice formation, most of all, for small aircraft. In some cases, the combination of a traditional icing protection system (electrical, pneumatic, and thermal) and the passive coatings can be considered as a strategic instrument to reduce the energy consumption. The effort of the present work was to develop a superhydrophobic coating, able to reduce the surface free energy (SFE) and the work of adhesion (WA) of substrates, by a simplified and non-expensive method. The developed coating, applied as a common paint with an aerograph, is able to reduce the SFE of substrates by 99% and the WA by 94%. The effects of both chemistry and surface morphology on the wettability of surfaces were also studied. In the reference samples, the higher the roughness, the lower the SFE and the WA. In coated samples with roughness ranging from 0.4 and 3 μm no relevant variations in water contact angle, nor in SFE and WA were observed.

**Keywords:** superhydrophobic; coating; anti-icing; spray-coat; aeronautical

#### **1. Introduction**

Typically, the presence of tiny pieces of ice or supercooled liquid water in the clouds, which remain liquid below zero and suddenly turn to ice after the impact with the aircraft surfaces, are the main sources of ice deposition during a flight. The presence of ice on surfaces alters the airflow over the wing and tail, and then reduces the lift force that keeps the plane in the air. This potentially causes aerodynamic stall, a condition that can lead to a temporary loss of control of the aircraft. In order to prevent or reduce the ice formation or alternatively to remove the ice once it was formed, anti-icing and de-icing systems are usually adopted. Currently the thermal and pneumatic types represent the most largely employed systems. In details, the thermal system melts the ice accretion or prevents the ice from forming by the application of heat on the protected surface of the wing. The heat is generated by the hot air "bled" off the jet engine into piccolo tubes routed through wings, tail surfaces, and engine inlets. The spent bleed air is then exhausted through holes in the lower surface of the wing. Similarly, the electro-thermal systems use resistive circuits buried in the airframe structure to generate heat when a current is applied. The heat can be generated continuously to protect the aircraft from icing

(anti-ice mode), or intermittently to shed ice as it accretes on key surfaces (de-ice), the latter being generally preferred due to the lower power consumption. Instead, the pneumatic de-ice boot consists of a rubber sheet bonded to the leading edge of the airfoil: when the ice builds up on the leading edge, an engine-driven pneumatic pump inflates the rubber boots, and then the ice is cracked and should fall off the leading edge of the wing. Additionally, some innovative systems are being developed, e.g., the electro-mechanical de-icing systems which use a mechanical force, generated for instance by actuators, to knock the ice off the flight surface, or the weeping wing system, which releases a glycol-based chemical onto the wing surface through small orifices which causes the detachment of ice. However, these methods can be inefficient, environmentally unfavorable, expensive and time consuming [1]. Thus, it would be advantageous if surfaces could passively prevent the ice formation and facilitate ice removal. In this contest, the superhydrophobic coatings, applied as usual paint, owing to their extraordinary water repellency, and not requiring additional energy consumption, can be viewed as excellent candidates for icephobicity in this area [2,3]. This is because, ultrahigh θ*rec*, typical of nanostructured superhydrophobic surfaces, means low ice adhesion strength [4,5]. Subsequently, if supercooled liquid water freezes when it impacts the solid surface of an aircraft, the resulting ice can be shed, taking advantage of external forces, such as aerodynamic forces, to overcome ice-surface adhesion forces [6]. Evidently, the efficiency of these coatings as an anti-icing passive method largely depends on the different scenarios in which the ice may form [1]. In some cases, therefore, a combination of a traditional icing protection system and the superhydrophobic coatings could be seen as a strategic instrument able to assure a high efficiency in a wide range of environmental conditions.

In the last decade, a huge body of literature concerning the development of superhydrophobic and icephobic coatings able to delay the ice formation [7], or reduce the ice adhesion strength once formed, has been produced [8–15]. A shared strategy is to create surfaces with low surface free energy and, contemporary, micro-nanoscopic rough structures. In this regard, it was found that the icing probability was reduced to zero with particle diameters ranging between about 30 and 70 nm [16]. Some authors achieved this goal by using etching and/or high temperatures [9,11], or complex fabrication techniques [12–17], which are often limited by the substrate type and geometry that can be successfully coated [18].

The effort of the present work was to develop a superhydrophobic coating for metallic substrates with a simplified and non-expensive method, which could be employed as a usual paint able to prevent/reduce the formation of ice, especially on small aircraft [19]. The newly formulated superhydrophobic coating consists of nanostructured layers able to generate hierarchical micro/nano-structured roughness, and reduce the surface free energy, which as previously declared, are the two main factors useful to making a superhydrophobic surface. The effect of the substrates' roughness on the surface properties of coated and uncoated samples was also studied, and then correlated to the wettability of samples, in order to identify and separate the morphological and the chemical contributions.

This work is being developed in the framework of the Clean Sky 2 ongoing SAT-AM (More Affordable Small Aircraft Manufacturing) project, whose main goal is the investigation of new technologies for a future small aircraft able to fly with low fuel consumption, low noise levels and needing low quantities of raw material in its life cycle. The reference vehicle for this project is the M28 designed and manufactured by Consortium Partner Polskie Zakłady Lotnicze (PZL), Mielec (PL). It is a commuter category 19 passenger, twin-engine high-wing cantilever monoplane of all-metal structure, with twin vertical tails and a robust tricycle non-retractable landing gear, featuring a steerable nose wheel to provide for operation from short, unprepared runways where hot or high-altitude conditions may exist. The M28 is best suited for passenger and/or cargo transportation and is certified under EU CS-23 and USA FAR 23 requirements.

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

The developed superhydrophobic coating was prepared starting from the formulation described in a previous work [19] and slightly modified in order to improve its durability. In this regard, other details about the performed changes cannot be disseminated in this manuscript. Once prepared, several layers of the developed coating were applied with an aerograph on substrates representative of the M28 air intake. The dehumidified air was pressured at three bar; whereas the other application parameters, such as the layers' number, the curing temperature and the distance between the aerograph and samples were optimized in order to guarantee uniformity in the coating's thickness and weight.

The substrate samples, 5 <sup>×</sup> 5 cm2 in dimensions, made of stainless steel 12H18N10T (same composition of the M28 air intake), and having roughness ranging from 0.7 to 4.6 μm, were provided by METROL (a partner of the Clean Sky project). The effect of the substrate roughness on the final properties of the coating in terms of wettability and adhesion of the coating on the substrate were studied.

Roughness of substrates before and after the application of the coating was measured by employing a SAMA SA6260 surface roughness meter. Roughness measures, performed according to the ISO 4288 [20], were reported as Ra, which represents the arithmetic average of the absolute values of the profile height deviations from the mean line. ID samples was correlated to the correspondent roughness values in Table 1.


**Table 1.** ID samples with roughness ranging between 0.7 and 4.6 μm.

Contact angles (CA) were calculated according to Young's equation [21] (see Figure 1):

$$
\gamma\_{\rm LV} \cos \theta = \gamma\_{\rm SV} - \gamma\_{\rm SL} \tag{1}
$$

where θ is the contact angle, γ*LV* and γ*SV* are the liquid and solid surface free energy, respectively, whereas γ*SL* is the solid/liquid interfacial free energy. The schematic illustration of the equilibrium among involved forces according to the Young equation is shown in Figure 1.

**Figure 1.** Schematic illustration of an ink drop on a solid substrate with the surface free energy and the contact angle as described by Young's equation [21].

The contact angle measurements were performed at room temperature in compliance with the ASTM D7490–13 [22] standard, using water (H2O) and diiodomethane (CH2I2). The Surface Free Energy (SFE) was calculated according to the Owens Wendt (OW) method [23,24], for which the surface energy of a solid is the sum of two components, a dispersive and a polar one. The dispersive component theoretically accounts for the Van der Waals and other non-site specific interactions between the surface and the applied liquid, whereas the polar component accounts for the dipole-dipole, dipole-induced dipole, hydrogen bonding, and other site-specific interactions [25]. The polar and dispersive components of the reference liquids are listed in Table 2. Tests were carried out by depositing ten drops of 3 μL of each liquid on the sample surface.


**Table 2.** Surface tension components of the reference liquids [26].

According to the OW method [23], the interfacial solid/liquid energy can be evaluated as:

$$\gamma\_{\rm s1} = \gamma\_{\rm s} + \gamma\_1 - 2\left(\gamma\_1^{\rm d}\gamma\_{\rm s}^{\rm d}\right)^{1/2} - 2\left(\gamma\_1^{\rm F}\gamma\_{\rm s}^{\rm F}\right)^{1/2} \tag{2}$$

which, combined with the Young equations, gives the following:

$$\gamma\_1(1+\cos\theta) = 2\left(\gamma\_{1\_1}^{\mathrm{d}}\gamma\_{\mathrm{s}}^{\mathrm{d}}\right)^{1/2} + 2\left(\gamma\_{1\_1}^{\mathrm{p}}\gamma\_{\mathrm{s}}^{\mathrm{p}}\right)^{1/2} \tag{3}$$

Therefore, the polar (γ*<sup>p</sup> s*) and dispersive (γ*<sup>d</sup> <sup>s</sup>*) components of the solid SFE are given as solution of the following non-linear system:

$$\begin{cases} \frac{\gamma\_{\mathbb{I}\_1}(1+\cos\vartheta\_1)}{2} = \left(\gamma\_{\mathbb{I}\_1}^{\mathrm{d}}\gamma\_{\mathrm{s}}^{\mathrm{d}}\right)^{1/2} + \left(\gamma\_{\mathbb{I}\_1}^{\mathrm{P}}\gamma\_{\mathrm{s}}^{\mathrm{P}}\right)^{1/2} \\\ \frac{\gamma\_{\mathbb{I}\_2}(1+\cos\vartheta\_2)}{2} = \left(\gamma\_{\mathbb{I}\_2}^{\mathrm{d}}\gamma\_{\mathrm{s}}^{\mathrm{d}}\right)^{1/2} + \left(\gamma\_{\mathbb{I}\_2}^{\mathrm{P}}\gamma\_{\mathrm{s}}^{\mathrm{P}}\right)^{1/2} \end{cases} \tag{4}$$

where γ*<sup>d</sup> l*1 , γ*<sup>p</sup> l*1 , γ*<sup>d</sup> l*2 , γ*<sup>p</sup> <sup>l</sup>*<sup>2</sup> and <sup>θ</sup>1, <sup>θ</sup><sup>2</sup> are the SFE components of the two reference liquids and the corresponding contact angles between the solid and the liquids, respectively.

Although the SFE's components can be known through the literature (Table 2), the contact angles have to be experimentally evaluated through several measures. That determines a stochastic nature of the contact angles and therefore of the solutions of the equation system (4).

So, to assess the SFE problem, a third reference liquid, i.e., the Formamide (HCONH2), was introduced in order to reformulate the initial problem in terms of the following nonlinear programming problem:

min γ*d S*,γ*p S* γ*d l*3 γ*d s* 1/2 + γ*p l*3 γ*p s* 1/2 <sup>−</sup> <sup>γ</sup>*l*<sup>3</sup> (<sup>1</sup> <sup>+</sup> cos <sup>θ</sup>3) 2 (5) *s*.*t*. ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ γ*d l*1 γ*d s* 1/2 + γ*p l*1 γ*p s* 1/2 − γ*l* <sup>1</sup> (1+cos <sup>θ</sup>1) <sup>2</sup> = 0 γ*d l*2 γ*d s* 1/2 + γ*p l*2 γ*p s* 1/2 − γ*l* <sup>2</sup> (1+cos <sup>θ</sup>2) <sup>2</sup> = 0 E[θ1] − σθ<sup>1</sup> ≤ θ<sup>1</sup> ≤ E[θ1] + σθ<sup>1</sup> E[θ2] − σθ<sup>2</sup> ≤ θ<sup>2</sup> ≤ E[θ2] + σθ<sup>2</sup> E[θ3] − σθ<sup>3</sup> ≤ θ<sup>3</sup> ≤ E[θ3] + σθ<sup>3</sup> (5a)

where a Gaussian distribution of the contact angles was assumed in order to statistically characterize them through the mean value (E[θ1], E[θ2], E[θ3]) and the standard deviation (σθ<sup>1</sup> , σθ<sup>2</sup> , σθ<sup>3</sup> ).

The stated optimization problem allowed determining the θ1, θ2, θ<sup>3</sup> contact angles that identify the SFE components, which minimize the third Equation (5), taking into account the measure variability and the constitutive Equations (4) of the considered liquids.

Applying the described method to the reference liquids reported in Table 2, the constraints (5a) were linearized simplifying the complexity of the problem that has been solved in Matlab/Simulink environment. The developed Matlab tool provides (as output) the values of the three best contact angles θ1,θ2,θ3, along with the values of the polar and dispersive components of the solid's SFE.

As an additional crosscheck, the determined optimal contact angles, namely, θ1,θ2, instead of the mean values of the experimentally measured angles measured with Water and Diiodomethane, were introduced in the widespread linearized equations [19–26]:

$$\begin{cases} \frac{\gamma\_{\overline{\prime}\_{1}} \left( 1 + \cos \overline{\vartheta}\_{1} \right)}{2 \sqrt{\gamma\_{\overline{\prime}\_{1}}^{d}}} = \frac{\sqrt{\gamma\_{\overline{\prime}\_{1}}^{p}}}{\sqrt{\gamma\_{\overline{\prime}\_{1}}^{d}}} \cdot \sqrt{\gamma\_{s}^{p}} + \sqrt{\gamma\_{s}^{d}}\\ \frac{\gamma\_{\overline{\prime}\_{2}} \left( 1 + \cos \overline{\vartheta}\_{2} \right)}{2 \sqrt{\gamma\_{\overline{\prime}\_{2}}^{d}}} = \frac{\sqrt{\gamma\_{\overline{\prime}\_{2}}^{p}}}{\sqrt{\gamma\_{\overline{\prime}\_{2}}^{d}}} \cdot \sqrt{\gamma\_{s}^{p}} + \sqrt{\gamma\_{s}^{d}} \end{cases} \tag{6}$$

and the required polar and dispersive components of the solid SFE, were obtained as the graphic solution of the equations system (6) through the intercept point between the two straight line (*r*<sup>1</sup> and *r*2) identified by the following angular coefficients and points:

$$r\_1: \quad m\_1 = \frac{\sqrt{\gamma\_{I1}^p}}{\sqrt{\gamma\_{I1}^d}}; \quad P\_1 = \left(0 \quad ; \frac{\gamma\_{I1} \left(1 + \cos \overline{\theta}\_1\right)}{2\sqrt{\gamma\_{I1}^d}}\right) \tag{7a}$$

$$r\_{\mathsf{T}}:\quad m\_{\mathsf{T}} = \frac{\sqrt{\mathsf{\gamma}\_{\mathsf{I}\mathsf{2}}^{p}}}{\sqrt{\mathsf{\gamma}\_{\mathsf{I}\mathsf{2}}^{d}}};\quad P\_{\mathsf{T}} = \left(0:\frac{\mathsf{\gamma}\_{\mathsf{I}\mathsf{2}}\left(1+\cos\overline{\theta}\_{\mathsf{2}}\right)}{2\sqrt{\mathsf{\gamma}\_{\mathsf{I}\mathsf{2}}^{d}}}\right) \tag{7b}$$

Once calculated the polar and dispersive components, the required solid SFE was assessed, as:

$$
\gamma\_s = \gamma\_s^p + \gamma\_s^d \tag{8}
$$

whereas the Work of Adhesion (*W*A) can be calculated as:

$$W\_{\rm A} = \gamma\_{\rm l} (1 + \cos \theta) \tag{9}$$

or equally as:

$$\mathcal{W}\_{\rm A} = 2\left(\gamma\_s^{\rm P}\gamma\_1^{\rm P}\right)^{1/2} + 2\left(\gamma\_s^{\rm d}\gamma\_1^{\rm d}\right)^{1/2} \tag{10}$$

Finally, the Surface Polarity (SP) was calculated as:

$$SP = \frac{\mathcal{V}^p}{\mathcal{V}^p + \mathcal{V}^d} \tag{11}$$

The Young Equation (1) is valid for smooth surfaces or substrate having very low roughness, but when the roughness is too high to be neglected, the equation should be corrected taking into account the actual surface morphology. Two different models can be useful to describe the wetting of textured surfaces, i.e., the Wenzel and Cassie–Baxter (Figure 2a,b respectively) models. In details, the Wenzel model [27,28] proposed a correction factor "*r*" for contact angle on rough surfaces, which is equal to the ratio of rough interfacial area over flat interfacial area under the droplet. Substituting the roughness ratio factor "*r*" in the Young's Equation (1), one can obtain the Wenzel's equation [28], i.e.,

$$r\cos\theta \*= r\cos\theta\tag{12}$$

where θ*\** and θ are the contact angles of a droplet on a rough surface and on the corresponding smooth surface, respectively. Wenzel's model assumes no air-trapping under droplet, which may not necessary be true.

**Figure 2.** (**a**) Wenzel and (**b**) Cassie–Baxter models useful to describe the wetting of textured surfaces.

On the contrary, the Cassie–Baxter model [29] allows to estimating the contact angle on rough surface with air-trapping, through the following equation:

$$\cos\theta \mathbf{\*} = -1 + \varphi\_{\mathbb{S}}(1 + \cos\theta) \tag{13}$$

where ϕ*<sup>S</sup>* is the area ratio of liquid-air interface to the whole interface; θ\* and θ are the contact angles with and without considering air-trapping.

Cutting and tape test were carried out, according to the ASTM D 3359-09 standard [30], by making a grid incision with a specific cutter on the coated samples and then by applying an adhesive tape to cover the cut area. The test adhesive was vigorously removed and the involved area was inspected. According to the ASTM D 3359-09 [29], the test passes if the area involved in detached flakes of coating at intersections is less than 15%.

Sample morphology, measurements of nanoparticle size and of coating thickness, were carried out using a field emission SEM Tescan Mira3 (Tescan Orsay Holding, Brno-Kohoutovice, Czech). Samples were observed in as-realized conditions, without metallization on the observed surface.

Surface 3D morphology and the corresponding roughness were assessed by using the laser profilometer Ametek Talyscan 150. The scanning area for 3D surface reconstruction was 40 <sup>×</sup> 40 mm<sup>2</sup> (800 profiles with 5 microns resolution), acquired with a scanning velocity of 10,500 microns/s.

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

#### *3.1. Substrates' Roughness of Substrates*

Roughness of substrates measured with the SAMA SA6260 surface roughness meter, before and after the application of the coating, were compared in Figure 3. It highlights that, in general, the application of the coating reduces the original roughness. This especially for samples 1–6, for which the values of the coated samples' roughness vary in the range 0.4–0.6 μm, in spite of the substrates' roughness ranging from 0.7 to 1.5 μm, whereas, for samples 7–9, the application of the coating reduces the original roughness from 4–4.5 to about 3 μm.

**Figure 3.** Roughness values of samples before (triangles) and after (circles) the application of the coating.

#### *3.2. Application of Coating*

The aerograph setting parameters were varied and optimized in order to guarantee uniformity in terms of thickness and weight of the applied coating. Here, as an example, two different setting parameters, labelled as *Procedure 1* and *Procedure 2*, are described; the corresponding schematics are shown in Figure 4. They differ mostly in the samples orientation with respect to the aerograph spray and in the drying temperature, whereas the layers' number, the coating formulation, and the substrates were the same. After the application, measures of contact angle with water at room temperature were carried out. It was found that Procedure 1 is able to produce samples with contact angle of 124◦, whereas Procedure 2 gives surfaces with a water contact angle of 155◦ (Figure 5). This difference is probably due to the different morphologies because of the diverse application procedures.

**Figure 4.** Schematic drawing of deposition process according to (**a**) Procedure 1 and (**b**) Procedure 2.

#### **Procedure 1**


#### **Procedure 2**


Hence, *Procedure 2* was adopted to apply three layers of the developed coating on sample 1–9 with dimensions of 5 <sup>×</sup> 5 cm2. The Samples were weighed before and after the application of the coating, in order to assess a rough estimation of the coating's specific weight, which was plotted as a function of the sample roughness in Figure 6. It was found that the coating's specific weight was about 11 <sup>±</sup> 1 g/m2.

**Figure 5.** Pictures related to the water contact angle of coated samples prepared according to (**a**) Procedure 1 and (**b**) Procedure 2.

**Figure 6.** The coating's specific weight as a function of the substrates' roughness.

SEM cross-section images of sample 2 (*R*a = 0.8 μm) and sample 7 (*R*a = 4 μm), employed to measure the coating's thickness, are shown in Figure 7a,b, respectively. It seems that the two samples display some differences in both coating thickness and morphology. Sample 7 in Figure 7b, in fact, has a coating thickness higher than sample 2 shown in Figure 7b (see values in Table 3), with a mean value of 190 μm vs. 127 μm. Moreover, sample 7 coating appears to be more uniform in thickness and morphology than the sample 2 coating, and with less amount of large and localized voids. Hence, it seems that the roughness is able to improve the coating's uniformity.

**Figure 7.** SEM images acquired in cross-section mode of samples 2 (**a**) and 7 (**b**). Some measured thicknesses are reported.

**Table 3.** Coating's thickness measured by the SEM acquired in cross section mode.


#### *3.3. Contact Angle Measurements*

The measures of the contact angles acquired with H2O and CH2I2 were reported as function of the substrate roughness in Figures 8 and 9, respectively. Some representative drop pictures caught during the measures were overlapped. It should be noted that, for the reference samples, the water contact angle is almost constant with a substrate roughness up to 1.8 μm beyond that, as expected [31], it increases as the roughness increases, reaching a value of about 100–115◦ for the highest value of measured roughness, i.e., 4.6 μm. The increasing of the contact angles with the roughness of metallic substrates can be ascribed to the changes in morphology. Other authors [31,32] also observed this behavior with maximum values of the achieved contact angles at about 120 ◦C [31]. Indeed, it was found that in order to achieve higher values of contact angles without modifying the surface chemistry, it is necessary to scale down the surface roughness into the sub-micron range, by creating specific and regular pillar-like structures having circular, square, triangle, crossed or fractal-like shapes [33]. On the contrary, the combination of substrate roughness and low SFE can be seen as a synergic approach able to achieve contact angle mean value of 158◦ (Figure 8). Additionally, it was also found that both the original and the actual sample roughness (Figures 3 and 8) do not affect the achieved contact angle, since, in spite of the roughness value, the contact angle value does not change.

The same uniformity in the achieved values of the contact angles was observed for measurements performed with CH2I2 in coated samples (Figure 9). Here, the mean value of the contact angles in coated samples is 148◦, regardless of the substrates (0.7–4.6 μm) and the actual roughness (0.3–3.2 μm) of the samples.

**Figure 8.** Water contact angles of reference and coated samples as a function of the substrates' roughness. Corresponding pictures are also shown.

**Figure 9.** CH2I2 contact angles of reference and coated samples as a function of the substrates' roughness. Corresponding pictures are also shown.

#### *3.4. Surface Free Energy and Work of Adhesion*

The solid SFE and the WA of both uncoated and coated samples were plotted in Figure 10a,b, respectively, as a function of the substrates' roughness.

**Figure 10.** SFE and WA of the uncoated (**a**) and coated (**b**) samples as a function of the roughness.

It was found that the solid SFE of the reference samples varies between 44 and 54 mN/m, for roughness values ranging from 0.7 to 2 μm, whereas for roughness of 4 μm the SFE becomes 38 mN/m (Figure 10a). Similarly, the WA varies between 108 and 122 mN/m, assuming a value of 69 mN/m for the sample at 4 μm in roughness. For samples having roughness higher than 4 μm, results cannot be achieved without taking into account the actual surface roughness, e.g., through Equations (12) and (13), or more complex relationship, i.e., those described in ref. [33]. Indeed, it was found that roughness higher than 4 μm the Equations' system (4) did not give a solution, so no intersection of the two curves can be observed in graphs of Figure 11b. On the contrary, for lower values of the contact angles, the Equation system (4) reached solutions, as shown in Figure 11a as an example.

**Figure 11.** Graphical solution of the SFE equations' system with (**a**) and without (**b**) solution.

The authors decided not to further investigate samples at high values of roughness, since values higher than 4 μm are out of interest for aeronautical applications.

On the other hand, the flat trend for SFE and WA was seen for coated samples (Figure 10b), since all values scatter around 0.38 mN/m for the SFE and around 6.1 mN/m for the WA. This result reflects the uniformity of the water contact angle values measured in spite of the original and actual roughness values: the latter changing between 0.3 and 3 m (Figure 3). These trends can be better highlighted if results are reported as a function of the solid SFE (Figure 12). For the reference samples (Figure 12a), it was found that, as expected, the higher the contact angle, the lower the SFE and the WA. Since this result was achieved by varying the roughness of substrates alone, it can be concluded that it can be ascribed to the surface morphology (MORPHOLOGICAL EFFECT). On the contrary, for coated samples, no changes in terms of water contact angle and WA, along with almost no variations in the solid SFE (Figure 12b) can be observed, in spite of the actual roughness of substrates. Therefore, it should be concluded that, for the analyzed samples, the chemical modification of the surfaces prevails on the morphological aspect (CHEMICAL EFFECT).

**Figure 12.** Water contact angle and WA as a function of (**a**) the surface free energy for the reference and (**b**) the coated samples.

The geometrical and chemical effects on the wettability of surfaces were also observed by Kam et al. [31]. They found that the maximum contact angle achieved with nano-microstructured metallic surfaces was 110◦, whereas the combination of the surface geometry and its chemical modification with silanes allowed to additionally modify the wettability of surfaces, since the contact angles reached values higher than 140◦.

Finally, in Figure 13 the polar and dispersive components of the surface free energy, and the WA assessed for the reference and the coated samples were compared. It highlights that the WA was reduced by 94%, the solid SFE by 99%, and the SP by 100%, since after the application of the coating, it changes from 34 to 0%.

**Figure 13.** Polar and dispersive components of the SFE and the WA for the references and the coated samples.

#### *3.5. 3D samples' Morphology*

The SEM images of coated samples 2 and 7 (Table 1), acquired at different magnifications (125–500,000×) are displayed in Figures 14 and 15, respectively. It highlights that at the lowest magnification (125×), both samples exhibit characteristic sponge-like structures [15,19,32,34–36], with the formation of air pockets, which enhance the superhydrophobicity. Samples differ in the density of the sponge-like structures. In fact, it appears that sample 7, having higher roughness (i.e., 2.8 μm) generates air pockets that are more densely packed than the less rough sample 2 (*R*a = 0.4 μm). In spite of this difference, both samples exhibit contact angles of 158◦ (see Figure 8). At the highest magnification (500,000×) (Figures 14 and 15), the hierarchical micro- and nano-structures appears, highlighting for both samples the dimension of the nanoparticles employed to formulate the coating, namely, about 30 nm in diameter. Therefore, in spite of the intrinsic roughness and the macroscopic density of the sponge-like structure, the generated air pockets made of multiscale roughness (microand nano-meter in dimension) guarantee the superhydrophobicity, i.e., contact angles of 158◦.

**Figure 14.** SEM images at different magnifications of a sample having an original roughness of sample 2 (**a**) 125×, (**b**) 750×, (**c**) 2500×, (**d**) 500,000×.

The surface 3D profilometry of coated samples 2 and 7, having actual roughness of 0.4 and 2.8 μm, respectively, were reported in Figure 16a,b; whereas Figure 16c,d show the corresponding references, having roughness of 0.8 and 4.0 μm, respectively. In the insert, the pictures representative of the water contact angles were also shown. It highlights that all samples exhibit a macroscopic needle-like morphology deriving from the specific textured roughness, and there are no differences in the needle-like morphology [36] for each series of samples (compare Figure 16a with Figure 16c, and Figure 16b with Figure 16d). The reference with low roughness (see Figure 16c) does not exhibit the Cassie state, the surface area in contact with the water droplet is high, and consequently the water contact angle is low (50◦). When the reference roughness increases to about 4 μm (see Figure 16d), the relative distance between the surface irregularities and the related height allow the contact angle to increase to 115◦. The increase of the water contact angle due to the increased roughness alone can be ascribed to the MORPHOLOGICAL EFFECT. By comparing coated samples with the corresponding references, it highlights that, although the needle-like morphology and the corresponding roughness do not change, the water contact angle increases from 50◦ and 115◦, respectively, to 158◦, as a consequence

of the chemical modification of the surface (CHEMICAL EFFECT). Finally, by comparing the two coated samples (see Figure 16a,b), it must be noticed that, in spite of the different actual roughness, i.e., 0.4 and 2.8 μm, for the samples 2 and 7, respectively, both of them exhibit the same water contact angle, namely, 158◦. Hence, in conclusion, it seems that the chemical effect prevails on the morphological effect.

**Figure 15.** SEM images at different magnifications of a sample having an original roughness of sample 7 (**a**) 125×, (**b**) 750×, (**c**) 2500×, (**d**) 500,000×.

**Figure 16.** Surface 3D profilometry obtained with the laser for coated sample 2 (**a**) and sample 7 (**b**); in (**c**) and (**d**) the corresponding references were also reported. In the insert, the pictures representative of the related water contact angles are shown.

#### *3.6. Cutting and Tape Test*

Figure 17 shows images related to the cutting and tape test performed on coated sample 7. In detail, Figure 17a shows the picture after the cut and Figure 17b shows magnification of the surface after the test. The shape of the water droplet in Figure 17c highlights that the superhydrophobicity was preserved after the test too. It must be noted that, according to the ASTM D 3359-09 [30], the test was passed, since quite no detached flakes of coating could be observed. This result was classified as 5B [30].

**Figure 17.** Cutting and tape test results on coated sample 7, (**a**) surface cut, (**b**) surface magnification, (**c**) surface wettability after test.

#### *3.7. Low Temperature Wettability*

A preliminary test about the wettability at low temperature of the developed coating was performed by freezing (at −27 ◦C) a few water droplets placed on both reference and coated sample 2 (pictures in Figure 18). It is interesting to note that the shape of the water droplets, and then of the wettability, does not change with freezing, and consequently, the coated samples display a reduced surface area in contact with the substrate if compared with the reference.

**Figure 18.** Pictures of samples acquired at −27 ◦C. Reference with water droplet before (**a**) and after (**b**) the icing at −27 ◦C; coated sample with water droplets before (**c**) and after (**d**) the icing at −27 ◦C.

#### **4. Conclusions**

In the present work, the authors studied the wettability of a new developed superhydrophobic coating applied with a layer-by-layer approach as a common aeronautical livery, that can be potentially employed as a passive anti-icing system for aeronautical applications.

It was found that the wettability of the uncoated samples decreases as the roughness increases achieving a value of about 115◦ for *R*a = 4 μm. The increase of the water contact angle due to the increased roughness alone can be ascribed to the morphology changes (MORPHOLOGICAL EFFECT). Instead, high values of contact angle, i.e., about 158◦, were achieved only after the application of the coating, regardless of the substrate roughness, ranging from 0.4 to 3 μm (CHEMICAL EFFECT), so highlighting that the chemical modifications prevail on the substrate morphology.

The coating's application reduced the SFE and WA by 99% and 94%, respectively, with respect to the reference. Finally, a preliminary test about the wettability of the developed coating at low temperatures showed that the reduced wettability of coated samples was preserved also at −27 ◦C. In this regard, future research will be addressed to characterize the developed coating in severe environmental conditions simulating real flights.

**Author Contributions:** Conceptualization, F.P.; Data curation, F.P.; Formal analysis, F.P.; Funding acquisition, L.D.P.; Investigation, F.P., D.D. and F.M.; Methodology, F.P.; Project administration, A.C. and L.D.P.; Software, G.C.; Supervision, F.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by H2020 Clean Sky 2 Framework (CS2), 807083—AIRFRAME ITD (GAM AIR 2018), Topic Identification Code: JTI-CS2-2015-CPW02-AIR-02-07, project number and acronym: 699757/SAT-AM (More Affordable Small Aircraft Manufacturing) and "The APC was funded by H2020 Clean Sky 2 Framework.

**Acknowledgments:** Giovanni Pucci (Sapienza University of Rome, Italy) is kindly acknowledged for his technical support.

**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/).

## *Article* **Assessment of Aircraft Surface Heat Exchanger Potential**

#### **Hagen Kellermann \*, Anaïs Luisa Habermann and Mirko Hornung**

Bauhaus Luftfahrt e.V., Willy-Messerschmitt Straße 1, 82024 Taufkirchen, Germany **\*** Correspondence: hagen.kellermann@bauhaus-luftfahrt.net

Received: 20 November 2019; Accepted: 12 December 2019; Published: 19 December 2019

**Abstract:** Providing sufficient cooling power for an aircraft will become increasingly challenging with the introduction of (hybrid-) electric propulsion. To avoid excessive drag from heat exchangers, the heat sink potential of the aircraft surface is evaluated in this study. Semi-empirical correlations are used to estimate aircraft surface area and heat transfer. The impact of surface heating on aircraft drag is qualitatively assessed. Locating surface heat exchangers where fully turbulent flow is present promises a decrease in aircraft drag. Surface cooling potential is investigated over a range from small regional aircraft to large wide body jets and a range of surface temperatures. Four mission points are considered: Take-off, hot day take-off, climb and cruise. The results show 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. Also, a clear trend favouring smaller aircraft with regards to the ratio of available to required cooling power is visible.

**Keywords:** aircraft thermal management; hybrid electric propulsion; surface heat exchanger

#### **1. Introduction**

Research for next generation commercial aircraft is driven by ambitious goals to reduce the aircraft's environmental impact such as the European Commission's Strategic Research and Innovation Agenda (SRIA) [1] that targets a 75% reduction in CO2 emissions by the year 2050 compared to the year 2000. A big contributor to achieve those targets is the propulsion system. Novel propulsion concepts with intercoolers [2], topping cycles [3] or bottoming cycles [4] are currently under investigation to reduce the specific fuel consumption. Another promising approach seems to be a higher electrification of the on-board systems or even the propulsion system. Examples are the more electric aircraft [5] or (hybrid-) electric propulsion systems [6]. Their electric components generate waste heat that needs to be rejected in an efficient way. Many concepts result in higher thermal loads of the systems. Conventional cooling concepts require ram air and heat exchangers, which are placed in the airflow path and thus generate drag [7]. Another option is to use existing aircraft surfaces for heat transfer from the inside of the aircraft to the ambient [8]. These structurally integrated heat exchangers may be beneficial for both weight and drag of the Thermal Management System (TMS) because no additional components such as the ram air heat exchanger are required and no components are installed in the flow path. Additionally, heat rejection to the aircraft's boundary layer may lead to drag reductions [9]. The aim of this paper is to investigate the heat sink potential of available aircraft surfaces.

Wang et al. present a good overview of the application of surface heat exchangers in aircraft up to the year 1999. In the beginning, the development of surface heat exchangers was driven by the cooling demand of piston engines with increasing power densities. In the 1920ies and 1930ies they were mainly used in racing aircraft. In some aircraft such as the "Supermarine S.6" surface heat exchangers covered surfaces of multiple components such as wings, fuselage and floats. In military aircraft, leading edge steam radiators were successfully tested. However, despite the proven thermodynamic performance

the technology was not put into practical applications due to hazards such as machine gun fire. When gas turbines started to replace piston engines in aircraft, the engine cooling problem vanished and with it surface heat exchangers. However, academic research on surface heating continued. The results indicate that heating aircraft surfaces might not only serve for heat dissipation but also as means of boundary layer control. It is commonly agreed that heat addition to a laminar boundary layer increases instabilities and therefore may lead to an earlier transition, thus increasing drag [8].

More recent studies showed a growing interest in surface heat exchangers again due to the increased cooling demand from the aforementioned technologies. Especially new engine concepts such as Ultra High Bypass Ratio Turbofans and open rotors with very compact gas generators and mechanical transmission have increased oil heat loads. Sousa et al. investigated a surface cooler with fins inside a turbo fan engine bypass as air cooled oil cooler (ACOC). Numerical calculations in combination with experiments were conducted. They showed that the surface cooler was capable of rejecting 76% of the take-off oil heat load [10]. Surface air cooled oil coolers (SACOC) are investigated by multiple EU-funded projects such as SHEFAE [11], SHEFAE 2 [12] and SACOC [13]. Sakuma et al. carried out investigations of the effects of varying SACOC geometries in the context of SHEFAE 2. They found that two 200 mm long heat exchangers could reject the same heat as one 900 mm long one while maintaining the same pressure drop allowing for area and weight optimization of the SACOC [14].

Recently, Liu et al. conducted numerical studies to describe pressure loss and heat transfer of different aircraft surface heat exchanger fin configurations including continuous, segmented and staggered fins. They found the continuous configuration to have the most advantageous heat transfer to pressure drop ratio [15]. Part of the wing surfaces were used for heat dissipation of a hybrid electric aircraft with a TMS utilizing fuel as working fluid. The results showed that steady state cooling of the electric propulsion system is possible in most operating points, however the aircraft only had 20% hybridization [16].

While there is a good amount of literature on surface coolers in aircraft applications, most investigate research questions tailored to one specific engine or aircraft or try to optimize the surface heat exchanger geometry. In contrast, this study aims to generally predict the thermodynamic potential of the aircraft surface for a range of differently sized aircraft. The goal is to quickly assess the feasibility of using a TMS with surface heat exchangers for a hybrid electric configuration.

For that purpose, a thermodynamic model of a surface heat exchanger covering existing aircraft surfaces is developed. A scalable geometric model of a tube and wing type aircraft with podded engines is derived with a semi-empirical approach. It is used to analyse the impact of aircraft size and available portions of the total surface area on the potential cooling power. Various sensitivities of the model including surface temperature, incoming radiation and component geometries are considered. In addition, drag increments resulting from non-adiabatic boundary layers are assessed.

The ambient conditions differ at each operating point. The study evaluates steady state heat transfer performance in pre-defined sets of Mach number (*Ma*), altitude (*alt*) and ISA temperature deviation (*dTISA*). They reflect typical operating conditions of commercial aircraft namely Take-Off (TO), Hot Day Take-Off (HTO), Climb (CL) and Cruise (CR), which are relevant sizing points for the TMS. The quantification of the potential of the aircraft's skin as heat sink can be used by future projects on advanced propulsion concepts to account for the total amount or a fraction of the system's waste heat removal.

#### **2. Aircraft Correlations**

The study aims to estimate the surface heat sink potential of a range of aircraft covering most of the commercial aviation market. Therefore, data for aircraft ranging from small regional aircraft up to large wide body jets are used as basis for the correlations. Besides correlations for the wetted surface area (*Awet*), the data is analysed with regard to propulsive power as it will be an indicator for the size of future hybrid electric power trains and thus the expected required cooling power (*Qrq*). Most data are obtained from Reference [17]. They provide aircraft data up to the year 2000 for aircraft from different manufactures including Airbus, Boeing, Fokker and Bombardier. Additional data especially for newer aircraft are extracted from documents provided by the manufacturers [18,19].

#### *2.1. Aircraft Component Geometries*

*Awet* of the aircraft consists of the surface areas of multiple components. This study is strictly limited to the tube and wing aircraft configuration and thus the components considered as possible locations for surface heat exchangers are:


For wing, horizontal and vertical tail the data at hand contains the exposed area (*Aexp*) that is, for a wing the area given is the base area outside the fuselage. In a first order approximation, *Aexp* is doubled to calculate *Awet*. More accurate semi-empirical methods to calculate the wetted area of bodies with an airfoil cross section for example in Reference [20] exist, however, for an initial potential assessment, it seems more reasonable to choose the simplest method possible. For the fuselage and nacelles, *Awet* is also not directly available. Instead, length and diameter (in case of the nacelle the maximum diameter) are included. *Awet* of these components is estimated by using the geometric model of a cylinder. This approach overestimates the area for the nacelles, because a cylinder with the nacelle length and the maximum diameter as constant diameter has a larger *Awet* than the actual nacelle with a variable diameter. For the fuselage, the overestimation of the lateral surface area is reduced by the fact that an open cylinder model is used but the fuselage is actually a closed body. For a quick estimation of the order of magnitude of the error from these geometric simplifications, a point validation is conducted using available data from an A320 sized aircraft model [21]. Table 1 shows the comparison between the simplified *Awet* (*Asim*) and the actual *Awet* (*Aact*) as well as the relative deviation of *Asim* from *Aact*:

**Table 1.** Comparison of simplified *Awet* with actual *Awet*.


The deviation for the total *Awet* is less than 10%. For each component, the expected direction of deviation is confirmed that is, for fuselage and nacelles the simplifications lead to an overestimation of *Awet* whereas for all the other components the methods underestimate *Awet*. The largest deviation is present for the fuselage with 15.8%. Overall, the deviations are considered acceptable for the scope of this study, because the aim is to find basic correlations among a wide range of aircraft rather than developing precise calculation methods for one specific aircraft.

#### *2.2. Surface Area Correlations*

Four possible aircraft design parameters are identified as potential variables to correlate with the total wetted surface area:


All four are defined in the conceptual design stage of an aircraft and affect the overall aircraft design. For all four variables, correlations with *Awet* are found using least squares polynomial fits. To assess the quality of each fit, the coefficient of determination (*r*2) is used. The best fits that is, the ones with the highest *r*<sup>2</sup> value for all four variables are linear fits of the log-log scaled data and are shown in Figure 1. The corresponding fits are summarized in Equation (1) with the coefficients given in Table 2.

$$
\log\_{10} A\_{\text{uvet}} = a \cdot \log\_{10} \mathbf{x} + c \tag{1}
$$

**Table 2.** Coefficients for log-log surface area fits.


**Figure 1.** *Awet* correlations with data from References [17–19].

The correlation of *Awet* with *Rdes* has a higher variance than the other three. *Awet* correlates very well with *MTOW*, *MPL* and *nmax* (*r*<sup>2</sup> > 0.95). For this study, the *MTOW* correlation is chosen because it has the highest *r*<sup>2</sup> value and *MTOW* is the most general and robust aircraft parameter to compare against. It could also be used for retro fitted cargo aircraft, which is not possible for *nmax*. However, the other correlations might be useful for a first *Awet* assessment prior to the *MTOW* calculation in the conceptual design phase. The correlations are limited to their source data that is, they may not be used outside the range of the source data. They may be used for future aircraft that is, hybrid electric aircraft if no significant change in the respective correlation is expected due to for example technology changes. From the heat transfer modelling (cf. Section 3) it becomes apparent that solely knowing the total wetted area of the aircraft (*Awet*,*tot*) is not sufficient even for the simple correlations that are used in this study. The distribution of *Awet*,*tot* among the component groups mentioned in Section 2.1 is investigated. No correlation is found with any of the x-parameters from Table 2. The share of each component of *Awet*,*tot* (*Awet*,*i*) is rather constant. Therefore, a mean is applied and the results are listed together with the standard deviation (*σ*) for each mean in Table 3.

**Table 3.** *Awet*,*i*/*Awet*,*tot* for each component.


#### *2.3. Propulsive Power*

To put the available cooling power (*Qav*) in perspective with the required cooling power (*Qrq*) an estimation of the expected waste heat is necessary. The quantity of waste heat of a future propulsion system will depend on many factors, especially the propulsive power (*Pprop*), the transmission efficiency (*ηtrans*), which includes all losses from shaft power (*Pshaf t*) to *Pprop*, the Degree of Power Hybridization (*HP*) [22] of the drive train and the overall electric efficiency (*ηec*). Calculation of the exact heat loads over the entire mission are part of a detailed iterative design process. For a first estimation, simple methods are applied to estimate the waste heat during take-off, which is likely to be one of the most critical mission points with regards to cooling requirements. Starting from the take-off thrust (*FTO*), which is available in the data set, *Qrq* is derived:

$$Q\_{rq} = (1 - \eta\_{trans}) \cdot P\_{prop} \cdot H\_P \cdot (1 - \eta\_{cc}) \tag{2}$$

with *Pprop* as:

$$P\_{prop} = F\_{TO} \cdot \upsilon\_{TO} \tag{3}$$

With *vTO* being the take-off velocity, which is calculated based on Sea Level (SL) conditions with *dTISA* = 0 and *Ma* = 0.2 as representative *MaTO*. The *FTO* values are obtained from another linear fit of the log-log scaled data over *MTOW*. The resulting fit (Equation (4)) is shown in Figure 2. It has an *r*<sup>2</sup> value of 0.983.

$$
\log\_{10} F\_{TO} = 0.913 \cdot \log\_{10} \text{MTON} + 0.895 \tag{4}
$$

**Figure 2.** *MTOW*–*FTO* correlation with data from References [17–19].

#### **3. Surface Heat Transfer**

In this section, the applied heat transfer models are described. The sensitivity of the methods to changes in geometry is tested and the impact of surface heating on drag is assessed.

#### *3.1. Modeling*

Flat plate models with uniform temperature distribution are used for heat transfer calculations for all components. Correlations for the local Nusselt number (*Nux*) from References [23–25] are applied to calculate local heat transfer coefficients (*αx*), which requires a local discretization of the geometry in flow direction. The trapezoidal shaped components (wing and tail planes) are also discretized in span wise direction to account for the different flow lengths and corresponding Reynolds numbers (*Rex*). Incoming solar radiation is accounted for in the overall heat balance by means of a material absorption coefficient and an incoming radiation power (*Prad*) on all surfaces that are exposed to the sun. Unless stated otherwise an absorption coefficient of 0.25 typical of white paint and *Prad* of 1362 W/m2, which is the constant value outside earth's atmosphere [26] are assumed. For each component, half of *Awet* is considered to be exposed to *Prad*. Those are conservative assumptions since *Prad* has a slightly lower value even at the highest flight levels than the above-mentioned value outside the atmosphere. Detailed descriptions of the used convection correlations and heat balances can be found in Reference [16]. The used 2D methods are less precise than for example 3D Computational Fluid Dynamics methods but they are sufficient for a first quantification of the surface cooling power in the conceptual design stage.

#### *3.2. Sensitivities*

The aforementioned local discretization of the heat transfer calculation depends on *Rex*. Section 2.2 focuses on correlations of the total and component wise *Awet*. However, to calculate *Rex* more knowledge of the geometry is required. For example, two fuselages with the same *Awet* have different *Rex* distributions if their slenderness ratios (Λ) differ. To account for these effects, the geometric model of the components is refined. For cylindrical components (fuselage, nacelle) the sensitivity of Λ is studied:

$$
\Lambda = l/d\tag{5}
$$

With length (*l*) and diameter (*d*). Wing and tail components are modelled as single section trapezoids with no leading edge sweep. Their geometries, specifically the span-wise chord distribution can be fully defined with the help of their *Aexp*, aspect ratio (*AR*) and taper ratio (*λ*) [20]. The following sensitivity studies are conducted around TO conditions. It is one of the most critical conditions for TMS design, because of the low air flow velocities, high ambient temperatures (*Tamb*) and large cooling demand (*Qrq*) due to maximum propulsive power. Unless otherwise specified, the values in Table 4 are assumed for the wing sensitivity studies. The values are not specific to any aircraft but generally lie inside the range of the given data. *Rex*,*<sup>c</sup>* is the critical Reynolds number and *Tsur f* the average surface temperature. In a real cooling application, the surface temperature would most likely not be uniform but have a gradient in the direction of a hot side flow underneath the surface. However, in this first approximation an average *Tsur f* is assumed for simplification.

**Table 4.** Wing sensitivity study parameters.


#### 3.2.1. Transition Location

Flat plate heat transfer correlations distinguish between laminar and turbulent flow. They rely on the knowledge of a critical location (*xc*) where transition occurs. Usually *xc* is defined by *Rex*,*<sup>c</sup>* which according to Reference [24] is between 1 × <sup>10</sup><sup>5</sup> to 3 × <sup>10</sup><sup>6</sup> depending on free stream turbulence and surface roughness. Detailed transition modelling is a complex research area and beyond the scope of this work. However, a sensitivity study with varying *Rex*,*c* and *Ma* ranging from slow taxiing *Ma* = 0.01 to representative TO *Ma* = 0.2 is conducted. The results of the theoretically available cooling power (*Qav*) as well as the relative *Qav* compared to the *Qav* at the lowest *Rex*,*<sup>c</sup>* for each *Ma* (*QRe*,*x*,*c*,*min*) are displayed in Figure 3.

**Figure 3.** Transitional Reynolds number sensitivity.

An increased *Ma* results in increased *Qav* because of the increased effects of forced convection. Shifting *Rex*,*<sup>c</sup>* to higher values, results in a decrease in *Qav*. Turbulent flows favour heat transfer more than the structured flow in laminar regions because of the increased particle mixing within the boundary layer. With increased *Rex*,*<sup>c</sup>* the portion of *Aexp* with laminar flow increases. For very low *Ma* increasing *Rex*,*<sup>c</sup>* beyond 2 × 106 results in laminar flow on the entire surface. A further increase in *Rex*,*<sup>c</sup>* has no additional effect. The transition point has a large influence on *Qav*. For 3D wings, transition is more complex than defining an *Rex*,*<sup>c</sup>* and assuming instantaneous transition. This study does not accurately account for real transition effects. The results in Section 4 assume fully turbulent

flow areas downstream the transition location. Therefore, the results of this study cannot directly be used for concepts with enhanced laminarity such as natural laminar flow (NLF) wings. Covering these advanced aerodynamic concepts is part of future work.

#### 3.2.2. Wing Aspect Ratio

*AR* is varied from 6 to 18—a range that includes all aircraft used for the correlations in Section 2 and also leaves margin for possible future aircraft with increased *AR*. The results of *Qav* as well as the relative *Qav* compared to the *Qav* at the lowest *AR* for each *Ma* (*QAR*,*min*) are displayed in Figure 4. For better comprehension of the trends in Figure 4, the effect of increasing *AR* on the local *α<sup>x</sup>* distribution for the lowest and largest *Ma* are illustrated with heat maps in Figure 5.

**Figure 4.** Wing aspect ratio sensitivity.

**Figure 5.** Local heat transfer coefficient for different Mach number–wing aspect ratio combinations.

*Qav* increases with *Ma*, because the forced convection increases, due to increasing *Re*. Depending on *Ma*, *Qav* increases or decreases with increasing *AR*. More specifically for the lowest *Ma* of 0.01, *Qav* decreases with increasing *AR*. For all other *Ma* used in the study *Qav* increases with *AR*. Two counteracting effects are the reason:


Tripling the aspect ratio results in ±8% *Qav* depending on *Ma*. The sensitivity is too weak for the expected precision of this study that aims to determine the order of magnitude of the surface cooling power. Hence, it is not regarded in the following studies.

#### 3.2.3. Wing Taper Ratio

*λ* is varied between 0.1 and 1.0. The same *Ma* range as in the previous sections is applied. Variations in *Qav* do not exceed ±2% with slight advantages for the non-tapered wings (*λ* = 1.0). The aforementioned effect of increasing flow length is positive for heat transfer of tapered wings near the wing tip but negative near the root, which leads to its equalization after integration over the entire span. As with the *AR* sensitivity, the effect is too small to be further considered in this work.

#### 3.2.4. Fuselage Slenderness Ratio

For a fuselage with *Awet* = 1600 m2, Λ is varied from 5 to 15 within the same *Ma* range as the previous sensitivity analysis. Regardless of *Ma*, the change in *Qav* from the lowest to the highest Λ value is around −8%, again due to the increasing flow length with increasing Λ. The effect is also within the expected precision of this study. For further investigations, Λ = 12 is used, which is conservative as it is one of the highest Λ values found in today's aircraft for example, for the Airbus A340-600.

#### *3.3. Drag*

For any aircraft component, which contributes to the aircraft's drag, local surface temperature can influence the aerodynamics of air passing the component surface at a certain velocity and with certain fluid characteristics. The two main occurring effects depending on the fluid's initial state are:


As the skin friction coefficient (*c <sup>f</sup>*) is significantly smaller in laminar than in turbulent flow, total skin friction drag (*Df*) of a surface can be decreased by moving the transition location downstream that is, by increasing the laminar length. During the last centuries, laminar flow control approaches have been studied intensively as a means to decrease drag. As such, surface temperature alteration can be employed to decrease the growth rate of unstable disturbances in the fluid and thus, to repress transition from laminar to turbulent flow [27]. The application of this method was shown in experiments by for example, References [27,28]. Two different approaches apply [29]:


In the two-dimensional case, Tollmien-Schlichting instabilities, which dominate the laminar boundary layer, are mitigated by cooling of the near wall boundary layer. In accordance with theory, flat plate experiments showed that the cooling of a surface leads to an increase of *Rex*,*<sup>c</sup>* and a

downstream movement of *xc* [28]. The effect is reversed when the surface is heated: the destabilizing effect of the temperature increase in the boundary layer dominates and *xc* moves upstream [27].

However, the stabilizing effect of cooling can also be utilized when a portion of the surface is heated at strategic locations. For a two-dimensional case, it was shown that heating a portion of a surface where stable laminar flow is present (preferably the leading edge) followed by a cool that is, unheated, "relaxation" surface downstream can lead to a preferable downstream movement of *xc*. The heated wall has to be situated in the region where Tollmien-Schlichting waves start to develop in the laminar boundary layer. The temperature of the near wall boundary layer is increased and when the fluid reaches the cooler wall further downstream, the temperature of the boundary layer is higher than the wall temperature. The boundary layer is cooled down and the growth rate of the unstable disturbances is decreased. The transition point moves downstream. If the surface is heated in an unstable flow region, the effect is reversed [27,29,30].

In three-dimensional airflows, however, cross-flow instabilities determine the boundary layer. Dovgal et al. showed that in this case, a temperature increase of the near wall boundary layer fosters cross-flow instabilities no matter if the whole surface or only a part of the surface is heated. The transition location moves upstream resulting in an increased *Df* [27,30]. Thus, for any three-dimensional aircraft component, localized and global surface heating in the laminar flow region facilitates laminar to turbulent transition and increases *Df* .

In contrast, when the boundary layer is fully turbulent, different mechanisms govern the flow: Heating of the near wall boundary layer reduces the turbulent *Df* . Kramer et al. conducted wind tunnel experiments and flight tests in 1999. They found that an increase of the near wall boundary layer temperature leads to a decrease of *Rex*, which in turn leads to a reduced local skin friction force [31]. For a body similar to a fuselage, they showed that the heating of the fore body leads to a higher drag reduction than the heating of the aft body, whereas the heating of the whole body has the highest drag reduction potential. The findings are supported by a numerical evaluation of the effect of heating on the turbulent boundary layer flow over slender and bluff fuselage-like bodies conducted by Lin and Ash in 1986 [32]. The following theoretical deviation of *Df* as a function of wall heating for a smooth flat plate is based on the deviation proposed by Reference [31].

The length Reynolds number is defined as:

$$Re\_{\mathbf{x}} = \frac{\rho \cdot \boldsymbol{\upsilon} \cdot \mathbf{x}}{\mu} \tag{6}$$

For *Rex* = <sup>10</sup><sup>6</sup> − 108, the turbulent *<sup>c</sup> <sup>f</sup>* for a flat plate of length *<sup>x</sup>* can be expressed with *<sup>K</sup>* = 0.036, *m* = 6 by Reference [33]:

$$c\_f = \frac{K}{\text{Re}\_x^{\frac{1}{m}}} \tag{7}$$

Total skin friction drag of a flat plate with the length *x* and total area *A* for a turbulent boundary layer is defined as [33]:

$$D\_f = c\_f \cdot \frac{1}{2} \cdot \rho \cdot v^2 \cdot A = \frac{0.036}{2} \cdot \frac{\rho}{\rho^{\frac{1}{6}}} \cdot \frac{v^2}{v^{\frac{1}{6}}} \cdot x^{\frac{1}{6}} \cdot \mu^{\frac{1}{6}} \cdot A. \tag{8}$$

Assuming a constant heated surface temperature (*Th*) along *x*, the temperature ratio of unheated air (*Tu*) and *Th* is defined as:

$$TR = \frac{T\_h}{T\_u} \tag{9}$$

*Aerospace* **2020**, *7*, 1

Applying the ideal gas law leads to *ρ* = *f*(1/*T*) and the dynamic viscosity of air can be simplified to *μ* = *f*(*T*). Thus:

$$\frac{\rho\_h}{\rho\_u} = \frac{T\_u}{T\_h} = \frac{1}{TR} \tag{10}$$

$$\frac{\mu\_h}{\mu\_u} = \frac{T\_h}{T\_u} = TR\tag{11}$$

$$\frac{\text{Re}\_{\text{x,l}}}{\text{Re}\_{\text{x,u}}} = \frac{\rho\_{\text{h}} \cdot v\_{\text{h}} \cdot \mu\_{\text{u}}}{\rho\_{\text{u}} \cdot v\_{\text{u}} \cdot \mu\_{\text{h}}} = \left(\frac{1}{TR}\right)^2\tag{12}$$

When the wall is heated (*TR* > 1), *Rex* decreases with increased temperature. In consequence, *c <sup>f</sup>* increases. However, the change in *ρ* has a larger effect on *Df* than the change in *c <sup>f</sup>* :

$$\frac{D\_{f,h}}{D\_{f,u}} = \frac{\rho\_h^{\frac{5}{6}}}{\rho\_u^{\frac{5}{6}}} \cdot \frac{\mu\_u^{\frac{1}{6}}}{\mu\_h^{\frac{5}{6}}} = \left(\frac{1}{TR}\right)^{\frac{5}{6}} \cdot TR^{\frac{1}{6}} = \left(\frac{1}{TR}\right)^{\frac{2}{3}}\tag{13}$$

and therefore if *Th* > *Tu* → *Df* ,*<sup>h</sup>* < *Df* ,*u*. The higher the wall temperature compared to the ambient temperature, the higher the drag decreasing potential. All simplified relations are depicted in Figure 6.

**Figure 6.** Theoretical impact of wall heating/cooling on a smooth flat plate turbulent boundary layer density, skin friction coefficient, skin friction drag force and boundary layer 99% thickness compared to an unheated wall. Valid for *Rex* <sup>=</sup> <sup>10</sup><sup>6</sup> <sup>−</sup> <sup>10</sup>8.

Wall heating not only has an impact on skin friction drag but also effects the pressure drag. The turbulent boundary layer velocity profile thickens, because [34]:

$$\delta = \frac{0.37x}{Re\_x^{\frac{1}{5}}}.\tag{14}$$

For a flat plate, the pressure gradient is zero at all locations. For a slender body (fuselage) or lifting surface (wing, tail planes), however, the pressure gradient varies in stream wise direction. Therefore, for a three-dimensional curved body, the heating of the wall has an effect on the (not-separated) pressure drag as shown by Lin and Ash. The heating of the wall increases the turbulent displacement thickness (*δ*∗), which in turn leads to a slight increase in pressure drag [32]. In addition, the boundary layer shape factor is increased. Thus, the adverse pressure gradient is increased, causing an earlier

flow separation [32,35]. The effect of wall heating on pressure drag is small compared to the effect on skin friction drag [32].

In summary, in regions in which the boundary layer is laminar, an increased temperature leads to an earlier transition and, thus, to an increase in total *Df* . To make use of the beneficial effect of wall heating on the turbulent drag force, the surface of aircraft components should preferably be heated only in regions in which a fully turbulent boundary layer is present. This means that for example, the fuselage nose (cockpit area) or wing leading edge (slats etc.) should not be used for heat disposal. For aircraft concepts that unite different technologies, which emit excessive heat and aim at an increased laminar flow control, detailed studies have to be conducted, compromising excessive heat disposal and drag reduction approaches.

#### **4. Surface Cooling Potential**

For the following studies the simple correlations derived in Sections 2 and 3 are combined to estimate *Qav* depending on *MTOW*. The calculated *Awet* is reduced for each component to account for more realistic cases with unusable surface area in each component. These reductions are based on observations and estimations from drawings in manufacturer's documents such as in References [18,19].

#### *4.1. Area Reduction Assumptions*

In Section 3.3 it was shown that heating surfaces underneath laminar flow has a negative effect on aircraft drag. Therefore, areas at the front of each component are avoided as locations for surface heat exchangers. Independent of the size of the aircraft, the first 4 m of the fuselage are not used because cockpit, sensors and nose landing gear bay are located here. In addition, the contraction of this part is responsible for the overestimation of *Awet* of the fuselage in Section 2. The rear 15% of the fuselage length are also not used because of the tail plane attachments, the auxiliary power unit (APU) and again the contraction that lead to an overestimation of *Awet*. For the remaining fuselage middle section, a stripe of 0.5 m width is spared on both sides to account for the windows. Another 10% is subtracted from the total middle section area to account for passenger doors, cargo doors and landing gear doors as well as sensors and air openings. The wing leading edge and trailing edge (20% chord length each) cannot be used as a heat sink due to slats, flaps and other control surfaces. Only the forward 50% of the nacelle length is used to account for possible thrust reversers installed in the back. The rear 33% of the horizontal and vertical tail plane's chord length are not used because of the installed control surfaces. In the following, all remaining surfaces are employed for heat rejection.

#### *4.2. Cooling Potential for Typical Operating Points*

*Qav* is investigated in multiple typical operating points: TO, HTO, CL and Cruise CR. The atmospheric conditions (*Ma*, *alt* and *dTISA*) of each operating point are listed in Table 5. The design space includes *MTOW* over the entire range of the database used in Section 2.1 as well as *Tsur f* ranging from 320 K to 400 K. *Tsur f* has to be lower than the maximum allowed operating temperature of electric components, which for motors can be up to 180 ◦C [36] but are significantly lower for batteries. The actual *Tsur f* depends on the installed drive train and the hot side of the cooling system. This study shows *Qav* for a wide range of *Tsur f* in Figure 7.


**Table 5.** Investigated operating points.

**Figure 7.** *Qav* in multiple operating points for aircraft equipped with surface heat exchangers.

Figure 7 can be used to estimate *Qav* of any tube and wing aircraft with known *MTOW*. For example, the A320 sized aircraft from Section 2.2 has an *MTOW* of 71,000 kg. Assuming *Tsur f* of 360 K an estimated *Qav* of approximately 250 kW in HTO—the most critical condition—results. In contrast, the same aircraft with the same *Tsur f* would be able to reject about 7 MW of heat in CR. In all operating points, *Qav* increases with *MTOW*, because *Awet* increases. The slope of *Qav* decreases with *MTOW*, because the *MTOW*–*Awet* correlation is weakly logarithmic and because aircraft with higher *MTOW* have increased flow lengths on all surfaces, which results in lower *α<sup>x</sup>* towards their rear ends (cf. Section 3.2). *Qav* also increases with increasing *Tsur f* due to the higher temperature difference to the ambient. In the HTO case, *Qav* is about five times as large for *Tsur f* = 360 K than the value corresponding with *Tsur f* = 320 K over the entire *MTOW* range. The high sensitivity is due to a relatively high ambient temperature (*Tamb*), resulting in an increase of heat transfer driving temperature difference (Δ*T*) from Δ12 K to Δ52 K (roughly factor five). For CL and CR the relative *Tsur f* sensitivity is not as strong because *Tamb* is lower. Over the entire *MTOW* range, *Qav* is about twice as large for CR and CL compared to TO. The ratio even increases when comparing CR and CL to HTO. The reasons for this large difference are the lower *Tamb* in CL and CR compared to TO and HTO as well as the higher *Ma* that increases convection. The difference in *Qav* between CL and CR is approximately 10% over the entire *MTOW* range for *Tsur f* = 320 K. The difference is less for higher *Tsur f* and hardly noticeable for the largest *Tsur f* of 400 K. CR has a lower *Tamb* than CL which results in a larger Δ*T*. The relative difference between Δ*TCL* and Δ*TCR* decreases with increasing *Tsur f* . Also, for heat transfer, the total *Tamb* is relevant and due to the increased *Ma* in CR it is not smaller by the same ratio compared to CL as the static *Tamb*. The increased flight speed should additionally result in higher *Nu* in CR but the effect is reduced by the lower *ρamb*. More elaborate studies on the dependence of forced convection on flight conditions can be found in [16].

#### *4.3. Hot Day Take-Off Performance*

Results from the previous section indicate that HTO is the condition with minimum *Qav*. Additionally, the propulsive power is usually at its maximum during TO, which means *Qrq* is at its maximum as well. For an aircraft application, the most relevant metric is the ratio of *Qav* to *Qrq* (*CQ*). The following study is conducted for HTO conditions. *Qrq* differs depending on the electric architecture and the mission profile amongst others. For the first part of this study a fully power hybridized aircraft (*HP* = 1) is assumed with estimated values for the efficiencies required by the methods shown in Section 2.3 listed in Table 6. The effect of varying *MTOW* over the entire range of the database (cf. Section 2.1) as well as *Tsur f* ranging from 320*K* to 400*K* is investigated. The results are shown in Figure 8a. For the second part of the study different *HP* values are assumed and the required *Tsur f* to achieve *CQ* = 1 that is, a *Qav* that matches *Qrq* is investigated. Figure 8b shows the results.

**Table 6.** Values for the estimation of *Qrq*.


**Figure 8.** Comparison of *Qav* and *Qrq* for hybrid electric aircraft in hot day take-off conditions. (**a**) Ratio of *Qav* to *Qrq* for different *Tsur f* . (**b**) Required *Tsur f* to achieve *CQ* = 1.

A first observation is that *Qav* and *Qrq* are within the same order of magnitude during HTO. Within the used parameter ranges values above and below unity exist for *CQ*. *CQ* is decreasing linearly with *Tsur f* and hyperbolically with *MTOW*. Smaller aircraft have a favourable *CQ*. This is mainly due to the increased flow length on all surfaces of larger aircraft but also due to the weakly logarithmic behaviour of the *MTOW*–*Awet* correlation. For the smallest considered aircraft, *CQ* ranges between 0.2 and 1.8 depending on *Tsur f* . In contrast, the *CQ* range for the largest considered aircraft is between 0.1 and 1.0. In Figure 8b the required *Tsur f* during HTO to achieve *CQ* = 1 is depicted for different *HP*. *Tsur f* grows linearly with *HP* because *Qrq* increases proportionally to *HP*. *Tsur f* grows logarithmically with *MTOW*, which is expected from Figure 8a: Smaller aircraft have an advantage over large aircraft with regards to potential cooling via existing aircraft surfaces. Taking the A320 sized aircraft from Section 2.1 (*MTOW* = 71,000 kg) as an example again with *HP* = 1, Figure 8b shows that *Tsur f* of about 370 K would be required during HTO to provide enough cooling power for the waste heat load of the drive train. Heating up the surface to an average *Tsur f* of 370 K is going to be challenging in an application with low grade waste heat potentially involving artificial measures such as vapour compression cycles to increase the temperature at which heat is rejected. However, such systems add weight and need power, which might diminish the benefits from a surface cooling system on aircraft

level. A comparison of a conventional cooling system with a surface cooling system on aircraft level will be performed in future studies. The results shown in Figure 8 are for steady state cooling in the most adverse conditions. A dynamic model might reveal that requiring steady state cooling during TO is unnecessary because thermal inertia of components and fluids can cope with temporarily high heat loads that is, *Qrq* > *Qav*. The dynamic behaviour of surface cooling systems will also be part of future work. The feasibility of using surface heat exchangers for cooling highly depends on *Qrq* and the requirements for *Tsur f* are more relaxed for *HP* < 1. For the aforementioned example reducing *HP* to 0.5 results in a 30 K decrease in required *Tsur f* to about 340 K. Thus, surface cooling might be a viable option for aircraft with lower *HP* or can be used in combination with a conventional cooling system for aircraft with large electrification to reduce heat exchanger size and drag.

#### **5. Conclusions and Outlook**

The potential of using the existing aircraft surfaces as heat sink for the waste heat of a (hybrid-) electric drive train was investigated. First, empirical correlations were derived to predict an aircraft's wetted area (*Awet*) from its maximum take-off weight (*MTOW*). The database included aircraft ranging from small regional aircraft to large twin aisle aircraft. The chosen correlation was a fit of the log-log scaled data that had a coefficient of determination (*r*2) of 0.986. To assess the ratio of available cooling power to required cooling power (*CQ*), a simple estimation of the waste heat based on take-off thrust was used. Heat transfer from wetted surfaces was modelled via flat plate correlations. To apply them, the total *Awet* was divided into five component groups: fuselage, wing, nacelles, horizontal tail and vertical tail. The mean of the relative area share was calculated for each component.

Sensitivities of the heat transfer model were studied. The flow transition had a considerable impact on the predicted heat flow. The applied methods in this work did not include accurate transition prediction. The detailed analysis of the heat transfer potential of surfaces with large laminar shares are part of future work while the results of this study may be used for concepts where turbulent flow dominates. Other sensitivities investigated were wing taper ratio and aspect ratio as well as fuselage slenderness ratio. Their impact was too small to be further considered because it was below the expected uncertainty level from the modelling simplifications. A qualitative assessment of the impact of surface heating on the aircraft's drag was performed. When heat is added to a laminar flow region an increase in skin friction drag is expected. The opposite is true for fully turbulent flow regions where heat addition reduces skin friction drag. A quantification of the expected effects is part of future work. Combining the findings for the drag with the flow transition sensitivity of the heat transfer leads to the conclusion that surface heat exchangers should only be installed in fully turbulent flow regimes to avoid a negative impact of surface heating on the aircraft aerodynamics.

Additional area reductions to account for unusable surface area for example, windows, landing gear doors and cockpit were applied and available cooling power (*Qav*) were calculated for a range of *MTOW* over the entire database and average surface temperatures (*Tsur f*) between 320 K and 400 K. *Qav* was evaluated in four operating points: Take-off (TO), Hot Day Take-off (HTO), Climb (CL) and Cruise (CR). *Qav* was largest in CR with about 7 MW for an A320 size aircraft and a medium *Tsur f* of 360 K. The most critical operating point was HTO with *Qav* of only 0.25 MW for the aforementioned aircraft and *Tsur f* .

*CQ* was calculated in HTO. The smallest aircraft showed an advantage over larger aircraft with *CQ* values ranging from 0.2 to 1.8 depending on *Tsur f* compared to 0.1 to 1.0 for the largest aircraft.

The results of this study may be used to quickly assess the feasibility of a surface cooling concept for a (hybrid-) electric aircraft. Future work will include more detailed models for surface heat transfer. Instead of assuming an average *Tsur f* , surface heat exchangers with a hot side flow will be modelled. These models that can also be used in a dynamic simulation will allow a more detailed sizing of the thermal management system. To assess heat transfer more precisely in a 3D flow field, numerical methods will be developed Those methods may go beyond the scope of a conceptual aircraft analysis and are part of more in-depth studies later in the design process.

The impact of adding heat to the boundary layer has only been qualitatively assessed in this work. Numerical methods will help to quantify the effect. Together with improved drag predictions, mass and power estimations the concept will be compared to a similar aircraft with a conventional cooling system to quantify its benefits. In addition, structural integration of surface heat exchangers may be a challenge. The concept will be evaluated with regard to maintainability.

**Author Contributions:** Conceptualization, methodology, simulation, analysis and writing of all aspects of the research except for Section 3.3, H.K.; Conceptualization, methodology, simulation, analysis and writing of Section 3.3, A.L.H.; Supervision, M.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received funding as part of SynergIE, a research project supported by the Federal Ministry for Economic Affairs and Energy in the national LuFo V program. Any opinions, findings and conclusions expressed in this document are those of the authors and do not necessarily reflect the views of the other project partners.

**Acknowledgments:** The authors would like to thank Arne Seitz for his continued support of the research, his critical review and many fruitful discussions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**


#### *Aerospace* **2020**, *7*, 1

#### **Greek Symbols**


#### **References**


© 2019 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/).
