Abstract
We propose to model ice shedding trajectories by an innovative paradigm that is based on cartesian grids, penalization and level sets. The use of cartesian grids bypasses the meshing issue, and penalization is an efficient alternative to explicitly impose boundary conditions so that the body-fitted meshes can be avoided, making multifluid/multiphysics flows easy to set up and simulate. Level sets describe the geometry in a nonparametric way so that geometrical and topological changes due to physics and in particular shed ice pieces are straight forward to follow. The model results are verified against the case of a free falling sphere. The capabilities of the proposed model are demonstrated on ice trajectories calculations for flow around iced cylinder and airfoil.
1. Introduction
Ice accretion on aerodynamic bodies is a serious and not yet totally mastered meteorological hazard due to supercooled water droplets (liquid water droplets at a temperature below the dew point) that impact on surfaces. Ice accretion is a multiphysics phenomenon [1] including fluid mechanics, heat transfer, and solid mechanics. Ice accretions have several negative effects, especially performance degradations for both aircraft and wind turbine. For aircraft, sudden performance degradations due to ice accretion cause several incidents and accidents each year [2]. Performance degradations include substantial reduction of engine performance and stability, reduction in aircraft maximum lift and stall angle, and increase of drag. For wind turbine, ice accretion on blades is a major concern in northern climate. When adverse meteorological conditions occur, ice accretes at the outer part of the blade with an approximately linear increase with time [3]. With growing ice accretion, the blade drag increases, diminishing the power output of the turbine and eventually causing a complete loss of production [4]. Other negative consequences include overloading due to delayed stall [5], increasing fatigue of components due to imbalance in the ice load, and damaging or harm caused by uncontrolled shedding of large ice chunks [6, 7].
In practice, ice accretion can be minimized by deicing systems [8] or prevented by anti-icing systems. By reducing the adhesive shear strength between ice and surface, de-icing systems remove ice formed on the protected surfaces following a periodic cycle. This cycle is defined such that intercycle ice shapes remain acceptable from a performance point of view. Because they are activated periodically, de-icing systems require less energy than anti-icing systems. Still, because of the enormous amounts of thermal or mechanical energy that would be required, complete removal of ice formation on large structure is not economically feasible. Only areas sensible to ice accretion, such as wing or blade leading edge, benefit from ice protection systems.
Actual concerns about greenhouse gases will lead to changes in the design of ice protection systems, thus renewing interest in de-icing simulation tools. To save fuel burn, aircraft manufacturers are investigating “new” ice protection systems such as electrothermal or electromechanical de-icing systems to replace anti-icing systems. One of the drawbacks of de-icing device is the ice pieces shed into the flow. The knowledge of ice shedding trajectories could allow assessing the risk of impact/ingestion on/in aircraft components located downstream [9]. When they leave the aircraft surface, ice pieces become projectiles that can hit and cause severe damage to aircraft surface or other components, as aircraft horizontal and vertical tails, or aircraft engine [10]. Aircraft certification authorities, such as FAA, have specific requirements for large ice fragment ingestion during engine certification. Control surfaces or wing flaps are also sensitive to ice shedding because they can be blocked by ice fragments. Aircraft manufacturers rely mainly on flight tests to evaluate the potential negative effects of ice shedding because of the lack of appropriate numerical tools [11]. The random shape and size taken by ice shed particles together with their rotation as they move make it difficult for classical CFD tools to predict trajectories.
For the wind turbine industry, new environmental preoccupations mean building more wind turbines in less favorable sites, especially in northern countries like Canada. Wind turbines are going to be built in sites prone to atmospheric icing, and the use of de-icing systems may reduce production losses. Even if no ice protection systems are installed, vibration and aerodynamic forces can cause ice shedding. Ice shed trajectory computations are then needed to reduce damage risk.
Ice shedding phenomena is an emerging topic for researchers working in numerical simulation of ice accretion. Until recently, only a few research works have paid attention to it. In one of the first numerical works published on de-icing system, De Witt et al. [12] studied the three-dimensional transient heat transfer in a multilayered body that is ice covered. This work focuses mainly on modeling the phase change in the ice layer with the movement of the solid-liquid interface as the latent heat of fusion is absorbed or released. In an article by Henry [13], criteria based on water film thickness is used to predict ice break-up. Once the piece of ice break up, it simply leaves the computational domain, and no calculation of its trajectory is done. Stresses in accreted ice caused by aerodynamic forces have been studied by Scavuzzo et al. [14], and an ice failure criterion [15] based on both normal and shear interface stresses has been used in a finite elements code to simulate the electroimpulse deicing process of aircraft wings.
Kohlman and Winn [16] proposed a trajectory simulation method to compute the trajectories of ice particles, represented by square plates of uniform thickness, into a uniform velocity field. The lift and drag were obtained from empirical correlations, and initial displacement and rotation were imposed to initiate trajectory calculations. Santos et al. [17] used a similar method but, this time, the trajectories were calculated into a nonuniform flow field around a wing. The ice fragment, a square flat plate, was released from a point in front of the leading edge of a wing. Initial position and velocities were varied, and the probabilities of ice impact at a location two chords downstream of the leading edge were obtained. Another method, based on a modified water droplet trajectory code, was used by Chandrasekharan and Hinson [18] to track trajectories of an ice disk and two ice slabs. The drag of the disk was assumed to be that of a sphere and the drag of the slabs to be that of a flat plate normal to the flow. Recently, Papadakis et al. [9] presented a statistical approach to perform trajectory computations for ice fragments shed from the wing and fuselage surfaces of a business jet. They carry an experimental study of aerodynamic loads around a potential ice fragment and derived empirical correlations. Those correlations have been used into a methodology based on trajectory calculation and probabilistic approach to identify areas where ice fragments most likely strike the aircraft.
On one side, some numerical codes predict ice break-up, and on the other side different numerical codes compute ice shedding trajectory. The numerical simulation of a full unsteady viscous flow, with a set of moving bodies immersed within, shows several difficulties for grid-based methods. Drawbacks income from the meshing procedure for complex geometries and the regriding procedure in tracing the body motion. A new approach capable of solving the complete de-icing problem is proposed in this paper. The approach is based on cartesian grids, penalization, and level sets.
The objective of the present paper is to demonstrate the approach capability to predict ice shedding trajectory, starting just after the ice break-up. First, a vortex method [19] is proposed to simulate the interaction of an incompressible flow with rigid bodies. For our method, we followed the innovative approach proposed by Coquerelle et al. [20, 21]. They propose an efficient and accurate technique to simulate unsteady incompressible viscous flows using hybrid vortex methods [22]. Hybrid vortex methods are based on the combination of Lagrangian mesh-free schemes and Eulerian grid-based schemes [22]. The Navier-Stokes equations are formulated in terms of the vorticity formulation, and the vorticity field is numerically determined by a particle discretization. The two schemes solve different terms of the equation, for example, the Lagrangian scheme is used to solve the nonlinear advective part of the equation while the Eulerian scheme is used to solve on the grid the diffusive part of the equation and the velocity term. A penalization method [23] is used to enforce the no-slip boundary condition inside the solid wall boundaries. The bodies around which the flow is computed are defined using the so-called penalization method or Brinckman-Navier-Stokes equations in which the bodies are considered as porous media with a very small intrinsic permeability. Level set functions are used to capture interfaces and compute rigid motions of the solid bodies [24].
After presenting the model equations, the numerical procedure is exposed. Then, numerical trajectories results are presented for sedimentation of a 2D cylinder on a flat plate, flow around an iced cylinder and an iced airfoil.
2. Fluid-Solid Flow Model
2.1. Physical Model
Given a computational domain , we consider incompressible flow in around rigid solids . A schematic representation of a computational domain composed of two solids is sketched in Figure 1.
We assume that the density is constant in the fluid and in the solids. The fluid-solid interaction problem can be modeled by the incompressible Navier-Stokes equations In the above system is the velocity, is the kinematic viscosity, is the density, is the gravity vector, and is the pressure. Depending on the location in the computational domain, the velocity is either the fluid velocity or the solid velocity . The idea is to extend the velocity field inside the solid body and to solve the flow equations with a penalization term to enforce rigid motion inside the solid [21].
Given a penalization parameter and denoting by the characteristic function of the solid , the model equation is the following: coupled with the incompressible mass conservation equation
The velocity of the moving rigid motion is obtained by averaging translation and angular velocities over the solid, with the following equation: where is center of gravity, , its inertia matrix, and is the coordinate at calculation point inside the domain. The rigid body follows the trajectories of the flow with the advection velocity . Therefore, characteristic function can be obtained by solving the advection equation
In this paper is computed from a level set function satisfying the same advection equation
We initialize as the signed distance function to the boundary of , because is negative outside and positive inside the solid: where is the Heaviside function.
It is important to notice here that, because is a rigid body motion, one can guarantee that remains a signed distance for all time. To summarize, each body-fluid interface is captured by a level set function. These level set functions are moved in rigid motion by advection.
2.2. Ice Shedding Law
Ice accreted remains attached to surface until the adhesive bond between ice and substrate fails or until a fracture occurs in the ice. The ice adhesive shear strength may depend on substance, surface roughnesses, and coating on the surface side. The ice adhesive shear strength depends also on atmospheric conditions present during the accretion process: temperature, liquid water content, air velocity, and so forth [25]. Generally, the adhesive shear strength for rime ice is lower than for glaze ice.
Usually, the aerodynamic forces alone are not sufficient to detach accreted ice from airfoil. To deice airfoil, either the adhesive shear strength of the ice is reduced almost to zero by heating and melting the accreted ice, or mechanical forces break the ice and give an impulsion sufficient to detach ice from surface. In the two cases, shear stress inside the ice and, to a lesser extent, normal stress play a critical role in shedding.
The choice of an ice shedding law dictates the initial conditions for the ice trajectory calculation. At the end of the break-up process, an ice piece will have a definite shape, a translational velocity, and a rotational velocity that becomes input for the trajectory calculations, see [9, 14]. Although some researchers have proposed models [26], these models are not commonly accepted and used to predict ice break up by the aircraft icing community.
In this paper we concentrate on ice trajectory calculation, and a method similar to the one proposed in [16, 17] is used to initiate trajectory calculations. In these papers, a statistical approach is used, and ice fragments are launched from several location points near an aircraft body. From the several trajectories, the probability of having an ice fragment in a given area near the aircraft is devised. In the present paper, an initial displacement, together with an initial velocity and initial rotation velocity equal to zero, is imposed on the ice piece to make it leave the body surface.
2.3. Numerical Method
In this paper, the penalized Navier-Stokes equation is defined in terms of a vorticity formulation, and the vorticity field is numerically determined by a particle discretization [19]. Let us consider the penalized Navier-Stokes equation in the vorticity formulation by applying the curl operator to (2) with The density is computed from the fluid and solid densities, respectively, and , using the following equation: Equation (8) is required to compute the variation of the pressure. Usually, in vorticity formulation, this is not necessary. Indeed if the density is constant the term vanishes. To avoid an explicit computation of the pressure, (8) can be reformulated using (1), in the following way:
Following Cottet and Koumoutsakos [22], the hybrid vortex methods are based on the combination of Lagrangian mesh-free schemes and Eulerian grid-based schemes on the same flow region. The Vortex-In-Cell (VIC) scheme is an example of hybrid vortex method: the non linear advection is computed by tracking the trajectories of the Lagrangian particles through a set of ODEs, whereas an Eulerian grid is adopted to solve efficiently the velocity field, the diffusive term, and the penalization term. Given , the material derivative, by expanding the penalization term, (11) becomes where is the 1D Dirac delta function and . We recall that with level-set functions, is the unit normal to the interface oriented towards the solid. The last two terms on the right of (12) play a significant physical role in the model. The first term clears the vorticity difference within the bodies, whereas the second member represents a vorticity generation term that is localized on the solid boundaries and allows the no-slip condition to be imposed. The penalized vorticity equation means that the rate of change of the vorticity, advected by the fluid in a Lagrangian frame of reference, is governed by the diffusive effects, the stretching effects, the production of bound vorticity, and the vorticity cancellation within the solid bodies.
Let us introduce the level-set Vortex-In-Cell (VIC) algorithm. The domain is meshed using a uniform fixed cartesian grid. We denote by the time step such as and , , are grid values of the level set functions, velocity and vorticity, respectively. In vortex methods, the rate of change of vorticity is modeled by means of discrete vortex particles, such that the solution of (12) is localized only in the rotational regions of the flow field. This is the most important advantage of the vortex methods, that is, the computational efforts are naturally addressed only to specific flow field zones. The vorticity field is represented by a set of particles where is the number of particles, the particle location, and , are, respectively, the volume (constant due to the incompressibility) and the strength of a general particle . is a smooth distribution function such that which acts on the vortex support. The interpolating-remeshing scheme is a fundamental tool for the accuracy of the whole method.
We propose a splitting algorithm to solve (8). Each time step is solved using four substeps.
Substep 1 (advection). During the first substep, the vortex elements are advected by the local flow velocity Grid vorticity above a certain cut-off is used to create particle at grid point locations [27], as shown in Figure 2(a). Using (14), particles are displaced with a fourth-order Runge-Kutta time-stepping scheme (Figure 2(b)). From the new vortex particles's locations, the vorticity field is remeshed on the grid by means of an interpolation procedure (Figure 2(c)).
(a)
(b)
(c)
The interpolation order is directly linked to the number of moments preserved by the new particles' distribution by comparison with the former one. In this paper, we used the following third-order interpolation kernel introduced by Monaghan [28]: where is the distance from the point to interpolate. In this 1D example, the influenced stencil consists of four surrounding grid points. In the 2D case, the scheme takes into account the sixteen closest grid points around the particle to interpolate/remesh. In Figure 2(c) the particle interpolating/remeshing scheme is illustrated.
Substep 2 (diffusion). Once Substep 1 is done, the remaining parts of the governing equation (11) can be solved. The equation to solve for the viscous contribution is then Equation (16) is approximated onto the grid by means of an Euler explicit scheme, while the Laplacian is evaluated, with a second-order accurate standard centered finite differences five points stencil.
Substep 3 (pressure gradient contribution). The equation to solve for the pressure gradient contribution is
To solve (17), density values are obtained from (10). Grid values of and are used to compute and centered finite differences at are used to compute .
Substep 4 (penalization). Finally, the penalization term is evaluated using The discretization and the integration of the penalization term affect the choice of the penalization parameter (the larger the parameter , the better the quality of the penalization in practice in our simulations is fixed to ). An Euler explicit time discretization of (18) does not allow to use . An implicit Euler time discretization is therefore used for the penalization term in the Navier-Stokes equation: The vorticity field at is evaluated on the grid by taking the curl of the velocity and computing the derivative though the second-order centered finite differences approximation, while is evaluated using (4). This method is unconditionally stable.
Because the equations are not written in primitive variables, special treatments are needed to recover the velocity field and to impose the boundary conditions. Since the incompressible velocity field is divergence-free, from the vector field theory, we can define a vector potential such that The vector potential is a 3D extension of the so-called stream function . This potential vector is imposed to be solenoidal and given the updated vorticity field, the stream function field is computed by solving the linear Poisson equation, on the cartesian grid with boundary conditions on .
In (19), is computed using (21) and (22) with the vorticity resulting from Substep 3.
In vortex methods, boundary conditions are explicitly used only to solve the Poisson equation and are enforced on the nonprimitive variable . In the present method, on the upstream and downstream boundaries (right and left boundaries on Figure 1), a Neumann condition, , is enforced. On the lower and upper boundaries, a Dirichlet condition sets the flow mass rate through the domain.
In our simulations, a Fast Poisson Solver is adopted (Fishpack90 library [29]).
For the present method, following [30], the forces can be evaluated without any information on the pressure field but only through the knowledge of the velocity and vorticity fields.
3. Numerical Simulations
The proposed approach results are validated against available literature results for the sedimentation of a 2D cylinder on a flat plate. Then, ice trajectory prediction capabilities are verified for two test cases at low Reynolds number. The first one, an iced cylinder, is representative of an unsteady flow around a bluff body. The second one, an iced airfoil, is more representative of the flow around an aerodynamic body, although with the attached ice shape, it behaves initially like a bluff body. For sake of aerodynamic force computations verification, the fluid and solid densities are assumed to be the same in the ice trajectory calculations. Inertia and gravity effects are thus negligible.
3.1. Sedimentation of a 2D Cylinder on a Flat Plane
We consider the case of a 2D cylinder in a square cavity, falling under gravity on a flat plane. This test is used to verify the method by a comparison with Glowinski [31] and Coquerelle et al. [20]. The dimension of the cavity is . The viscosity is 0.01. The density inside and outside the cylinder is, respectively, 1.5 and 1. The cylinder has a radius of 0.125, no roughness, and is initially located at the point . It accelerates under gravity, set to , then settles to a steady velocity, due to the friction forces, and eventually hit the bottom of the cavity and stops.
In Figure 3, the evolution of velocity component along a line, , cutting the middle of the cylinder is drawn at time 0.3. The cylinder has almost reached its final velocity. Between and the velocity is the one inside the solid cylinder. Fluid is accelerated around the cylinder due to mass conservation principle. On the side of the cavity, fluid sticks to the wall, and thus the velocity is 0. To impose this boundary condition on the cavity walls a penalization technique has also been used.
In Figure 4, the barycenter computed position and the evolution of its vertical velocity as a function of time are plotted. The solution is in agreement with Coquerelle et al. [20] and Glowinski [31] solutions. Contrary to literature results, our cylinder simply stops when it hits the wall instead of bouncing. This is to be expected because no collision model is used in our simulation. Nevertheless, result comparisons indicate that inertia, gravity, and drag forces are correctly implemented in our code.
(a) Evolution of the y coordinate of the cylinder barycenter
(b) Evolution of the component V of the velocity
3.2. Iced Cylinder
We consider first an incompressible flow around two static rigid solids, the cylinder and the ice piece , in . A schematic representation of the computational domain composed of these two solids is sketched in Figure 1. The domain is a rectangular box . The diameter of the circular cylinder is one unit, the free stream velocity is one velocity unit, and the Reynolds number is defined by . The whole computational domain is meshed by an equispaced Cartesian orthogonal grid . The choice of this mesh spacing is justified by a grid sensitivity study presented further in this section.
The ice piece, solid , is defined by the area between the curve of (23) and the line of equation followed by a rotation of 45 degrees around cylinder center:
For the subsequent simulations, the flow field is computed by solving the Poisson problem with a homogeneous Neumann condition, , at downstream and upstream boundaries, while a Dirichlet condition of the stream function, , is imposed on top and bottom boundaries. In particular the potential flow around the cylinder is considered, and the value of the associated stream function is enforced, that is, This is equivalent to imposing a symmetry condition on top and bottom.
The undisturbed stream passing through the cylinder is used as initialization. The Reynolds number of the flow is . Following [32], the time step for calculations is evaluated with the aim of solving with accuracy the diffusive phenomena. Since for detecting the principal diffusive scales with efficiency and accuracy the Von Neumann number (VNN) has to be in the order of unit, the nondimensional time step is given by In present computations , then .
Figure 5 shows the distance level set function used for the penalization. The black line corresponds to the isoline 0, that is, it represents the fluid solid interfaces.
A grid sensitivity analysis using Mesh 1: , Mesh 2: , Mesh 3: , and Mesh 4: has been performed to select the appropriate grid. Figure 6 shows the evolution in time of the lift coefficient, , around both the cylinder and the ice chunk for the four different meshes. Numerical solutions obtained with Mesh 3 and Mesh 4 are in very good agreement for the lift coefficient amplitude. A slight frequency difference is still visible, causing a shift in the maxima and the minima position between Mesh 3 and Mesh 4, visible for time greater than 80. For the following study, Mesh 3 has been selected to perform ice shedding trajectory calculations.
Using Mesh 3, the evolution of the lift coefficient, , and drag coefficient, , around cylinder and ice piece is presented in Figure 7(a). The flow is well established around the two static solids at about nondimensional time . On the drag coefficient curve, there is two maximum values, indicating two distinct vortex shedding time. This is different from the usual drag coefficient curve obtained for a single cylinder in a flow.
(a) Evolution of lift and drag on both cylinder and ice piece
(b) Distribution on the drag period of the 17 shedding start times (every 24 time steps)
Once the flow is established, the ice piece is released in the flow. Because of the vortex shedding induced by both the cylinder and the ice piece, several times of shedding start-up are investigated. More precisely seventeen starting times are equally distributed (every 24 time steps) along an established period (Figure 7(b)) to observe the effects of the vortex shedding on the ice trajectory. For the present paper, the 17 numerical simulations of ice shedding will be termed as shed 1, shed , and shed 17, each one corresponding to the starting time 1, , and 17 of Figure 7(b).
Once ice shedding starts, the ice piece freely evolves in the incompressible flow contained in , whereas the cylinder remains static, modeling a simplified situation of ice shedding. For the static body , is fixed at any time. The ice piece velocity, , and angle of rotation are computed according to (4) and injected in the penalization term. Recall that the ice piece is assumed to be a rigid solid body without any deformation. Figure 8 shows the vorticity contours of the flow corresponding to shed 13. In this figure, green contours correspond to positive vorticity whereas blue contours correspond to negative vorticity. Notice that once ice shedding starts, the time step does not remain constant and depends on the ice moving velocity .
(a) Initial solution
(b) Solution after 2800 time steps
(c) Solution after 5600 time steps
(d) Solution after 8400 time steps
As it is usually the case for flow around bluff bodies, the ice body creates its own vortex shedding. The ice initially moves downstream, for the first 2800 time steps. Then it gets trapped by the flow in one of the vortices behind the cylinder and start moving downward until around 5600 time steps, then upward around 8400 time steps, almost stopping its downstream motion.
In Figure 9, we investigate the effects of the vortex shedding on the ice trajectories. For the 17 simulations, the ice piece barycenter trajectories are drawn in Figure 9(a). As expected, the ice trajectory is quite different from one solution to another, especially behind the cylinder, around . This random aspect of the trajectories calculations is also observed by Papadakis et al. [9], and they use a statistical approach to take into account the randomness of the shedding process around an aircraft. Farther downstream, after , trajectories tend to get closer. In all the cases, no ice trajectory goes below . Details on the top trajectory (shed 3) in red, the bottom trajectory (shed 13) in blue and a middle trajectory (shed 10) in green are drawn in Figure 9(b).
(a) Ice piece barycenter trajectory for the 17 sheddings
(b) Typical ice shedding trajectories: in red shed 3; in green shed 10; in blue shed 13
(c) Evolution of the angle of rotation for the 17 sheddings
Vortex shedding has also a significant influence on the ice angle of rotation. Due to the low ice density used in this verification, inertia forces are negligible, and ice piece align quickly with the flow. The ice piece rotates and tends to keep its round side upstream, in a low drag orientation. The present method enables us to plot the rotation together with the displacement on Figure 9(c).
3.3. Iced NACA 0012
The computational domain shown in Figure 1 can be adopted for the ice shedding test over an infinite wing, by replacing the circular cylinder with an iced NACA 0012 airfoil. For this numerical test, the ice piece geometry (see Figures 10(a) and 10(b)) is numerized from Papadakis et al. [9] and scaled to be placed on a NACA 0012 airfoil of chord unity. The ice piece geometry with its horn shape is representative of a large glaze ice accretion. This kind of ice shape is more susceptible to break up due to aerodynamic forces. Such a large ice shape tends to have high pressure forces on the upstream side and low pressure forces on the downstream side, creating a strong bending moment at airfoil surface. In Figure 10(c), the initial displacement of the ice shape, just after breakup, is shown. This kind of initial displacement may be generated by mechanical devices, such as inflatable boot.
(a) Ice piece on a NACA 0012 airfoil
(b) Zoom on the ice piece geometry
(c) Shedding start-up (translation imposed for the first time step of shedding)
In such case the boundary conditions for the Poisson problem are modified as follows.(i)Dirichlet condition in order to enforce the mass flow at upstream (left boundary). (ii)Neumann condition at downstream (right boundary). (iii)On the bottom and top boundaries Dirichlet boundaries conditions are prescribed ( and , resp.) in order to enforce a far-field undisturbed flow condition.
First the flow is computed around the static iced NACA 0012 airfoil. The domain is a rectangular box of size and the size of the cartesian grid used is . The Reynolds number of the flow, based on airfoil chord , , is 1000, and the dimensionless time step is . For this case, the distance level set function used for the penalization is plotted in Figure 11(a) for initial flow and in Figure 11(b) for the start of trajectory calculation.
(a) Before shedding
(b) After impulsion
The periodic flow is well established at about , see Figure 12. Because of the ice shape, the drag is higher than on a clean airfoil, and a positive lift coefficient is generated, even if the airfoil's angle of attack is 0.
Figures 13 and 14 show, respectively, the vorticity and the velocity for this established solution. In blue, the vorticity is negative, and in green, the vorticity is positive. The ice shape creates massive flow separation on the upper side of the airfoil, and large vortices of opposite vorticity sign are shed. As expected in the case of massive separation, there is a wake of vortices behind the airfoil. The velocity is accelerated outside of the separation area on the upper side of the airfoil. The separation area is indicated by the negative velocity component area behind the ice shape. This negative area extend all over the airfoil upper side.
Then ice shedding is initiated. To model ice shedding and be sure that the ice piece will not collide with airfoil, we simply give an impulsion to the ice chunk, as if a deicing boot device has been used. The presented method can handle collisions between bodies, but it would involve defining an additional physical model. For sake of trajectory verification, the problem is kept as simple as possible. The impulsion is given by the vector , see Figure 10(c). The distance level set function, after shedding, is presented in Figure 11(b), and the level set changes with the ice piece motion.
Details of the first moments after ice shedding are shown in Figure 15. On the left, vorticity contours are plotted and, on the right, velocity components. The ice shape starts to move in the normal direction from the airfoil, then, it starts to rotate to align itself with the flow field around airfoil. Obviously, the vorticity field and the velocity field change drastically between initial time (Figures 15(a) and 15(b)) and final time (Figures 15(e) and 15(f)). The vorticity pattern, for example, in Figure 15(c), shows that the ice piece is slower than airflow because a secondary wake is created around the shape.
(a) Solution after 100 time steps
(b) Solution after 100 time steps
(c) Solution after 500 time steps
(d) Solution after 500 time steps
(e) Solution after 1000 time steps
(f) Solution after 1000 time steps
Finally, a sample trajectory is studied in more detail in Figure 16. Because the flow field around the clean airfoil is simpler than around a cylinder, the ice shape simply flows downstream in Figure 16(a). The details of the initial movement of the shape gravity center are plotted as a red line. The shape slowly rotates and moves downstream, as Figure 16(c) shows. After a quick rotation between time 50 and 60, because of the low inertia of the ice piece, the angle of rotation becomes almost constant. The ice shape is outside of the boundary layer, even if this one is thick due to the low Reynolds number of the flow.
(a) Ice chunk trajectory
(b) Zoom on leading edge
(c) Evolution of the angle of rotation
4. Conclusions
The proposed approach based on cartesian grids, penalization, and level sets is promising. The originality of the numerical tools used enables to take easily into account the topology changes. The fluid-solid flow model has been described in detail together with the numerical method used to solve the problem. The trajectory calculation capabilities have been verified against other numerical results for the sedimentation of a cylinder on a flat plate. The approach enable us to calculate ice shedding trajectory around aerodynamic shape and bluff body. Verification of the results for bluff body calculation shows that trajectories depends strongly on flow condition at the shedding time for unsteady problem. Ice piece trajectories distribute themselves almost randomly around an area. For the rotation angle however, the ice chunks tend to orient themselves with the flow such as to reduce the drag. Further verifications of the approach for trajectory calculations around an iced airfoil show first that large separation occurs on the upper side of the airfoil, making the flow unsteady at the low Reynolds number studied. Secondly, ice trajectories are much simpler due to the fact that once ice separates from the airfoil, large vortex disappears from the flow.
The next step is to add high Reynolds number capabilities to the flow solver, more representative of wind turbine or aircraft physics. Model developments are also needed to study ice break-up and predict the correct initial conditions for trajectory calculations.