Detailed atmospheric ice accretion surface measurement using micro‐computed tomography

Surfaces exposed to atmospheric cold temperature and humid environments are prone to ice accretion. Airplanes, electrical power transmission cables, and wind turbines are typical examples for which icing has to be considered. The measurement of the resulting ice shapes is a challenging process. While macroscopic characteristics of the ice geometry can be observed using photography and optical scanning techniques, microscopic measurements are difficult to conduct because grooved surface partially occludes the geometry of chasms. To overcome this optical inaccessibility, we propose a method to carry out detailed high‐resolution measurements of the accretion surface with micro‐computed tomography. This approach provides a unique visualization of the empty spaces in the feather region. The information obtained by this technique can improve the understanding of ice accretion physics and its computational modeling.


| INTRODUCTION
Atmospheric ice accretion occurs on several engineering applications. Electrical power lines suffer from additional weight and loads due to ice accumulation (Farzaneh, 2008). Simulations show that the loss of performance of wind turbines can be as large as 25% when exposed to icing conditions (Yirtici et al., 2019). Aircraft surfaces are also sensitive to atmospheric icing events (Gent et al., 2000), which may lead to major safety risks.
The shape and appearance of the ice accretion on engineering components are highly dependent on the atmospheric conditions. The environment temperature, duration of the event, and cloud characteristics in terms of water concentration (Liquid Water Content, LWC), and the size of water droplets (Median Volume Diameter, MVD) have a strong influence on the characteristics of the accretion formed during an icing encounter. Cloud measurements have been conducted to evaluate their composition during in-flight encounters. Standard envelopes of LWC-MVD for aircraft icing conditions are provided by EASA (EASA, 2016). Similar standards describing icing conditions for wind turbines, where the blades can suffer from considerable accretion, can also be found.
Atmospheric ice accumulation is usually classified into two common categories (Gent et al., 2000). All the parameters mentioned above influence the type of ice, but its appearance is particularly sensitive to temperature and water concentration in the cloud. Accretion taking place at temperatures close to 0 C and higher LWC values is known as glaze ice. This type of ice has a transparent appearance and grows in a freezing process in which a water film forms over an exposed surface. Colder temperatures and lower LWC values lead to the formation of rime ice, where water droplets freeze almost instantly after impact on a dry surface. This type of ice has a milky-white opaque appearance due to the number of air pockets entrapped between the frozen droplets.
Ice accumulation is a complex phenomenon resulting in shapes with features in several length scales. To assess the characteristics of the shape, features in a broad range of scales (from 10 −2 to 10 −5 m) have to be measured. Attributes like ice layer thickness and ice accumulation limits can be considered as macroscopic features of the accretion since they are usually in the range of several millimeters (McClain et al., 2017). On the other hand, surface roughness (Liu et al., 2020b) and ice feathers chasms in the inception phase, as well as pore sizes (Crabeck et al., 2016), are features which can be as small as tenths of millimeters (<0.1 mm) and can be considered as microscopic features. Ice porosity and its internal structure have been analyzed in such scales (down to a few μm) using micro-computed tomography (Crabeck et al., 2016). This technique can be applied to casts to produce a high-resolution scan of the ice surface.
Icing wind tunnel operators have developed different ways of measuring the ice thickness. Cut sections and pencil tracing have been the most used methods to capture ice characteristics (Lee et al., 2014). More complex approaches like ultrasonic-based techniques have proved to be useful (Liu et al., 2017) to determine the type of ice and its thickness, depending on the acoustic attenuation coefficient of the sample. Regular photography has been used to quantify roughness elements (Anderson et al., 2003) and their diameter, while photogrammetry has been proposed to measure the roughness in frost due to cold-soaked conditions (Miyauchi et al., 2019).
High-speed imaging combined with a digital image projection technique has been implemented to assess ice shapes on power cables (Veerakumar et al., 2020) and a wind turbine airfoil model (Gao et al., 2019) with deviations of only 0.15 mm during the calibration process. A similar approach has been used to characterize not only ice but also water film thickness on an airfoil during ice accretion (Liu et al., 2020a). Optical (Gong and Bansmer, 2015b;Broeren et al., 2018) and mid-IR laser scanning methods (Gong and Bansmer, 2015a) have been developed to improve the accuracy of the measurements. A comparison of different scanning methods is provided by (Neubauer et al., 2019). All the scanning methods tested can properly scan surfaces created by overlaying different sandpapers on top of a 2D ice shape.
The techniques mentioned so far have improved the capability to assess 3D ice geometries when compared to pencil tracing. Nevertheless, limitations of these techniques are ice accretions presenting feather features close to each other or deep gaps between scallops features. Occlusion problems cannot be overcome by these approaches. As mentioned by Lee et al. (2014), computed tomography can reconstruct gaps between ice structures, giving insights usually ignored by other means. We present here a further step to overcoming this challenge by showing that micro-computed tomography can be used to create a detailed reconstruction of the ice accretion surface.

| METHODOLOGY
An airfoil model (Velandia and Bansmer, 2019) was mounted in the test section of the Braunschweig Icing Wind Tunnel (Bansmer et al., 2018), where icing encounters can be reproduced. Two different ice types were tested reproducing glaze ice and rime ice. The accretion obtained for these cases was measured and reconstructed. We created a high-fidelity mold of the accretion to preserve the characteristics of ice and scan it using microcomputed tomography. The mold and cast procedures have proven to be a good way to document ice accretion geometries (Broeren et al., 2018).
We conducted measurements for three cases: A, B, and C. All tests were performed at 50 m/s and a LWC of 0.43 g/m 3 . For case A, rime ice formation, the temperature of the tunnel was set to −15 C, the MVD was 40 μm and the duration of the test was 120 s. Cases B and C exhibit glaze ice formation. The temperature of the tunnel in these cases was −5 C and the MVD was 50 μm. The duration of the test was 120 s for case B and 40 s for case C. Considering these boundary conditions, it is possible to estimate the freezing fraction of the accretion process at the stagnation point. Following the procedure suggested by Bond and Anderson (2004), the freezing fraction at the stagnation point of the model was evaluated as 1 for case A and 0.404 for cases B and C.

| Mold and cast
A mold of the accretion was produced by extracting the model from the test section, creating a silicone mold, and using this mold to make a polyurethane cast.
• Iced model extraction. The airfoil was taken to a climate chamber, where the temperature was set to −10 C. This way, ambient conditions were controlled during the manipulation of the accretion and ice melting was avoided. The model was placed on its side so that the accretion was not in contact with any surface. • Creation of silicone mold. A cylindrical container was placed on the leading edge of the model. The location of the cylindrical shape can be inferred from the blue dummy shown in Figure 1a. The height of the container was aligned with the spanwise direction of the model. Low viscosity two-component platinum silicone was used to generate the molds. To avoid air pockets due to air entrapment in the final mold, the components were placed in a desiccator and degassed for 2 min using a vacuum pump. The components were carefully mixed, and the degassing process was repeated once the mixture was homogeneous. The silicone was poured into the cylindrical container and left to cure for 16 hr. To achieve good sealing between the cylindrical container and the model surface, ice was melted near to the contact zones between the container and the airfoil, as can be seen in the region marked with an asterisk in Figure 1c,d. A portion of the ice was also melted in the highest part of the sample to ensure a clean profile contour in the silicone mold. This clean section was used to reconstruct the original shape of the model and compare it to the iced one. The negative of an ice layer on the silicone mold can be seen in Figure 1b. One of the regions melted to improve the container's seal is shown in the bottom part of the black rectangle. • Creation of polyurethane cast. A polyurethane cast was made around the silicone model. A cylinder covering the model was produced using a two-component liquid polyurethane mixture. The process used to remove air from the silicone was repeated for the polyurethane. After curing, the obtained cast is a highfidelity copy of the accretion. The mold and cast were not separated from each other. By keeping them together, we ensured the structural stability of the silicone (thus reducing deformation) and the high contrast of the interface between the materials for the scanning step. Furthermore, due to the size of some of the ice feathers, removing the mold from the cast could lead to the breaking of the smallest structures of the original shape.
Two major drawbacks can arise from the outlined mold and cast procedure. First, ice surfaces, when molded, can suffer from shrinkage so that the cast has not the same geometry as the original shape. This effect will be further discussed in Section 3.1. Second, air can be entrapped in the mold and cast procedure, despite the degassing conducted. This can affect the reconstruction wind tunnel to a high-resolution digitized ice shape for case a. (a) Shows the ice accretion and (b) shows a sketched of a dummy of the negative silicone mold. The 2D image obtained from the micro-computed tomography and the corresponding region of interest (ROI) to be analyzed is shown in (c). The resulting high-resolution slice is presented in (d). The region marked with an asterisk (*) was melted to make a good seal between the profile and the silicone container. The dotted red line represents the approximate location of the clean profile of the ice geometry. The identification of air, in our samples, is addressed in Sections 2.3 and 3.2.

| Micro-computed tomography
The samples were scanned using high-resolution micro-computed tomography (MicroXCT-400, Zeiss [Xradia], Germany). The non-destructive scanning allowed us to generate a volumetric view of the internal sample. A schematic of the principle of work of this technique can be seen in Figure 2. As a nondestructive method, the XMT system leads to a volumetric view of the internal sample structure. In contrast to conventional scanners, XMT systems provide an optical magnification. The resolution achieved is determined by the position of the source and the detector and the optical lens used. A sample located on a stage is exposed to an X-ray beam emitted by the source. The dependence between the material absorption, its density, and the shape of the sample results in a radiography image. A scintillator converts the radiography into a spectrum of visible light. After the magnification, a CCD-Detector takes a picture of the projected image as shown in Figure 2. Multiple reconstructions are used to build a volumetric representation of the internal structure of the sample. This volume is then sliced into several 2D images. Each image shows one horizontal slice with a height determined by the resolution of the scanning process. The 3D internal structure can be analyzed by stacking all the 2D images obtained from the measurement. The ice accretion surface is represented by the interface between the silicone and the polyurethane. The sample rotated on a stage and several X-ray images (source voltage 100 kV and source current 100 μA) were taken from several angles (rotation step 0.0693 , exposure time 8,000 ms). These images were reconstructed using a back filtered projection algorithm and corrected for beam hardening and center shift using the software XM Reconstructor (Xradia). A 0.39X objective was used for the scanning reaching a final voxel size of 27 3 μm 3 . The resolution was the same in all directions, which means that each horizontal slice is 27 μm thick. The total volume of the samples was about 165 cm 3 for all cases. Nevertheless, the region of interest analyzed is smaller for each sample (around 35 cm 3 ).

| Image processing
The result of the scanning process was a stack of images showing the accretion contour. The images were cropped to a region of interest containing only the interface between the two materials. This step was intended to reduce the data for processing and remove areas that are not relevant to the reconstruction of the accretion. Figure 1c shows the region of interest (ROI) as a white rectangle and demonstrates that the cropping conducted does not influence the measurement of the ice layer.
Two different image processing approaches were conducted and compared to obtain the interface between the two materials ( Figure 3): • Method I. The region of interest was filtered using a Gaussian smoothing kernel (kernel size is 9 × 9) with a standard deviation of 1.65. The images were binarized using a global threshold (set at 0.5). The results of the binarization process were fed to a Canny detection algorithm (Canny, 1986), which is based on gradient intensity, to obtain the accretion contour. • Method II. The region of interest is denoised using a 3 level Daubechies (third order) wavelet decomposition approach, which has proved to remove noise from CT images while preserving the main structures (Mohammadi and Leventouri, 2019). An adaptive threshold (Bradley and Roth, 2007) was used to binarize the image. This method overcomes possible inaccuracies related to intensity variations from beam hardening or ring artifacts. Again, the same edge detection algorithm based on gradient intensity was applied to the binarized image to obtain the interface between materials.
This step of the method also allowed to identify air inclusions. Figure 3f shows the interface between silicone and polyurethane for one sample where no degassing was conducted, as well as a plot of the image intensity level on the section L2-L2', crossing an air inclusion. These inclusions can be detected because of their lower intensity in the image. The absence of dark regions in the F I G U R E 2 XMT principle (adapted from Landis and Keane (2010)) interface between the two materials confirms that no air was entrapped during the mold and cast procedure.

| Optical scanning
Optical scanning was used to set a reference measurement for case C. The cast was placed in a GOM ATOS ScanBox. 1 This technology is based on a triangulation method (Galanulis et al., 2020). A blue light striped pattern is projected on the object to be measured. Two cameras scan the geometry from two different perspectives. The deformation of the stripped pattern and the position of the cameras allow a detailed 3D reconstruction of the surface scanned. This type of measurements is extensively used in the quality control and inspection process. The scan produced by this technique was compared to the results of case C.

| RESULTS
Before analyzing the scanned surfaces, it is necessary to address the shrinkage of ice structures first. A comparison of the proposed method with another established technique is also presented here.

| Shrinkage estimation
Simplified ice structures were created by freezing water accumulation to estimate the impact of shrinkage on the size of the final cast. Surface conditions in the wind tunnel were simulated using a metallic plate as a substrate. The compared ice structures can be seen in Figure 4. The horizontal and vertical extent of each structure was measured. The difference between the ice and cast is plotted in Figure 4e. Ice structures were about 0.1 mm larger (4.47% on average) than their respective cast. The shrinkage resulting from mold and cast of larger ice accretion has been previously reported (Lee et al., 2014).
This shrinkage can be explained by two main effects. First, ice can sublimate before the silicone is poured. We estimated the shrinkage of an ice hemisphere of 1 mm of diameter using Stefan's law (Holman, 2009) and assuming ice has a diffusion coefficient of D = 6 × 10 −6 m 2 /s (Douglas and Mellon, 2019) in air. Under these conditions, the hemisphere would shrink 0.012 mm per hour. In half an hour, which is the time of preparation and pouring of the silicone mixture, both horizontal and vertical lengths would have experienced a reduction of 0.012 mm. Second, shrinkage is one of the most typical mold and cast defects. Assuming both effects account for the mean shrinkage of these samples, cast shrinkage F I G U R E 3 Image processing methods and comparison. A slice of the scan for case a is shown in (a). (b) Presents the image intensity (gray line) on the L1-L1' section in (a) as well as the result of the filtering process using methods I (black line) and II (red line).
The filtered images are shown in (d) and (e), respectively. The edge detection between polyurethane (p) and silicone (s) for methods I (black line) and II (red line) can be observed in (c). (f) Exhibits an example of air entrapment for a case without degassing between the two materials on the section L2-L2' was estimated to be 0.088 mm on photographed structures. As a consequence, the overall uncertainty of the proposed measurement technique can be estimated. The shrinkage of the silicone mold, as the predominant phenomenon, equally distributes over the entire mold. Therefore, the topology of the surface is preserved except for a mean displacement in the order of 0.1 mm, which in our case results in a decreased ice thickness measurement. The alteration of the surface roughness is mostly attributed to the sublimation of ice, which is in the order of 0.01 mm.

| Surface reconstruction from microcomputed tomography images
To begin the reconstruction process, the region of interest was inspected for the presence of air. The amount of polyurethane and silicone was estimated by marking each pixel according to its image intensity value. Thresholds for both materials were identified, and any pixel below these intensities was marked as air. Figure 3f exhibits how air entrapment can be identified. Case A, in the region of interest, was composed of 64.1% of polyurethane and 35.9% of silicone while for case B, these percentages were 66.5 and 33.5%, respectively, showing that no air entrapment was found.
Since no air entrapment was identified in the samples, the ice surfaces were well represented by the interface between the two materials. Two methods were tested and compared to estimate the location of the interface (see Figure 3). First, the original images ( Figure 3a) were filtered to reduce the noise but maintained the location of the interface appropriately. The result can be seen in Figure 3d,e.
This effect can be better examined by extracting the image intensity in the section L1-L1' (Figure 3b), which crosses a feather structure in case A. The intensity image in L1-L1' is noisy, which hinders the identification of the interface. Method I shows, even after smoothing, some oscillation of the intensity values, while method II presents a smoother profile where edge identification based on intensity gradient can be performed. The detection of this interface is also affected by the type of threshold used (global for method I and adaptive for method II). Figure 3c shows the location where both methods identify different materials. Method I estimates the interface at locations 39, 64, and 90 px, while method II estimates the interface at locations 43, 61, and 95 px. The mean difference of 3 px (approx. 0.08 mm) in the interface location provides an assessment of the sensitivity of the results to the processing method.
All cases were reconstructed using the method II. Figure 5a,b shows the ice surface reconstruction for cases A and B. Case A presents an accretion where the water droplets freeze upon the impact, and a lot of feather structures appear on the upper side of the profile. Case B shows an accretion with different rough zones and a smooth zone near the stagnation line of the profile (Velandia and Bansmer, 2019). Both accretions exhibit ice features that are difficult to measure. The height and spatial distribution of the roughness elements or the void pockets between feathers are challenging elements to evaluate in these geometries. The use of micro-computed tomography allowed to analyze some of these aspects, as shown in the following section.

| Insights into feather formation
A detailed patch of the feather zone in case A is extracted in Figure 5d. The contour is colored by height, which was measured perpendicularly to the profile. Our measurements have shown how ice feathers are independent structures at the beginning of the formation when they are isolated from the main ice accretion (near the position s = L3 0 ). Afterward, bridges are formed between these structures, creating bigger clusters with void regions (near the position s = L3). These void spaces between ice feathers can be well estimated by our approach. Furthermore, the geometric characterization of the feathers is possible, estimating the growth angle to α = 60 , and the maximum height was h = 3.6 mm. These insights are useful to modelers and physicists studying ice formations.

| Comparison with optical scanning technique
To assess how micro-computed tomography performs compared to optical scanning, the cast for case C was scanned using micro-computed tomography and an optical scanning device (GOM Atos Scanbox). We compared both scans and provide the signed distance (δ) between the reconstructions in Figure 5c. The contour shown corresponds to the reconstruction from the micro-computed tomography data. The coloring of each point is the signed distance to the optical scanning result. The differences lie in a range of ±0.2 mm. The contour shows that most of the points are beneath the plotted reconstruction, which can be explained by the smoothing of the surface due to the scanning technique. The feather region shows deep blue spots, representing the occluded regions that the optical scan was not able to measure. There is one redcolored ice element on the top border of the cast. This structure was damaged when separating the cast from the mold and was therefore not present during the second scanning process. This proves the sensitiveness of the casts to handling.

| CONCLUSIONS
We have presented a measurement procedure to conduct a detailed characterization of atmospheric ice accretion. The mold and cast approach proved to be able to capture and preserve detailed information of the ice geometry.
Minor shrinkage was observed when producing a cast of simplified structures. The casts were scanned by microcomputed tomography, and two methods were compared to recognize the interface between the mold and cast on the scans. The images were denoised using a wavelet approach and binarized by an adaptive threshold. A gradient-based edge recognition was applied to the binarized image. The generated reconstruction compared well to an optical scan. The results obtained by the proposed mold and cast and the scanning method provide geometric information of structures that have been difficult to characterize by other techniques so far. This information is important to understand and build numerical models of the ice accretion process.

F I G U R E 5 A 3D reconstruction of the interface between mold
and cast representing case a (a), case B (b) and case C (c). Case C (c) is colored by the signed distance between the micro-computed tomography and the optical scan measurements (δ in mm). A detailed high-resolution feather structures in the section L3-L3' of case a is shown in (d)