Development of an Analytic Convection Model for a Heated Multi-Hole Probe for Aircraft Applications

Ice accretion or icing is a well-known phenomenon that entails a risk for the correct functioning of an aircraft. One of the areas more vulnerable to icing is the air data measuring system. This paper studies the icing protection offered by a heating system installed inside a multi-hole probe. The problem is initially solved analytically, creating a tool that can be used in order to predict the heating performance depending on the flying conditions. Later, the performance of the real system is investigated with a heated five-hole probe prototype in a wind tunnel experiment. The measured results are compared with the predictions made by the analytical model. Last, the icing protection provided by the system is estimated with respect to flying altitude and speed. As a result, a prediction tool that can be used in order to make quick icing risk predictions for straight cylindrical probes is delivered. Furthermore, the study provides some understanding about how parameters like altitude and air speed affect the occurrence of ice accretion.


Introduction and Motivation
In flight control, the onboard measurement of air speed and angle of attack is necessary in order to estimate the drag and lift forces that determine the movement of an aircraft [1,2]. In the case of aircraft likely to operate under harsh conditions, such as Unmanned Aerial Vehicles (UAV), the measurement of these flight control parameters is usually performed by multi-hole probes [3,4]. The reason is their structural simplicity and robustness in comparison with other measurement techniques.
Multi-hole probes are pressure-based velocity measurement systems. They are often used in cases where the flow direction is unknown or presents large variation [5]. Multihole probes normally consist of a slender body containing a series of tubes called pressure channels. The channels extend parallel to each other until the head of the probe, where they are open to the environment. These openings are called the pressure ports. Given that pressure varies with altitude, probes used in aircraft applications are equipped with an additional series of reading ports, called the static ring, where static pressure is measured. The probes are calibrated before their usage in flight conditions and can resolve flow angles up to ±60 • with high precision [6].
Atmosphere conditions related to altitude can be determined by the U.S. Standard Atmosphere of 1976, which is an idealized, steady-state representation of the atmosphere that provides valid relations between these parameters [7]. From a pre-study, the UAV typical service ceiling is determined to be between 10,000 and 30,000 ft. According to U.S. Standard Atmosphere data [7], these altitudes correspond to typical temperatures of −4.8 • C and −44.4 • C. For these conditions, ice accretion is likely to take place.
Ice accretion, or icing, is a well-known phenomenon in aviation since it can have negative consequences for many exposed aircraft systems [8]. At temperatures below freezing, water may still remain liquid, as supercooled droplets. This state is unstable and these droplets may freeze abruptly when coming in contact with a solid surface [9]. The obstruction of one or more of the pressure ports in a probe can lead to the corruption of the system readings [10]. The normally chosen Icing Protection System (IPS) used to protect this type of system is thermo-electric protection [9]. This consists of heating the probe surface by circulating an electric current through a resistive element. The IPS can be classified into de-icing or anti-icing depending on whether they allow the ice to build up or not. If the resulting surface temperature is high enough to evaporate the impinging water, the system is said to be evaporative and corresponds to the de-icing category. If the temperature is only high enough to prevent the solidification, letting the liquid water flow over the surface driven by the flow aerodynamic forces, the system is called wet runback and belongs to the anti-icing category. Anti-icing systems require less heating power and are therefore lighter, but they cannot remove formed ice on the probe surface as de-icing systems can.
This paper characterizes the effect of adding an anti-icing system to a multi-hole probe in collaboration with the probe manufacturer Vectoflow GmbH. In the first part of the paper, the heat transfer problem between the heating system and the probe environment is studied in order to present an analytical model that makes system performance predictions depending on the aircraft flight conditions. Thereafter, the system is tested in the subsonic Wind Tunnel B (W/T-B) of the Chair of Aerodynamics and Fluid Mechanics of the Technical University of Munich (TUM-AER) in order to measure its real performance. For that purpose, a heated probe prototype is developed and instrumented in order to perform the necessary measurements during the experimentation. In the last part, the analytical model validity is discussed, the anti-icing protection provided is estimated and an outlook is given. The study structure is represented in Figure 1.

Theory and Analytical Model
In this section, the underlying analytic convection problem is mathematically explained and characterized. As mentioned, a heating system should be included in a multi-hole probe in order to protect it from icing. In order to keep ice sensitive areas protected, electrical anti-icing systems need to provide enough heat to the surfaces so that their temperature is maintained sufficiently over the freezing point [11]. In order to make estimations for this temperature with the heating system activated, the related heat transfer theory and the design of a prediction model are presented in the following.

Heat Convection Theory
The thermal problem that presents the heat exchange from an electrically heated surface to an airflow is a convection problem. The term convection is used to describe energy transfer between a surface and a fluid moving over this surface [12]. The convection heat flux q is defined as where T ∞ is the free stream temperature, T S the surface temperature and h the convection heat transfer coefficient. The total heat transfer rate across an isothermal surface of area A s can be obtained by integrating Equation (1) as With the definition of an average convection heat transfer coefficienth, Equation (2) can be rewritten as: Complexity regarding convection problems normally lies in determining the values of h orh, since these depend on three main groups: fluid properties, surface geometry and boundary layer regime [12]. As a consequence, heat convection problems can be solved analytically only for a limited number of elementary cases. For problems involving turbulent flow or separation-for example around cylinders, spheres or other curved bodies-the direct measuring of the heat transfer coefficient is still the main approach employed [13]. The experimental results are presented in the form of empirical relations that do not rely on theory but on measurements. These expressions are obtained using the method of dimensional analysis.
Dimensional analysis provides a method for computing sets of dimensionless parameters that describe a problem defined by a certain set of variables [14]. For its application, the variables affecting the phenomenon must be known beforehand [13]. Experience tells us that the convection coefficient h is a function of the velocity V, viscosity µ, density ρ and thermal conductivity k of the fluid and the characteristic dimension d: This expression may be rewritten after applying dimensional analysis in the form: where C, a and b are undefined coefficients. The three dimensionless terms in Equation (5) are the Nusselt number Nu = hd k , the Reynolds number Re = ρVd µ and the Prandtl number Pr = µc p k . Therefore, experimental testing results are usually generalized by establishing the relationship between the Nusselt and the Reynolds and the Prandtl numbers [15]: Each of the three dimensionless parameters has a physical significance. The Reynolds number represents the ratio between inertial and viscous forces. Its value can be used in order to predict the transition of the boundary layer regime [12]. The critical Reynolds number Re c at which the flow becomes turbulent depends on the geometric configuration. For example, for tube flows, it lies around Re c ≈ 2300, and for flows over a flat plate it occurs at Re c ≈ 5 × 10 5 [13]. The Nusselt number gives the ratio of the convective to conductive heat transfer across the fluid boundary. It provides a measure of the convection heat transfer occurring at the surface [13]. In the presence of turbulence, greater flow velocities in the vicinity of the solid surface promote large temperature gradients on the surface [12]. This translates into higher Nu values and consequently also to an increase of the h values. Hence, it can be concluded that the turbulent regime increases the convection heat exchange between surface and fluid.
The Prandtl number correlates three fluid properties and thus is a property of the fluid itself [13]. For air at 1 atm and 25 • C, its value is Pr = 0.707 and does not vary significantly at the conditions investigated within this work.
Furthermore, in order to make surface temperature predictions for the heated probe, it is necessary to elaborate a model for the specific geometry of the probe as discussed in Equation (6). In the following section, existing simple models are used for the derivation of the new probe model. The simple models have been derived experimentally for similar geometries and can be found in the literature.

Heat Convection Model for the Heated Probe
After the presentation of the general introduction to heat convection theory in the previous section, in this section, the heat convection model for a heated pressure probe is elaborated. The underlying problem consists of a straight probe under an axial inflow conditions. The probe has a cylindrical body and a hemispherical head shape. The airflow is assumed to be steady and incompressible. The pressure holes are neglected and a solid probe without pressure tubings is assumed in the following. The probe body is represented by a cylinder directly connected to the hemispherical head. Given that the heater is located at the front part of the probe, it is expected that the entire length of the probe will not be relevant for the formulation of the convection problem. Hence, heat is assumed to be dissipated mostly at the front part. For that reason, the cylinder length is defined as L, which represents the length of the probe that is effectively heated. Figure 2 shows the geometrical assumptions made in order to enable the definition of an analytical model.  The total amount of dissipated heat at the transformed probe q total is obtained by summing up the dissipated heat through the front hemisphere q hem and the dissipated heat through the cylinder q cyl : Hence, in order to describe the full convection problem, it is necessary to account for two different convection models: one for the front hemisphere and another for the cylinder. The two corresponding models and the resulting analytical model are presented in the following subsections.

The Front Hemisphere
The first presented model corresponds to the probe head represented by a front hemisphere of the probe diameter D. While the flow is well attached to the surface for the front half of the sphere, the adverse pressure gradient at the back half leads to separation and the formation of a wake behind the sphere. These differences cause the proportion of heat dissipated from the front and back halves to be different [16]. A schematic representation of this phenomena is shown in Figure 3. Eastop and Smith [17] present a heat convection model that differentiates between the front and back halves. This model is also discussed in Will et al. [18]. The model consists of the calculation of the average Nusselt number over a sphere. The estimation is performed by the sum of two terms: The first term accounts for the contribution of the front hemisphere and the second for the back hemisphere. For the studied case, it is desired to obtain a new equation from Equation (8) to estimate Nu D exclusively over the front hemisphere. Since the estimated value is an average, it would be sufficient to substitute the term corresponding to the back hemisphere by an additional front hemisphere term. This way, the mean value over the entire sphere corresponds to the mean over the front hemisphere. In Appendix A, additional information and an explanation concerning this assumption is discussed. The final expression gives:

The Curved Surface of a Cylinder in Axial Flow
The second model represents the cylindrical probe body. The case presented by the considered geometry is not exactly equivalent to a cylinder under an axial flow. The flow around the probe geometry is expected to go around the probe head and attach to the cylinder surface without any or small flow separation at the transition between the hemisphere and the cylindrical part. A representation of the expected flow streamlines is shown in Figure 4.  Wiberg et al. [19] introduce a series of experiments for a cylinder of diameter D in axial flow. One of the layouts consists of a cylinder with a circular disc located upstream in order to mimic a smooth flow at the cylinder ends. The circular disc has the diameter D disc = 1/3 D and is located a distance D upstream from the front of the cylinder under study. The results show that the disc upstream from the cylinder causes the flow streamlines to diverge from the cylinder axis before reaching the cylinder surface. As a consequence, a better attachment of the flow to the leading edge of the cylinder curved surface is obtained. This reduces the separation effect at the junction between the front face and curved faces of the cylinder and it becomes negligible [19]. The flow behavior for this case is depicted in Figure 5.This provides a scenario more similar to the one expected for the considered geometry that is presented in Figure 4. For this adapted configuration, Wiberg et al. elaborate an expression to characterize the heat convection over the cylinder curved surface [19]:

Resulting Analytical Model
With the determination of the two convection models for the two separated geometries, both models can be combined in order to define the resulting analytical model for the probe geometry. The model is based on fundamental convection theory, presented in Section 2.1, and defines a methodology for making T S predictions at the surface of the considered probe geometry represented in Figure 2. As specified in Equation (7), the global model is defined as the combination of two models. These two models are defined by the expressions in Equations (9) and (10). Each of them is used to determine the convected heat through each of the respective geometries using the corresponding terms for the Reynolds and Nusselt numbers in Equations (5) and (6). After adequately combining all the mentioned expressions, the result is a model that takes D, L, T ∞ , p, V and q total as inputs and returns T S as a single output. The thermophysical properties of air are obtained using the software "CoolProp" which is based on Helmholtz energy calculations after inserting the temperature and the pressure for the fluid [20]. The reference temperature used in all Reynolds and Nusselt number estimations is T ∞ . The computation proceeding is summarized in Figure 6.

Experimental Setup and Probe Assembly
In order to validate the heat convection model of the previous section, a heated probe prototype is manufactured, instrumented and tested in a wind tunnel. The agreement between the model predictions and the real performance of the system is evaluated by comparing the analytical model output to the acquired experimental data.
Like the rest of the probes manufactured by Vectoflow GmbH, the heated probe prototype is manufactured using the Powder Bed Fusion (PBF) method, being a good example of the direct tooling phase that Additive Manufacturing (AM) has experienced in the last years [21]. This phase corresponds to the application of AM in the production of finished parts. This production technique offers a higher degree of customization and the realization of more complex geometries in comparison with conventional means. Boerner and Niehuis [22] and Heckmeier et al. [23] make use of additive manufacturing advantages by employing Vectoflow probes on their studies.
The prototype consists of a straight five-hole probe with a static ring on its shaft and an axial cavity from the back of the probe with a heating element. The insertion of the heater inside the probe is represented in Figure 7. With minimisation of weight as one of the main design goals, the probe diameter is set to a feasible minimum of D = 8 mm and the probe length to L = 153 mm. A 4 mm long cavity connecting the probe surface to the the heater axial cavity is added to the design in order to facilitate the bonding of the heater by offering a way for adhesive introduction during the heater mounting process. This cavity can be observed on the printed part in Figure 8. Additionally, the pressure channels are blocked with wax at the back of the probe in order to avoid the flow of air though them during experimentation. The heating system performance is evaluated measuring temperature on the probe surface under a series of different heater power intensities and airflow conditions. The experiments with the heated five-hole probe are conducted in the W/T-B of TUM-AER. The low-speed wind tunnel, which is of Göttingen type (closed-loop), has a cross section of h · b = 1.20 m × 1.55 m. Turbulence intensity lies below Tu = 1%. The incoming free stream velocity V is monitored with a standard Prandtl probe installed at the nozzle exit, acquiring the dynamic pressure p dyn . Furthermore, a temperature probe (PT100) is installed to acquire flow temperature data T ∞ . Hence, together with the output of the Prandtl and the temperature probes, the atmospheric pressure signal p s are monitored. The power supplied to the probe heating system is controlled by regulating an external voltage source. The test configurations are depicted in Table 1. The first four configurations have an identical flow velocity V while the heater power q is increased. In the last two configurations, the heater power is maintained and the airflow speed is stepwise increased.  The temperatures on the probe surface are read by six type K thermocouples [24] mounted at different positions along the length of the probe. The axial coordinate values for each of the measurement points are given in Table 2. The probe tip is defined as the origin (z = 0 mm). Therefore, the axial coordinate z is also referred to as the distance from the probe tip. The temperature measurement T 1 is located on the probe head. Then, T 2 , T 3 and T 4 are located between the probe head and static ring. Last, T 5 and T 6 are located after the static ring, with T 5 located very close to it. A better understanding of the exact position of the temperature measurement locations is represented in Figure 9. The measuring points are not located over a common axial plane over the surface, since temperatures are expected to show independence with the azimuthal coordinate θ due to axial symmetry. During the test, the probe is positioned in the wind tunnel test area aligned with the airflow. Figure 10 shows the final setup. The probe is located centered near the wind tunnel nozzle. For each test configuration, the wind tunnel is turned on until the desired air speed is reached according to the read dynamic pressure. After reaching thermal equilibrium, the airflow temperature as well as the read temperatures by the thermocouples on the prototype surface are acquired for each configuration.

Results
In this section, the acquired temperature test results are presented first. Then they are compared with predictions made by the developed analytical heat convection model in order to evaluate the agreement with the real system behavior. Finally, the expected system behavior under real application conditions is estimated. The result of this last step is the generation of icing prediction graphs with respect to flight altitude and speed.

Temperature Measurements
The final test results are presented in Table 3. According to an uncertainty evaluation of the measurement data, all measurements show a deviation lower than ±1 • C with a 95% confidence level. The test data are depicted in Figures 11 and 12 with data points, while in addition spline curve fits are added. The resulting temperature profiles are represented with respect to the axial coordinate z. For all cases, temperature increases from the probe head to the heater, reaching a maximum, and then decreases as the distance from the probe tip is further increased. Figure 11 shows how temperatures increase as q increases, while in Figure 12 temperature trends decrease as V is increased.

Comparison to the Analytical Model
The test results are used in order to evaluate the validity of the predictions made by the analytical model. The comparison between the analytical model output and the experimental test results is done by defining a representative temperature T test that approximates the profile mean temperature over the probe length considered by the analytical model. This length is defined as L = 5D = 40 mm for all test configurations. This length was chosen due to the positions of measurement points and represents the area most influenced by the heating element. Figure 13 represents this length over the probe geometry and the position of the available test reading points in order to determine the most adequate way to define T test . Since T 2 , T 3 and T 4 present an acceptably even distribution over the analytical model length, the value T test is defined as the mean temperature averaged with these three values. The comparison between T analytical and T test is given in Table 4. The formulas for the calculation of the error ∆T and the relative error δT with respect to T test are given in Equations (11) and (12). The relative error is computed with respect to the temperature difference to set T ∞ as the reference. The comparison between T analytical and T test is also represented in Figure 14, where the function T analytical = T test is represented by a line in black in order to show the agreement between the test results and the model.  From the comparison shown in Figure 14, it can be concluded that the analytical model and the test results present a good agreement since all points fall very close to the T analytical = T test line. The results shown in Table 4 show that the relative errors δT vary from −4.4% to +3.0%.

Evaluation of the Heating System Anti-Icing Capability
The air temperature and static pressure inside the TUM-AER wind tunnel test section cannot be adjusted to flight conditions. Hence, in order to translate the test results to typical flight elevation atmospheric conditions, expressions are found in relation to the dimensional analysis problem presented in Section 2.1. For this case, the selected dimensionless number to be conserved between different scenarios is the Reynolds number Re D . By preserving Reynolds similarity, the determination of equivalent air speeds is performed as: with air properties ν 1 , k 1 for scenario 1 and ν 2 , k 2 for scenario 2. According to the developed analytical model, Nu D conservation is a direct consequence of the Re D conservation. Furthermore, the dissipated heat at one scenario or another is independent of the air properties; this is also a quantity that remains equal between scenarios. This new consideration results in the following equations which allow the translation of equivalent temperatures between scenarios.
It is desired to study the efficacy of the heating system of the prototype as an anti-icing system. To do this, the temperatures at the probe head and static ring are predicted by using the data measured during the experimentation. The output of the model should be a graph that represents the predicted temperatures depending on the flying conditions, that is, the flight speed and the altitude, which can be related to certain pressure and temperature conditions according to the Standard Atmosphere data. The anti-icing evaluation graph is conceived as the maximum performance that the system can deliver. Therefore, only the data for the system working at maximum capacity is used. This means that from the results presented at the introduction of Section 4, only those from configuration 4, 5 and 6 are considered here. The head and static ring temperatures, T 1 and T 5 , from these configurations are represented in Table 5 in • C. The test data are translated to the Standard Atmosphere cases of 0, 1.5, 3, 4.5 and 9 km elevation using Equations (15) and (19). This is an extrapolation of the data taken in wind tunnel experiments and the results are represented in Figure 15a,b in • C.
It can be observed that, as the respective altitude is increased, the translated air speeds increase, displacing the data points to the right. For the 9 km case, the displacement is high enough that the data point matching configuration 6 falls outside the air speed range considered by the graph. For the probe head, this point is located at V = 87.8 m/s and has a value of −14.0 • C. For the static ring, the point falls at the same speed and its value is −5.6 • C.
Given the good agreement observed between the analytical model predictions and the test data, T analytical is the selected predictor in order to interpolate and extrapolate the head and static ring data shown in Figure 15a,b. For the probe head as well as for the probe static ring, a linear regression study is performed for each of the Standard Atmosphere cases based on the translated data. All these regression models together form the global predictive model. The model returns the completed anti-icing evaluation graphs which are presented in Figure 16a,b.
These graphs include the translated test data together with the predicted temperatures by the regression model. The predictions are displayed as contour plots, where the isotherms are drawn as black curves with identifying temperature labels in • C. The colored areas represent qualitative degrees of icing risk. Risk is considered to be high when the predicted temperature is below 0 • C and moderate when it is below +20 • C. The data are extrapolated to the left and right of the given data points to enable creation of the full contour line plot. However, this increase in uncertainty does not create an issue as the majority of regions to the right of the points lie in temperature ranges below 20 • C (moderate risk) and those to the left of the points are quite safe ranges of operation.

Discussion and Outlook
The heat transfer problem of a heated probe under axial flow was analytically and experimentally investigated. The analytical model is based on fundamental heat convection theory and uses the combination of two Nusselt number Nu determination models in order to estimate an average value on the probe surface. Furthermore, the system performance is investigated in Wind Tunnel B of the Chair of Aerodynamics and Fluidmechanics of the Technical University of Munich by performing temperature measurements on the probe surface for six different configuration sets of input parameters. The measurements are represented as temperature profiles along the probe surface. The measured temperatures are compared to the predictions made by the analytical model. The comparison with the analytical model shows relative errors δT between −4.4% to +3.0% for the different test configurations. Lastly, the capability of the probe heating system for maintaining the probe's protection from icing is evaluated. This is done by predicting the temperature at the areas of the probe where the pressure ports are located-the probe head and the static ring. The test data are expanded to a series of elevation cases that consider air properties according to Standard Atmosphere data. For each case, linear regression is used to estimate the expressions that relate the probe head and static ring temperatures with the analytical model output. The resulting anti-icing evaluation graphs are used to determine the probe icing risk when operating at certain elevations and speeds.
A future investigation step should include the repetition of the presented tests under icing conditions in an icing wind tunnel. By directly studying the performance of the system in icing conditions it would be possible to check the validity of the test results translation performed during the analysis presented in this paper. Coming back to the empirically derived Equation (8), a formulation for the behavior of the entire sphere is given. The two terms represent the two trends for the change in Reynolds number: The front hemisphere follows the Re 0.50 trend, and the back part the Re 0.92 trend. Hence, considering just the first term in Equation (8) is equivalent to substituting the back half term by zero. By doing so, we do not have half a sphere, but an entire sphere that is not convecting any heat to the environment through its back half. In order to actually get the appropriate value for the average Nusselt number for the front half, it is necessary to consider a sphere with two front halves. Therefore, to adequately consider the front half, we should take the same 0.42 · Re 0.50 expression for an imaginary back hemisphere.
According to the information presented so far, an assumed expression for the average Nusselt number of a front hemisphere should not only provide a higher value than the average value for the whole sphere for low Reynolds numbers, but also should return a similar value for Re ∼ 70, 000. In Figure A1, the expressions presented in Equations (8) and (9) are plotted in order to prove the completion of these conditions after the performed assumption. The represented data validates the assumed expression for Nu f ront given that the average Nusselt numbers have the same value at Re = 90,000, which agrees well with the expected Reynolds number Re ∼ 70,000.