Multiphysics Based Design Study of an Atmospheric Icing Sensor

A design study of an atmospheric icing sensor is presented in this article. The proposed design of the sensor is aimed to detect an atmospheric icing event, to determine ice type as well as to measure icing load, icing rate and melting rate. This sensor is constantly slowly rotating cylinder with four fins which measure atmospheric icing load and icing rate using rotational load measurement technique. By replacing the fins with capacitance measurement cards, this sensor will be able to confirm an icing event, determing icing type and measure melting rate. Combination of multiphysics tools are employed to optimize the various design parameters of the sensor. Finally a working model was manufactured to check the basic functionality of the sensor in the laboratory.

Usually, active anti/deicing systems are integrated with ice sensors.Important parameters to know about icing events are rate of icing, icing load, icing type and melting rate.The presently commercially available anti-icing and de-icing systems have not yet been proven reliable [8] and one reason of their failure could be associated with the inaccurate and incomplete information of an icing event.A comparison of different icing sensors was performed in order to develop a new ice measurement system to optimize the performance of Wind Turbines.This study proposes to combine different sensing methodologies to quantify a larger number of icing parameters, to get a complete picture [9].
Based upon the state of the art review [2], it was found useful to design an icing sensor which should be capable to detect an icing event, determine icing type, measure icing load, icing rate and melting rate.The following design constraints were identified as desirable in the proposed design: a) The sensor preferably slowly rotate around its own axis, to ensure uniform ice accretion characteristics and uniform loading on the drive motor, at the same time this would prevent icebound clogging of the device.b) The sensor will utilise multiple physical phenomena for measurement of various atmospheric icing parameters.As such, the proposed device will be a complete atmospheric icing measuring device.
Keeping in view the above mentioned aim and design constrains, the following objectives were defined: i.To measure atmospheric icing load and icing rate using the rotational load measurement technique.ii.To confirm icing event, determine icing type and measure melting rate using the capacitance measurement technique.
It was desired that the sensor should have constant slow rotational motion round its own axis in order to allow reasonably uniform ice accretion on its surface and to measure icing load and icing rate using rotational load measurement technique.To meet this initially a parametric design study was carried out by selecting potential geometric candidates for this new generation of icing sensor.The potential geometric candidates were then compared at different operating and geomtric constrains using CFD tools.Based upon the CFD results a working prototype of the sensor was manufactured, to inspect the basic functionality of the sensor in the laboratory environment.
Based upon the ISO 12494 [3] recommendations, a circular geometry was selected as a potential geometric candidate to meet the design constrains and design objectives.However in order to make sure that cylinder is relatively better than other symmetric geometric shapes exposed to same atmospheric icing environmental conditions, two other symmetric geometries hexagon and square were also selected for further analysis.In order to quantify a larger number of icing parameters, it was required that the sensor should utilize another physical phenomena e.g.capacitance for confirming an icing event, determining icing type and measure melting rate.For capacitive measuring technique, it was found reasonable that the sensor requires at least a pair of fins which should face each other.The fins should be installed at the natural corners of the geometry for improved structural and aerodynamic stability.Hence based upon these arguments, following geometric candidates were shortlisted for further investigations, i. Square with four regularly attached fins ii.Circle with four regularly attached fins iii.Hexagon with six regularly attached fins In order to meet the design constrain, three different geometric shapes were analysed using Multiphysics based CFD numerical solver (FENSAP-ICE [11]).Computational fluid dynamics based numerical modelling of atmospheric ice accretion on structures is a complex coupled process and mainly involve the airflow behaviour, droplet impingement & surface thermodynamics.First, it requires the airflow simulations, then super cooled water droplet behaviour is simulated to obtain the distribution of water impingement along the surface and finally the surface thermodynamic analysis are performed to estimate the ice growth.C-type structured numerical mesh was used, where mesh sensitivity analysis were carried out to accurately determine the boundary layer characteristics (shear stresses and heat fluxes).The one equation Spalart Allmaras turbulence model [12] was used as a compromise between acceptable computational cost and the required accuracy in simulating the turbulent flow.The initial sand grain roughness height (equivalent of surface roughness in CFD) for the surface was assumed 0.5 mm (0 mm means smooth wall) and also a variable surface roughness model of FENSAP-ICE based on the beads was used.Rotational effects were introduced by applying a tangential velocity on the rotating surface grid nodes.By describing the rate of rotation (RPM) in FENSAP-ICE, the components of tangential velocity were automatically computed for each surface node depending on its normal distance from the axis of rotation.Arbitrary Lagrangian-Eulerian formulation was used in FENSAP-ICE for the mesh displacement due to ice accretion with time.This approach adds the grid speed terms to the Navier-Stokes equations to account for the mesh velocity.At varying operating and geometric parameters it became possible to select optimum rotating geometric configuration that can be used as a rotating part of a proposed designed icing sensor.Different geometric and operating parameters were considered during this numerical study, To study the effects of geometric shape variations on ice accretion, three different cross sections (hexagon, circle and square) with fins, analysis were carried out at wind speed of 10 /,  = 35 and  = −10 0  for both rotating and non-rotating conditions.Figure 2 shows the velocity streamlines across each rotating geometric cross section.Results show a more streamline airflow behaviour for circular cross section, as compared to hexagon and square, particularly between fins.
Analysis were carried out to calculate the accreted ice growth and thickness for each cross section.Figure 3 shows the ice growth and thickness distribution along each cross section and results clearly highlights the uneven distribution of ice loads along the hexagon and square.Based upon these results it was found reasonable to reject square from the list of candidates.

𝑡 = 120 𝑚𝑖𝑛 10 0 𝐶
To study the effect of wind flow velocity variations on the rate and shape of ice accreting, analysis were carried out at two different wind velocities (10 / & 25 m/s) for rotating circular and hexagon cross sections.Figure 4 shows the accreted ice shapes and growth for both geometries.Analysis show that in case of hexagon the ice shapes get more irregular along the fins at high wind speed, which is mainly due to more irregular flow behaviour along fins.Based upon these results, it was found reasonable to reject hexagon from the list of candidates for CFD Simulations.

120𝑚𝑖𝑛 −10 0 𝐶
To understand the effects of fin size on resultant ice growth cfd analysis were carriedout for circular cross section.Two different sizes of fins 47 & 94  were used for this study.Fin size of 94mm was equal to one half of the longest diagonal of hexagon and 47mm was one fourth of the longest diagonal of hexagon.In order to compare the results, and keep adequate geometric similarity between different geometries, the fin size was kept constant for all.However to compare the effect of fin size, using optimal usage of resources only two sizes were compared.Figure 5 shows a considerable change in the ice growth with fin size.Decrease in fin size is considerably changing the flow behaviour along fins, which is resulting in development of flow recirculation zone along fin section.Analysis show a considerable increase in ice growth particularly along tip sections on fin in case of reducing its length.
To study the effect geometric rotation speed, the numerical analysis were carried out at an rpm of 6 and 30 (minimum rpm = 6 was selected because this was the minimum rpm at which motor was withdrawing stable current from the supply) of a circular cross section, at atmospheric air temperature of −10 0  and wind speed of 10 /.Figure 6 shows the ice growth for two different geometric rotational speeds.Increase in the  is increasing the tangential component of flow velocity, which leads to a change in the droplet behaviour and water run back, as Figure 6 also shows a development of ice hump near the fin root section.
This CFD based numerical study laid the foundation to understand and analyze the atmospheric ice accretion profiles on three different rotating geometric cross sections having fins.The results from this study indicated that the ice accretion on circular cross section (  = 162 ) rotating at an rpm of 6 with regularly attached four fins of dimension 94mm was relatively more uniform at different simulated geometric and operating parameters, than its two competitors, which were hexagon (ℎℎ = 162 ) with six fins and square ( ℎ = 162) with four fins.
The main purpose of this working model was to experimentally investigate and validate the numerical simulations finding.The other purpose of this working model was to experimentally validate rotational load measurement technique.Nevertheless, it was also developed to investigate the icebound clogging of the rotational unit with the stationary unit during the ice accumulation process.
The sensor should be spacious to install all the electronic cards, motor, and microcontrollers and should provide sufficient space for maintenance.Based upon the numerical simulation results and required space the CAD model of the sensor was prepared, see Figure 7.The dimensions of the rotational part were in agreement with the design found by numerical simulations.As required, this sensor would be divided into two units, rotational and stationary.The rotational unit would be used as measurement part and the stationary unit would be used as an electronics housing.
The rotating part of this sensor would have two parts (see Figure 7) which were rotating fins (capacitive cards) and rotational hub.The signals from the capacitive cards would be connected via a slip ring with a small microcontroller unit.This small microcontroller unit would be a part of central microcontroller which would also control the motor drive controller for maintaining a constant rpm.The motor would be located in the stationary unit.In order to select, it was important to calculate the static and dynamic loadings which the motor should withstand.
Rotational kinetic energy  was required to maintain a fixed rpm.This was calculated by Equation 1.The mass moment of inertia of the rotating system ′  ′ was found from the first CAD model of the sensor (see Figure 7 and Table 1).
In order to avoid dynamic effects of atmospheric ice, the additional work done on this motor could be represented by drag energy ..This was calculated using Equation (2).
Where ρ air is the density of air equal to 1.225kg/m 3 ,   is the velocity of air equal to 10m/s, A surface is the surface area of the exposed part which is equal to 2πRl for nonrotating geometries, with relative wind orthogonal to the surface, the surface area is considered to be half.However for rotating geometries, the total surface area need to be considered.Also as this calculation leading to find the maximum torque required by the motor therefore the overall geometry is considered.R is the maximum radius of the surface which is assumed to equal to 0.16m, l is the exposed length of the rotary part equal to 0.12m and c d is the coefficient of drag of the cylinder equal to 0.68.Therefore the total torque required to select the motor was Keeping in view the factor of safety under icing conditions the selected motor was 2 phase hybrid step motor ZSH 87/2, with a holding torque of 3.6Nm and detent torque of 0.05Nm.To electronically drive this step motor SMCI47-S drive controller was selected.The first working prototype of the sensor was manufactured using Zprinter 450 (a 3D Printing/Prototyping Machine at UiT) using ZP150 powder.From the CAD model .stlfile was developed which was then used in rapid prototyping machine.After the model was manufactured by Zprinter, the additional powdered ZP150 was required to remove from the part.This developed part was very fragile, hence after removing all the powdered material completely (using pressurized air blower, if required) the part had to be dipped in epoxy resin in order to strengthen the part.Once this part was strengthened using epoxy, it had be baked in the oven in order to harden the bond between epoxy and the hardened composite ZP150.

Part
Both techniques (rotational load measurement technique and capacitance measurement technique) were dealt separately and then coupled together in order to detect an icing event, determine ice type as well as measure icing load, icing rate and melting rate.This section is divided in two subsections.In Sec.3.1 it is aimed to measure icing load, icing rate and compare icing intensities using rotational load measurement technique.In Sec.3.2, it is intended to detect icing event, determine its type and measure melting rate using two different capacitive measurement techniques (time domain and frequency domain approach).
The experiments in this section are intended to measure atmospheric icing load and icing rate as a function of current demand.Furthermore, in order to compare the results with existing icing load sensor e.g.Ice Monitor icing intensities were calculated on both sensors.The experiments were conducted at Cryospheric Environmental Simulator, Shinjo, Japan.During these experiments it was aimed to prove use rotational load measurement technique to model the ice growth rate around a slowly forced rotational prototype sensor in order to calculate the icing load by measuring the variation in the current demand.The proposed experimental setup can be seen in Figure 8.The measurement setup and the accreted results can be seen in Figure 9.The prevailing conditions/parameters inside the tunnel are provided in Table 2 The calibration of the current sensor was done at a standard room temperature at varying rotating masses and rotational speed of the sensory system.In order to measure atmospheric ice accretion load and rate around a rotational geometry it was important first to mathematically represent the problem.The power input to a motor rotating at a constant rpm, can vary by changing the load on the rotating shaft of the motor.Based upon this Current-Inertia and Current-Mass relations were developed, see Equation (4) and Equation ( 5) respectively.For more detail on the derivation of these equations, see Mughal et.al. [13].
The focus of this first experiment (as discussed in Mughal et.al. [13]) was to compare the analytical, calibrated and experimental values of the parameters ,  and .During the experimentation, the atmospheric ice was only allowed to deposit on a rotating cylinder with four plates rotating at a constant rpm of 6.The base of the CES icing tunnel test section was used as a separation between rotating and stationary component so as to simulate the load and rate of accreted ice.These results were improved by introducing a forced rotation constant =0.7 (dependent upon the freezing fraction, rotation rate and thickness of the rotating plate).
From the results it was found that although there was a deviation between the experimental, calibrated and analytical values of the parameters  and  from due to natural constraints of lab based experimentation in the icing tunnel (size, scaling, electrical losses, etc.) but the results were linear.The graphical representation of the atmospheric icing load and icing rate can be seen in Figure 10 which are further tabulated in Table 3.
Furthermore, it was also found important to experimentally and analytically validate the CFD based numerical results of Section 2.2 and measure icing intensities on them.It was also found useful and compare the results with freely rotating IceMonitor.To meet this requirement it was found reasonable to do the following six step process, i.To estimate the collision efficiencies α 1 using Langmuir and Blodgett [14]  Following the above six step process α 1 were calculated for ISO standard (IceMoniter) using Langmuir and Blodgett [14] empirical relations.The results were then compared with the CFD based numerical model.The results were reasonably close with a deviation of 6%.Following the third step of six step process, ice accretion intensities were numerically calculated, the results are tabulated in Table 4. Based upon the slowly rotating icing intensity model of Makkonen in which he have proposed windward area ratio 'WAR' of 2/π for a freely rotating wire (for details see.Makkonen [15,16]) it was required in step four to determine wind area ratio of slowly forced rotating hexagon with six fins and slowly forced rotating cylinder with four fins.Through CFD based numerical simulation WAR for the slowly forced rotating hexagon with six fins was found to be 1/6 and for slowly forced rotation cylinder with four fins was found to be 1/2 (for details see Figure 2).As a rule of thumb this WAR varies from 0 → 1.The analytical ice accretion intensities on these geometries were then calculated using proposed models of Makkonen [15][16][17].By assuming freezing fraction α 3 = 1, relative  Finally an experimental setup similar with Figure 8 was established at CES, Shinjo, Japan to measure the intensities on all three geometries of interest.The graphical results of this experiment can be seen in Figure 11.The results are tabulated in Table 4.

Parameters
The capacitance of a material with polar characteristics could be measured by either using time domain or frequency domain approach [19].The constant   =   ′ −   ′′ is the complex dielectric constant which is dependent upon the excitation frequency, see [19].Water/Ice has strong dipolar effects particularly at low frequencies particularly as indicated in the experimental results of Evans [20], Stiles [21] and Kuroiwa[22].This dipolar orientation is generally associated with the relaxation phenomenon.Relaxation time  0 is a measure of the mobility of the dipoles that exist in the material to reorient themselves.In the frequency domain analysis, the relaxation frequency 'f c ' is indicative of the relaxation time, τ 0 = 1 2πf c ⁄ . Nevertheless this relaxation time constant τ 0 can also be detected without sweeping the excitation frequency using time domain approach by charging and discharging cycles and recording the time constant Using Multiphysics tools, some equations were derived by Mughal and Virk [23][24][25] and the results were compared with the experimental results of Evanes [20], Stiles [21] and Kuroiwa [22] which were found satisfactory.
This section is divided in two subsections.In the first subsection the time domai n approach is utilized to detect an icing event, determine ice type and measure melting rate.In the second sub section dc frequency sweep approach is utilized for the same.
Using this approach the time response of atmospheric ice was measured by continuous measuring the delay in the discharging time constant due to its presence on the electrode plate.The utilization of this technology for measuring atmospheric icing was new.Patented mutual charge transfer IC QT60240 [26] was considered suitable for the hardware implementation.For more details about the circuit design architecture see Mughal et.al [27].The final system was tested in the Cold Room Chamber of Campus Narvik, UiT.The system was inspected at room temperature and then temperature was reduced to -10oC.The set of experiments comprised of the combination of water and ice layer measurements formed/placed over the electrode panel.The set of measurements conducted were.Zero Crossover Timings [26] (zero crossover timing is not the same as the relaxation time.The greater the capacitance of measurand the lower will be the zero crossover timing) were recorded of different samples.The readings corresponding to different samples (hard rime ice generated in Cold Room, natural snow collected from outside campus and pure ice from freezer) placed over the electrode could be seen in Figure 12a and the result are tabulated in Table 5.By increasing the temperature from −10   → +20  C, the melting of ice layer/block occurred as a natural process.The τ ZeroCrossover measurements during the melting process were clearly observed, see Figure 12b.

>523 μs
In order to further investigate this technology, another series of experimentation were performed to measure salinity ratio SR, for details see [27].DC frequency domain approach was utilized using AD5933 IC [28].For more details about the circuit design architecture see Mughal et.al [27].The experimentations were done in the Cold Climate Chamber using following samples, i. Glaze ice frozen on the electrode plate at a temperature of -16C ii.The results of glaze ice at three different temperatures could be seen in Figure .There was no prominent peak observed in the dielectric constant variation (ε r ′′ = imaginary part) with the variation in frequency, hence it was not possible to compare the results of Stiles [29] and Kuroiwa [30].However the average values of the phase of the capacitance ((ε r ′′ /ε r ′ )) reflect observable trend.At the frequency of 6000Hz (relaxation frequency '  '), a change of phase is observed in the average value of capacitance phase.This value of   was found in agreement with the experimental results of Stiles [29].It was found that that the average value of this phase reached to a maximum value 20o by decreasing the temperature till−24  , hence more capacitance.The results were however were very noisy due to conductance of atmospheric ice ( Figure ).In another test, the results were compared with the natural occurring snow, outside the university campus and water samples.The results are shown in Figure 14.These results also indicated the same trend in the dielectric constant variation (ε r ′′ = imaginary part) with frequency due to very high conductivity values.However, it was found that the the average values of the phase of the capacitance ((ε r ′′ /ε r ′ )) reflect observable trend for all ice, snow and water, Figure 14.
The need of a slowly rotational atmospheric icing sensor to detect and icing event, determine icing type, measure icing load, icing rate and melting rate was identified in [2].To meet this need, two design constrains and two design objectives were emphasized.Based upon these design constrains, design objectives and keeping in view the ISO 12494 [3] recommendations, three different geometrical candidates were categorized.The selected geometric candidates were then analyzed numerically (using Multiphase CFD solvers) for reasonably uniform ice accumulation characteristics on their circumference, using coupled multiphase numerical analysis ( [10]).The optimized geometry selected from this numerical study was a circle with four regularly attached fins with a fin size of 94mm rotating at a constant rpm of 6 ( [10]).A working model was then manufactured using 3D Manufacturing Facility in order to validate the theoretical and experimental analysis.The working model was then tested to measure the accumulated icing load and icing rate based upon the rotational load measurement technique at the Cryospheric Environmental Simulator 'CES' in Japan ( [13]).Experimental testing of prototype icing sensor for rotational loading proves the associated physics for an adequate measurement of icing load and icing rate.It was observed that due to the rotation mechanism the distribution of ice was uniform and no obvious ground effects were visible.Also there were no indications of ice deposition between the sensor and the test section base and even on the sensor top surface.These ground effects were prominent on ice load monitor and non rotating polygon.Hence, if the rotating fins are assumed to be some electrode plate then reasonably uniform deposition of atmospheric ice on the plates may support the electronic measurements to some extent.
Using the same CES experimental facility, the accumulated atmospheric icing intensity were experimentally determined on two geometries, hexagonal prism with 6 fins and cylinder with 4 fins.The results were compared with analytical models of Makkonen ( [15,17]) and also compared with the freely rotating ISO Standarized IceMoniter ([6]).It was deduced that the analytical relations of Makkonen can still be considered as a generalized case for forced rotating geometries.However the collision efficiency  1 can be assumed an explicit function of windward area ratio  and forced rotation constant  ( = 1 for freely rotating geometries).It was also found that the experimental icing intensity on cylinder with four fins slowly rotating at rpm of 6, was reasonably close to freely rotating IceMonitor.Based upon the experimental observations, the following can be deduced that the constantly slowly rotating cylinder with four fins is still a better geometry to measure incremental mass moment of inertia due to the reasonably uniform ice accretion as predicted by CFD Simulations.It was also identified in this study that the six step approach would be very useful to carefully understand the ice accretion on different rotating and non rotating geometries.
Both time domain and frequency domain approaches were utilized to measure the capacitance variation in atmospheric ice during ice accumulation process in order to detect the icing event, determine its type and measure the melting rate.The experiments were performed at the cold climate chamber facility at the Institute of Industrial Technology, UiT ( [27]).From the viewpoint of measurement, time domain measurement technique is simpler than the frequency domain measurement technique.Preliminary results reflect that it was possible to detect an icing event using frequency domain measurement technique but it was not possible to measure the melting rate using frequency domain capacitance measurement technique.However time domain capacitance measurement technique reflected that it is more useful and direct method than frequency domain capacitance measurement technique because using this approach it was possible to detect an icing event, determing ice type and its melting rate during the ice accumulation process.
The complete sensor was located in the Cold Room Chamber.Before starting the experiment, the sensor was tested and it was found that there were electromagnetic radiation problem due to the interference of step motor circuitry and the pcb circuitry which leads to deliver some erroneous data.This problem was avoided by wrapping the motor in an aluminum foil and grounding it and wrapping the pcb in a transparent plastic foil.The cold room was then operated at a temperature of -10 o C and ice was accumulated on the electrode plate.Atmospheric was detected, type was identified and melting rate was measured using time domain capacitance measurement, whereas icing load and rate were measured using rotational load measurement.
The initial design study to develop an atmospheric icing sensor have been laid out which have the potential capability to detect an icing event, determine icing type, measure icing load, icing rate and melting rate.Based upon this initial study a robust prototype of the sensor could be developed which can be coupled with a new generation of Cold Region Weather Station and can also be coupled with an Intelligent Ice Protection System for onshore applications and offshore applications.
velocity v = 10m/s, α 2 = 1, liquid water content W = 0.5g/m 3 , collision efficiency α 1 from numerical model, we estimated the analytical values of ice accretion intensities using Equation (6) on the geometries of interest. =  1  2  3  (6) i. Water layer measurements over the electrode panel a.Small water puddles over the electrode surface b. large water puddles over the electrode surface ii.Ice placed over the electrode surface a. Pure ice from freezer b.Hard rime ice from Cold Room Chamber c. Natural Snow from ground outside university campus d.Transitions from ice to wet ice e. Pure ice to water transformation rate Figure iii.Glaze ice frozen on the electrode plate at a temperature of -20C iv. Figure v. Glaze ice frozen on the electrode plate at a temperature of -24C vi. Figure vii.Natural Snow from ground outside university campus, Figure 14.viii.Normal tap water, Figure 14, ix.Pure ice to water transformation rate.
[15,16]al relations on ISO Standarized icing sensor.ii.To use numerical simulations to determine ice accretion intensities on a ISO Standarized icing sensor and compare the results with the Langmuir and Blodgett empirical model[14].iii.If the results are reasonable (within range) then to do numerical simulations on the rotating geometries to estimate collision efficiency and ice accretion intensities numerically.iv.To numerically estimate the Windward Area to use the Makkonen Model from[15,16].v.To analytically estimate ice accretion intensities by using the reasonable estimates of W and α 1 .vi.To compare atmospheric icing intensities through numerical, analytical & experimental results.
The work reported in this paper was partially funded by the Research Council of Norway, project no.195153 and partially by the consortium of the project ColdTech-Sustainable Cold Climate Technology.