Analysis of Ice-Shaped Surface Roughness Based on Fractal Theory

，，，，*

1.School of Systems Science and Engineering，Sun Yat-Sen University，Guangzhou 510000，P.R.China；

2.School of Mathematics，Sun Yat-Sen University，Guangzhou 510000，P.R.China；

3.Icing and Anti/de-icing Key Laboratory of China Aerodynamics Research and Development Center，Mianyang 621000，P.R.China

Abstract: The shape of ice accretion on aircraft surfaces is crucial to icing wind tunnel tests.Currently，geometrical parameters of ice，such as height，angle，and location，are used to characterise the ice shape from a 2-D perspective.However，the surface roughness of ice-shape，which is crucial to aerodynamic analysis，is always ignored.In this paper，the fractal theory is used to characterise the ice roughness，and the corresponding characterisation method is explained.An aerofoil-icing test is conducted in a large icing wind tunnel to verify the feasibility and validity of the proposed method.In the test，the icing growth information of the aerofoil surface is collected using laser line scan technology.Then，the 3-D ice shape is reconstructed using the collected data.Subsequently，the 3-D ice shape is analyzed using fractal theory，where the profile curves at different positions of the ice shape are extracted.Additionally，the corresponding fractal dimension and joint roughness characterisation are calculated to summarise the linear regression equations of the fractal dimension.Then，the data points from profile curves are extracted to simulate the fractal interpolation functions of the ice.Correlation analyses show that ice accretion on the aircraft surface exhibits fractal features，and the fractal dimension is proportional to the joint roughness characterisation，which can be used as the assessment parameter of surface roughness of ice.Consequently，the fractal interpolation simulation of the iceshape curves represent an excellent approximation of the ice accretion on aircraft surfaces.The fractal characterisation of rough surfaces provides a new approach for scientifically quantifying 3-D ice features.

Key words：roughness；fractal theory；box-counting dimension；multiple regression analysis；characterisation

0 Introduction

Aircraft icing is a common phenomenon that can cause dangers during flights.Consequently，numerous studies have focused on icing predictions，post-icing aerodynamic characteristics，and aircraft icing-safety protection.In particular，analysis of ice formation mechanisms，ice types，and the ice evolution law on aircraft surfaces is a fundamental research field.Due to the influence of atmospheric conditions and flight parameters，ice accretion on aerofoils occurs in various shapes and sizes，with complex and irregular structures.Consequently，quantitatively evaluating the structural characteristics of ice can be challenging.Some studies focus on quantitatively describing the 3-D ice shape.Currently，2-D quantification methods are used.However，these methods rely on a single parameter，such as the angle and height of the ice horn，making the characteristics of the ice shapes difficult to be fully described.

Moreover，other studies have shown that roughness greatly affects the aerodynamics of the aircraft，the convective heat transfer of the surface，and the final form of the ice after icing［1-3］.Nevertheless，when studying the effect of ice shape fidelity on aerodynamic performance［4-5］，smooth simulations can adequately represent the ice accretion by neglecting the roughness of the ice surface.Therefore，establishing a systematic description of the surface roughness of complex ice shapes for aircraft icing research is required.

Since the 1990s，the effect of ice roughness has been extensively studied by experimental and numerical methods.Some studies have documented the ice growth processes using high-speed cameras.The features are manually extracted to obtain the roughness data from the photographic images［6-9］.According to the records of early wind tunnel tests，spectral analysis can be used to calculate roughness parameters，and the results can be used for subsequent research［10-11］.

The NASA Lewis Center developed the concept of effective roughness of icing surfaces［12］，which represents the uniform distribution of equivalent particle microelements on the ice surface.Shin et al.［13-14］proposed an empirical formula for surface equivalent particle roughness with respect to liquid water content，temperature，mean droplet diameter，chord length，and other factors.However，few test results are available and even these could be biased because of subjective factors of the researchers，making previous results inadequate.In China，little research has been conducted in this field，and most focused on the relationship between roughness and convective heat transfer，assuming the roughness parameters without systematic research［15-18］.Most of the ice-shaped surfaces simulated by these assumed parameters have been uniformly structured with undulating morphologies，while ice surfaces are disordered and uneven.Consequently，the simulated results have differed substantially from the experiments.

The fractal theory was proposed in the 1970s［19］，drawing much attention in the fields of economics，material engineering，chemistry，and various other disciplines［20］.As an inimitable branch of mathematics，fractals provide a theoretical foundation to analyse irregular objectives in natural and human society.It should be noted that surface roughness measurements and characterisation have long been research topics in tribological simulations［21-22］，and their research methods have certain implications to the scientific characterisation of iceshaped surface roughness.

Consequently，this study uses the fractal theory to analyze the characteristics of ice surfaces.Firstly，a laser scanning measurement system is used to scan the ice accretion in an ice wind tunnel，process the scanned data，and reconstruct the 3-D digital ice form.Secondly，the ice growth process is examined to verify whether the ice has a fractal structure.Subsequently，the fractal dimension and joint roughness coefficients of the surface profile curve derived from the digital model are calculated，and the correlation between the fractal dimension and joint roughness coefficients is explored.Finally，the fractal dimension of ice used to construct the fractal interpolation curves is determined to reconstruct the ice-shaped surfaces.

1 Fractal Theory Analysis for Ice Joint Roughness

1.1 Fractal definition

A fractal is defined to be a set for which the Hausdorff dimension strictly exceeds the topological dimension.However，this definition excludes some sets with evident fractal characteristics.Thus，the definition is improved and a fractal is defined to be a setFholding the following properties［19］：

（1）Fhas a fine structure，i.e.，detailed at arbitrarily small scales.

（2）Fhas many irregularities that can be difficult to describe by traditional geometrical language，both globally and locally.

（3）Fhas some forms of self-similarity，approximate or statistical.

（4）Usually，the fractal dimension ofFis greater than its topological dimension.

（5）In most cases，F(xiàn)can be defined in a simple way，for example，recursively.

1.2 Fractal dimension

The fractal dimension is a crucial feature parameter of fractals，usually a non-integer.In particular，the measured relationship between a fractal and the fractal dimension must obey a power law，i.e.，asδ→0，for constantcands，then

wheresis the dimension of the fractal setF，andcthe dimensional length of the set.If we take the logarithm of both sides of Eq.（1），then

Whenδis in a suitable range，the value ofscan be obtained by estimating the slope of the function expressed in Eq.（2），i.e.，lgMδ(F) and logδ.

Various methods can be used to calculate the fractal dimension.In particular，the box-counting method is simple，easy to implement and offers good suitability.It is one of the most widely used methods to calculate the fractal dimension.It uses a series of boxes or lattices whose lengths defineδthat covers the fractal set.The different sizes of lattices intersecting the fractal set，denoted asMδ(F)，can then be counted.When the lattice lengthδ→0，the logarithmic rate of increase in the number of covered lattices is the box-counting dimension.

1.3 Characterisation of fractal dimension of ice-shaped roughness

The first step to determine the ice surface roughness is to measure the ice shapes in an icing wind tunnel.Subsequently，an image processing method can be used to calculate the box-counting dimension of the ice profile at the cross-sections of ice accretion.The detailed process is shown in Fig.1.

Fig.1 Algorithm for determining ice-shaped surface roughness

The steps can be summarised as follows：

（1）Create 2-D ice profile prints of the crosssections of ice accretion，with pixels ofN×N，N=2k，k∈N+.Then，convert the image into a binary image，and develop a 0 and 1 pixel matrix.

（2）Divide the pixel matrix into two disjoint blocks，and each block is aω×ωmatrix，whereω=2i(i=0，1，2，…，k).Then，record the number of blocks containing element “1”，denote them asMω(F)，and repeat the procedure for all blocks to obtain a sequence of (ω，Mω(F)).

（3）Obtain a diagram by fitting the data(-lgω，lgMω(F)) using the least square method.Calculate the slope of the fitting line，and the boxcounting dimension will be obtained when the slope takes a negative value.

2 Measurement of 3-D Ice Shapes in Icing Wind Tunnels

2.1 Testing equipment and method

Icing experiments for a straight wing were performed in a 3 m × 2 m icing wing tunnel in China Aerodynamics Research and Development Centre.The span and chord of the wing were 1.98 m and 2.2 m，respectively.The main material was aluminium alloy.The operating condition were as follows：H=3 km，V=50 m/s，Ts=-22.5 ℃，α=3.5°，MVD=20.0 μm，and LWC=0.8 g/m3.

During the experiments，a visible laser line scanning system was used to recover the 3-D shape development of the continuous growth of ice accretion on the model.The system included a highspeed camera，a laser，and a high-precision rotating platform.The camera imaging resolution was 4 096 pixel×3 000 pixel，and the maximum frame rate under full frame was 68 frame/s.The diameter of the high-precision rotating platform was 102 mm with a rotational angle range of 360°，at a resolution of 0.01°.Moreover，the fastest speed was approximately 20（°）/s with positioning accuracy repeatability of ≤ ±0.005°.The laser output wavelength was approximately 635 nm，and the minimum line width was approximately 1 mm.The output power could be adjusted from 0 to 120 mW.

When measuring the iced surface，the laser produced a laser sheet that was projected onto a reflector，and the reflection was re-projected onto the iced surface of the model.The reflector was then rotated using a high-precision rotary table，controlling the laser sheet scan of the iced surface.This system did not require spatial movement during measurement，so it could be installed on any side frame of the icing wind tunnel test section.The test site is shown in Fig.2.

Fig.2 Equipment and test site

2.2 Testing results

During the experiment，a single scan was performed at four different time，i.e.，5，10，15，and 20 min.The span of the scanned model was 30 mm，and the ice growth on the wing was measured.The Meshmixer software application was used to reconstruct the digital ice surfaces using 3-D points clouds，as shown in Fig.3.

Fig.3 Single scan performed at four different time

Post-processing was performed to intercept the 2-D ice profile at cross-sections of the ice accretion.For example，the profile atX=40 mm andt=10 min is shown in Fig.4.

Fig.4 Profile curves (X=40 mm, t=10 min)

2.3 Fractal characteristics of growth ice

Fig.5 shows an image from a field measurement.At early stages of the icing process，super cooled droplets define the leading edge and waterfilm development.Owing to instability，the surface of the film begins to exhibit roughness［23-24］.As the unfrozen water flows back along the chord，it breaks up into individual rivulets.Thus，the roughness elements on the water film can be observed as nucleation points.

Fig.5 Leading edge ice time sequence

The rough layer thickens over time，which makes it difficult for the droplets to be collected on the leeward sides.Thus，the icing volume of wake regions decreases.Moreover，the increase in roughness is accompanied with a decrease in the thickness of the boundary layer and the roughness elements interacting with flows，enhancing the convective heat transfer.

At a later growth stage（t=15—20 min），the distribution of roughness elements on the wing remains practically constant.Moreover，the ice is superimposed in the previous locations，and the growth rate of roughness elements is becoming appropriate.

It should be noted that the actual objective hardly satisfies the first property of fractals，and is usually restricted to a particular scale.The characteristic scale for this paper is the millimetre scale.For the ice macroscopic roughness characteristics，the ice crystals at the millimetre scale are small structures，so the first property is satisfied.Visually，ice presents a complex and disordered form，both globally and locally.By magnifying the ice profile，the observed icing regions are uncompleted，with several fine peaks.Thus，ice formation has self-similarity，and its growth process is scale-independent.At the mid-growth stage，the constant impact of droplets increases the height of the elements roughness，forming the final ice shape.The growth ice is similar in structure to the final ice shape，which is an iterative process.Subsequently，the fractal dimension of the 2-D ice profile can be calculated using the boxcounting dimension method，as shown in Fig.6（X=40 mm andt=10 min）.Its box-counting dimension is 1.173 4，which is greater than its topological dimension.Thus，each fractal property can be verified，demonstrating that the ice on the wing can be considered as a fractal set.

Fig.6 Box-counting dimension of profile curves (X=40 mm,t=10 min)

3 Analysis of Ice Roughness and Fractal Dimension

Fig.7 shows the 2-D ice profile at the cross-sections of the ice accretion.Figs.7（b—e）show the interceptions at different locations（X1=10 mm，X2=20 mm，X3=30 mm，andX4=40 mm）and different time（t1=5 min，t2=10 min，t3=15 min，andt4=20 min）.

Fig.7 Profile curves at different chord wise positions

3.1 Ice roughness and fractal dimension calculation

In actual engineering applications，the arithmetical mean deviation of the profile（Ra），the tenpoint height of irregularities（Rz），and the maximum height of the profile（Ry），are the primary evaluation parameters for surface roughness［25］.In particular，Rais the absolute arithmetic mean of the contour deviated distance |yi| within a sample range of the part surface，which can be determined as

Rzis defined as the sum of the average of five contour peak heights and the average of five contour bottom depths.It can be determined as

Ryis defined as the distance within the sample length from the peak line（ypmax）to the bottom line（yvmax）.It can be determined as

The results of roughness parameters of profile curves are summarized in Table 1.In Table 1，Dis the box-counting dimension.

Table 1 Calculated results of roughness parameters of profile curves

3.2 Correlation analysis

The relationship between the roughness and the fractal dimension was studied.In particular，a correlation analysis was performed to determine if the two variables were correlated［26］.Pearson correlation is a typical measure of the degree of correla-tion between two variables，and this measure ranges between -1 and 1.A positive value implies a positive correlation，while a negative value implies a negative correlation.The Pearson correlation coefficients betweenRa，Rz，Ryand the box-counting dimensionDwere calculated，as shown in Fig.8.

Fig.8 Pearson correlation coefficients

Fig.8 shows that all coefficients are positive.Thus，the roughness is positively correlated with the fractal dimension.The correlation coefficient betweenRaandDis between 0.6 and 0.8，implying a strong correlation.The correlation coefficients for bothRzandRyare greater than 0.8，implying a very strong correlation.Consequently，fractal analyses can be further evaluated.

3.3 Multiple regression analysis

If the fractal dimension associated with a single roughness coefficient is used to express the rough surface，there could be a problem that the fractal dimension may not correspond to the roughness.For instance，the roughness parameters of profileX=40 mm andt=10 min decrease overall，although the fractal dimension increases.Consequently，the problem requires a multi-variable coupling analysis.

Regression analysis is a set of statistical processes for estimating the relationship between a dependent variable and one or more independent variables.ConsideringRa，Rz，andRyas independent variables and the box-counting dimension as a dependent factor，a multiple regression analysis using two regression models can be performed［27-28］.

3.3.1 Regress regression method

Regress regression methods are typically used in regression analysis.These methods can remove abnormal samples during the analysis to improve the fit of the regression equation.In this paper，multiple regression analysis was performed，and the residuals for each sample are shown in Fig.9.

The vertical lines in Fig.9 indicate the confidence intervals of the residuals for each sample，which are close and involve the zero value.The tenth sample does not contain zero values and is an outlier sample，which is removed from subsequent analyses.Thus，only normal samples are considered，and the regression formula of the box-counting dimension is expressed as

R2=0.939 4，F(xiàn)=56.864 4，P=5.517 6E-07 whereR2measures the goodness of fit and takes values between 0 and 1.The closer to 1 this value is，the better the fit.In this case，R2=0.939 4 implies that the fitted equation has a high fitting degree.Additionally，F(xiàn)is a parameter that allows evaluating whether the sample is appropriate for analysis，i.e.，it can verify eligibility.The probabilityPcorresponding toFrepresents the confidence level of the result.The smaller the value，the higher the confidence in the result.A value ofP=5.517 6E-07 implies a significance level ofα=0.05.

3.3.2 Stepwise regression method

Stepwise regression is a variable selection regression analysis method to screen out insignificant samples and reserve variables that greatly affect the outcome of the fitted equation.This method establishes the optimal regression formula with a minimum of variables，as shown in Fig.10.In Fig.10，the variablesX1，X2，andX3in the diagram correspond toRa，Rz，andRy，respectively.The red text and line indicate that the variables are rejected，the blue text and line indicate that the variables are retained，i.e.，onlyRaandRzare considered for fitting the fractal dimension equation

Fig.10 Stepwise regression analysis

The results clearly indicate that the fit is good，and the equation holds.The results of the fit for both models can be summarized as follows：

（1）The regression formulas for both models fit well，demonstrating that the linear relationship between the fractal dimension and the roughness parameters holds.

（2）Both models are practical and effective.In particular，the stepwise regression model screens out some variables that do not contribute significantly.Thus，this method helps to intuitively understand the variables that have a major influence.Moreover，if multiple independent variables are present，this method can improve data availability and rationality.However，removing the independent variables can cause other problems，such as Edgenuity.

（3）Due to the small sample size collected，the applicability of the obtained equation is limited.

4 Fractal Interpolation Simulation of Ice Shape

Based on the above mathematic model，once the correlation parameters for a particular ice profile have been calculated，a small volume of the original ice data can be combined with these parameters to interpolate an ice profile.Interpolated ice shapes have many applications，such as the numerical simulation of ice accretion on an aerofoil or the product cast of ice shape for wind tunnel tests，which require scaling and simplifying of the ice shape.Moreover，the interpolated ice shapes can be highly effective in generating digital ice shapes by simply using a number of original data points.Additionally，the degree of closeness and simplification from interpolation can be flexible and controllable to satisfy different engineering requirements.

In this paper，15 points acquired from the ice shape were used as the interpolation data for the fractal interpolation with IFS and linear interpolation［29］.For example，the profile atZ=225 mm andt=10 min was considered（Fig.11）.

Fig.11 Fractal interpolation for ice shape

In general，the generated ice profiles by using the two interpolation methods maintain the basic contours of actual ice shapes.After enlargement，linear interpolation connects adjacent interpolation points with straight lines，covering the irregular characteristics of the data.By contrast，fractal interpolation employs the self-similarity principle to restore the rough surface of ice effectively.In this paper，fractal interpolation outperforms linear interpolation in parts with bigger horn curvature.

5 Conclusions

The fractal dimension was considered to characterise the surface roughness of ice shapes.Based on the complex ice shapes obtained from large icing wind tunnel tests，the box-counting dimension was calculated by extracting their surface profile.Subsequently，the relationship between the fractal dimension and the joint roughness coefficient was investigated，and the fractal interpolation method was used to generate digital ice shapes.Ice formation has fractal characteristics and the fractal dimension can be used as one of the quantitative indicators to evaluate the integrity，complexity，and irregularity of ice shapes，with the integer part indicating the topological dimension of the fractal set and the fractional portion indicating the ability of the set to fill the space.The roughness quantitative evaluation method contains real information regarding the ice shape.The combined fractal dimension describes the surface structure of ice and fully reflects the ice morphology.Fractal interpolation can simulate the ice surface and display more realistic results than traditional interpolation methods.