# Parametric Studies of Flat Plate Trajectories Using VIC and Penalization

Authors: François Morency; Héloïse Beaugendre
Journal: Modelling and Simulation in Engineering
 (2018)

---

## Abstract

Flying debris is generated in several situations: when a roof is exposed to a storm, when ice accretes on rotating wind turbines, or during inflight aircraft deicing. Four dimensionless parameters play a role in the motion of flying debris. The goal of the present paper is to investigate the relative importance of four dimensionless parameters: the Reynolds number, the Froude number, the Tachikawa number, and the mass moment of inertia parameters. Flying debris trajectories are computed with a fluid-solid interaction model formulated for an incompressible 2D laminar flow. The rigid moving solid effects are modelled in the Navier-Stokes equations using penalization. A VIC scheme is used to solve the flow equations. The aerodynamic forces and moments are used to compute the acceleration and the velocity of the solid. A database of 64 trajectories is built using a two-level full factorial design for the four factors. The dispersion of the plate position at a given horizontal position decreases with the Froude number. Moreover, the Tachikawa number has a significant effect on the median plate position.


---

## Body

## 1. Introduction

Flying debris, such as ice fragment, can cause serious damage to structures. One way to mitigate the risks associated with their impacts is to use numerical tools to predict their trajectories. However, ice fragment trajectories are stochastic, as many parameters affect the flying path and the flow [1]. In particular, trajectories are sensitive to the size, shape, location, and initial velocity of ice fragments [2]. Consequently, probability distribution maps are built to show the likely path of such fragments, which can lead to the development of risk mitigation solutions. A large number of trajectories are needed to obtain statistically relevant maps in 3D. Mathematical models enable the simulation of multiple trajectories for various initial conditions, which then allows the computation of probability maps around a geometry [3, 4]. The motion of the ice piece is driven mostly by aerodynamic forces. Research works covering wind engineering have established that three dimensionless parameters play a role in the motion of flying debris, apart from the Reynolds number that may play a role in aerodynamic forces [5–7].For low fidelity models, motion equations are solved with a Runge-Kutta method in a Lagrangian frame of reference attached to the debris. The aerodynamic forces are computed using the static lift, drag, and moment coefficient correlations for the studied debris shape, neglecting the shape effects on the flow field. The two most studied shapes are the sphere [8] and the rectangular flat plate [6]. For the flat plate, the coefficient correlations depend only on the angle of attack, and additional aerodynamic terms must be added to take into account debris rotation [9]. For the sphere, the coefficient correlations depend only on the Reynolds number.High fidelity models involve a tight coupling of the flow field to the flying debris motion. These models use a computational fluid dynamic (CFD) to solve the Navier-Stokes equations around the moving debris in an Eulerian frame of reference [7, 10], and the forces and moments are computed from the flow solution. The computational costs of such models are higher than those for models uncoupled from the flow field, especially as the flow Reynolds number increases. However, high fidelity models allow the computation of the trajectory of any kind of shape and do not require experimental databases or correlations for aerodynamic forces and moments. Also, high fidelity models take into account Reynolds number effects and lift and drag fluctuations (10% to 20%) caused by vortex shedding [11, 12] in the case of flow separation around debris. Both models are complementary and need research developments to better assess the risk of impacts.The goal of the present paper is to investigate the relative importance of four dimensionless parameters: the Reynolds number, the Froude number, the Tachikawa number, and the mass moment of inertia parameter. A secondary objective of the paper is to study the effect of launching time on the trajectories. The focus is on the start of the motion, when aerodynamic forces on the ice fragment are high. A trajectory database is built for a 2D rectangular plate at Reynold numbers around 7000, but with other parameter values typical of inflight aircraft deicing. The plate makes an initial angle with the flow field such that flow separation occurs.The CFD solver uses an immersed boundary method (IBM) to model the flow around the plate geometry [17]. A Cartesian grid is used for the computational domain. The plate boundaries are modelled by adding continuous forcing terms directly to the flow equations using the penalization method [18–20].The paper is organized as follows. First, the penalized Navier-Stokes equations are presented in Section2, together with the numerical method based on a Vortex-In-Cell (VIC) scheme. The model used to compute the forces and moments on the plate, along with the appropriate governing equations for fluid-solid interactions, is also described. In Section 3, test cases are presented to validate computations of forces and moments against literature results. In Section 4, a database is created to study the influences of the four parameters and initial angle of attack on the trajectories of a 2D rectangular plate.


## 2. Mathematical Model and Numerical Method

The fluid-solid interaction model consists of a set of equations for the fluid flow and another set for the solid motion. The effects of the flow on the solid motion are imposed through aerodynamic forces and moments, while those of the solid on the flow are imposed by penalization.


### 2.1. Fluid Motion

The incompressible laminar Navier-Stokes equations model the flow into the computational domainΩ. A rigid solid Si, moving at velocity usi in the domain Ω, is defined with a level set function Φsi. In this work, Φsi is the signed distance function to Si, negative inside the solid and positive outside. The boundaries of Si are located where the sign of the level set function changes from positive to negative. Figure 1 illustrates the signed distance function around the plate. A Cartesian grid is used to discretize the domain Ω.Figure 1
Global level set function of one rigid body (rectangular plate) onto a Cartesian grid.
(a)
(b)The velocity field is calculated throughout the computational domain, around and inside the solid, as illustrated in Figure2. The penalization technique is used to impose the solid velocity vector usi in the area of the computational domain identified by the negative Φsi [18].Figure 2
The penalization technique: the flow field is extended inside the body, the instantaneous vorticity, and the streamlines around a plate.For a very large penalization parameter,λ≫1, and the Heavyside function, H(Φsi), that is, H(Φsi)=1 if Φsi<0 and H(Φsi)=0 otherwise, the penalized Navier-Stokes equations in a vorticity formulation are [10](1)∂ω∂t+u·∇ω=ω·∇u+ν∇2ω+λ∇×HΦsiusi-u(2)with∇·u=0inΩ,where u is the velocity vector, ν=μ/ρ is the kinematic viscosity, ρ is the fluid density, and ω=∇×u is the vorticity vector.The incompressible flow has a divergence-free velocity field. To impose boundary conditions and penalization inside the body, the velocity can be recovered from the vector potentialΨ:(3)u=∇×Ψ.The vector potentialΨ is a 3D extension of the so-called scalar stream function Ψ. This potential vector is imposed to be solenoidal; that is, ∇·Ψ=0. By definition, the stream function field is related to ω by(4)ΔΨ=-ω.


### 2.2. Aerodynamic Forces and Moments

The penalization term in (1), λ∇×H(Φsi)(usi-u), forces the fluid velocity inside the rigid body Si to be nearly equal to usi. The local acceleration imposed by the penalization is computed by(5)∂u∂t=λHΦsiusi-u.The penalization term can never be zero, because an acceleration is always needed to counteract the pressure and viscous forces from the Navier-Stokes equation. Thus,(6)u=usi-1λHΦsi∂u∂tand a large penalization parameter λ is needed to have u≈usi. The sensitivity of the solution to the choice of the penalization parameter is discussed in Section 2.4.The aerodynamic forces and moment are obtained by integration over the computational solid domainSi:(7)F=∫SiρλHΦsiusi-udx,T=∫SiρλHΦsir×usi-udx.In 2D flows, the moment becomes a scalar.


### 2.3. Solid Motion

For a fluid-solid interaction the forceF and moment T induce the solid motion:(8)∂usi∂t=FMa+ρsi-ρρsiG,∂θ˙∂t=TI,where Ma is the mass of the solid Si, I its moment of inertia, ρsi its density, L is the characteristic length, θ˙ the angular velocity, and G=0-g the gravity vector.Following [5–7], the dimensionless variables x′=(g/U∞2)x, y′=(g/U∞2)y, and t′=(g/U∞)t are introduced, U∞ being the freestream velocity. In the resulting dimensionless equations of motion, three parameters have significant effects on the solid motion: the Tachikawa number K (ratio of the aerodynamic force to the gravity force), the Froude number Fr (or dimensionless plate length Ln=1/Fr2), and the dimensionless mass moment of inertia parameter In of the solid:(9)K=ρAU∞22Mag,Fr=U∞gL,In=IMaL2,where A is the reference area of the solid.For a flat plate of thicknessh in 2D motion, Ma=ρsiAh, and the mass moment of inertia is(10)I=MaL2+h212.Thus, In=1/12(1+(h/L)2) and K=(ρU∞2)/(2ρsihg). The dimensionless mass moment of inertia parameter depends solely on the aspect ratio of the plate h/L.The motion equations in dimensionless form are(11)∂u′si∂t′=KF′+ρsi-ρρsi0-1,(12)∂θ˙′∂t′=KFr2T′Inwith(13)u′si=usiU∞,F′=F1/2ρU∞2L,T′=T1/2ρU∞2L2.For the particular case of a solid in motion in gas, most oftenρ/ρsi≪1. Therefore, the solid motion depends on three parameters, separately from the dimensionless forces and moments. Note that, contrary to previous works [5–7], the freestream velocity U∞ is used for dimensionless forces and moments instead of the relative velocity.


### 2.4. Fluid Structure Interaction with Penalized VIC Scheme

The domainΩ is meshed using a uniform fixed Cartesian grid. Equation (1) is solved using a classical VIC scheme [21], detailed in [10]. The grid values of Φsin, un, and ωn are evaluated at time tn=nΔt. Each time step Δt is solved using three substeps as follows:(1)
Advection:(14a)ω⋆-ωΔt=-u·∇ω.(2)
Stretching and diffusion:(14b)ω⋆⋆-ω⋆Δt=ω⋆·∇u+ν∇2ω⋆.(3)
Penalization term:(14c)ωn+1=∇×u⋆⋆+λΔtHΦsiusin1+λΔtHΦsi.In the substep(1), the vortex particles located at grid point are displaced with a fourth-order Runge-Kutta time-stepping scheme and remeshed on the grid. In substep (2), the Laplacian is evaluated on the grid with a second-order accurate standard five-point stencil. Then, the velocity u⋆⋆ is computed using ω⋆⋆ (ΔΨ⋆⋆=-ω⋆⋆ and u⋆⋆=∇×Ψ⋆⋆). The stream function is computed by solving the linear Poisson equation on the Cartesian grid, using the fast Poisson solver routine of the Intel® Math Kernel Library.To assess the force acting on the solid, the local acceleration is computed as follows:(15)∂u∂t=λHΦsiusin-u⋆⋆+λΔtHΦsiusin1+λΔtHΦsi,∂u∂t=λHΦsiusin-u⋆⋆1+λΔtHΦsi.The aerodynamic forces and moment are then defined by(16)F=∫SiρλHΦsi1+λΔtHΦsiusin-u⋆⋆dx,T=∫SiρλHΦsi1+λΔtHΦsir×usin-u⋆dx.A first-order explicit Euler scheme is applied for the time integration to update the solid linear and angular velocities. The solid velocities are then used to displace the solid boundary during the time intervalΔt. From the displaced solid position, Φsi is computed prior to the next time step. The solid boundary is defined by line segments and nodes. For each mesh node, the shortest distance between the node and the displaced line segments becomes the new Φsi value. A node is inside the geometry, negative Φsi, if a vertical line drawn from the node toward the upper limit of the computational domain crosses an uneven times the solid boundary line segments. The level set function is recalculated at each time step, thus avoiding the possible solid deformations associated with the advection of Φsi. In 2D, the computational cost of the recalculation is small compared to the cost of the fluid solution.The penalization term depends on theλ values. The case of a 2D cylinder in a square cavity, falling under gravity on a flat plane, is used to study the sensitivity of the solution to the choice of these values. This test case has been presented in detail and used for validation in a previous article [22]. The dimension of the cavity is [0,2]×[0,6]. 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 and is initially located at the point (1,4). It accelerates under gravity, set to g=-980. Figures 3(a) and 3(b) show the evolution of the cylinder vertical velocity Vy with time. The parameter values range from 0.1/Δt to 108/Δt with a tenfold increase between values. For values of 10/Δt, a close-up view on velocity values near t=0.3 is provided on Figure 3(b) to see the difference. There is 0.17% of difference between velocities computed with 104/Δt and 108/Δt. All the calculations in this paper are done with λ=108/Δt.Figure 3
Sensitivity of the cylinder falling velocity toλ values.
(a)
(b)


## 2.1. Fluid Motion

The incompressible laminar Navier-Stokes equations model the flow into the computational domainΩ. A rigid solid Si, moving at velocity usi in the domain Ω, is defined with a level set function Φsi. In this work, Φsi is the signed distance function to Si, negative inside the solid and positive outside. The boundaries of Si are located where the sign of the level set function changes from positive to negative. Figure 1 illustrates the signed distance function around the plate. A Cartesian grid is used to discretize the domain Ω.Figure 1
Global level set function of one rigid body (rectangular plate) onto a Cartesian grid.
(a)
(b)The velocity field is calculated throughout the computational domain, around and inside the solid, as illustrated in Figure2. The penalization technique is used to impose the solid velocity vector usi in the area of the computational domain identified by the negative Φsi [18].Figure 2
The penalization technique: the flow field is extended inside the body, the instantaneous vorticity, and the streamlines around a plate.For a very large penalization parameter,λ≫1, and the Heavyside function, H(Φsi), that is, H(Φsi)=1 if Φsi<0 and H(Φsi)=0 otherwise, the penalized Navier-Stokes equations in a vorticity formulation are [10](1)∂ω∂t+u·∇ω=ω·∇u+ν∇2ω+λ∇×HΦsiusi-u(2)with∇·u=0inΩ,where u is the velocity vector, ν=μ/ρ is the kinematic viscosity, ρ is the fluid density, and ω=∇×u is the vorticity vector.The incompressible flow has a divergence-free velocity field. To impose boundary conditions and penalization inside the body, the velocity can be recovered from the vector potentialΨ:(3)u=∇×Ψ.The vector potentialΨ is a 3D extension of the so-called scalar stream function Ψ. This potential vector is imposed to be solenoidal; that is, ∇·Ψ=0. By definition, the stream function field is related to ω by(4)ΔΨ=-ω.


## 2.2. Aerodynamic Forces and Moments

The penalization term in (1), λ∇×H(Φsi)(usi-u), forces the fluid velocity inside the rigid body Si to be nearly equal to usi. The local acceleration imposed by the penalization is computed by(5)∂u∂t=λHΦsiusi-u.The penalization term can never be zero, because an acceleration is always needed to counteract the pressure and viscous forces from the Navier-Stokes equation. Thus,(6)u=usi-1λHΦsi∂u∂tand a large penalization parameter λ is needed to have u≈usi. The sensitivity of the solution to the choice of the penalization parameter is discussed in Section 2.4.The aerodynamic forces and moment are obtained by integration over the computational solid domainSi:(7)F=∫SiρλHΦsiusi-udx,T=∫SiρλHΦsir×usi-udx.In 2D flows, the moment becomes a scalar.


## 2.3. Solid Motion

For a fluid-solid interaction the forceF and moment T induce the solid motion:(8)∂usi∂t=FMa+ρsi-ρρsiG,∂θ˙∂t=TI,where Ma is the mass of the solid Si, I its moment of inertia, ρsi its density, L is the characteristic length, θ˙ the angular velocity, and G=0-g the gravity vector.Following [5–7], the dimensionless variables x′=(g/U∞2)x, y′=(g/U∞2)y, and t′=(g/U∞)t are introduced, U∞ being the freestream velocity. In the resulting dimensionless equations of motion, three parameters have significant effects on the solid motion: the Tachikawa number K (ratio of the aerodynamic force to the gravity force), the Froude number Fr (or dimensionless plate length Ln=1/Fr2), and the dimensionless mass moment of inertia parameter In of the solid:(9)K=ρAU∞22Mag,Fr=U∞gL,In=IMaL2,where A is the reference area of the solid.For a flat plate of thicknessh in 2D motion, Ma=ρsiAh, and the mass moment of inertia is(10)I=MaL2+h212.Thus, In=1/12(1+(h/L)2) and K=(ρU∞2)/(2ρsihg). The dimensionless mass moment of inertia parameter depends solely on the aspect ratio of the plate h/L.The motion equations in dimensionless form are(11)∂u′si∂t′=KF′+ρsi-ρρsi0-1,(12)∂θ˙′∂t′=KFr2T′Inwith(13)u′si=usiU∞,F′=F1/2ρU∞2L,T′=T1/2ρU∞2L2.For the particular case of a solid in motion in gas, most oftenρ/ρsi≪1. Therefore, the solid motion depends on three parameters, separately from the dimensionless forces and moments. Note that, contrary to previous works [5–7], the freestream velocity U∞ is used for dimensionless forces and moments instead of the relative velocity.


## 2.4. Fluid Structure Interaction with Penalized VIC Scheme

The domainΩ is meshed using a uniform fixed Cartesian grid. Equation (1) is solved using a classical VIC scheme [21], detailed in [10]. The grid values of Φsin, un, and ωn are evaluated at time tn=nΔt. Each time step Δt is solved using three substeps as follows:(1)
Advection:(14a)ω⋆-ωΔt=-u·∇ω.(2)
Stretching and diffusion:(14b)ω⋆⋆-ω⋆Δt=ω⋆·∇u+ν∇2ω⋆.(3)
Penalization term:(14c)ωn+1=∇×u⋆⋆+λΔtHΦsiusin1+λΔtHΦsi.In the substep(1), the vortex particles located at grid point are displaced with a fourth-order Runge-Kutta time-stepping scheme and remeshed on the grid. In substep (2), the Laplacian is evaluated on the grid with a second-order accurate standard five-point stencil. Then, the velocity u⋆⋆ is computed using ω⋆⋆ (ΔΨ⋆⋆=-ω⋆⋆ and u⋆⋆=∇×Ψ⋆⋆). The stream function is computed by solving the linear Poisson equation on the Cartesian grid, using the fast Poisson solver routine of the Intel® Math Kernel Library.To assess the force acting on the solid, the local acceleration is computed as follows:(15)∂u∂t=λHΦsiusin-u⋆⋆+λΔtHΦsiusin1+λΔtHΦsi,∂u∂t=λHΦsiusin-u⋆⋆1+λΔtHΦsi.The aerodynamic forces and moment are then defined by(16)F=∫SiρλHΦsi1+λΔtHΦsiusin-u⋆⋆dx,T=∫SiρλHΦsi1+λΔtHΦsir×usin-u⋆dx.A first-order explicit Euler scheme is applied for the time integration to update the solid linear and angular velocities. The solid velocities are then used to displace the solid boundary during the time intervalΔt. From the displaced solid position, Φsi is computed prior to the next time step. The solid boundary is defined by line segments and nodes. For each mesh node, the shortest distance between the node and the displaced line segments becomes the new Φsi value. A node is inside the geometry, negative Φsi, if a vertical line drawn from the node toward the upper limit of the computational domain crosses an uneven times the solid boundary line segments. The level set function is recalculated at each time step, thus avoiding the possible solid deformations associated with the advection of Φsi. In 2D, the computational cost of the recalculation is small compared to the cost of the fluid solution.The penalization term depends on theλ values. The case of a 2D cylinder in a square cavity, falling under gravity on a flat plane, is used to study the sensitivity of the solution to the choice of these values. This test case has been presented in detail and used for validation in a previous article [22]. The dimension of the cavity is [0,2]×[0,6]. 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 and is initially located at the point (1,4). It accelerates under gravity, set to g=-980. Figures 3(a) and 3(b) show the evolution of the cylinder vertical velocity Vy with time. The parameter values range from 0.1/Δt to 108/Δt with a tenfold increase between values. For values of 10/Δt, a close-up view on velocity values near t=0.3 is provided on Figure 3(b) to see the difference. There is 0.17% of difference between velocities computed with 104/Δt and 108/Δt. All the calculations in this paper are done with λ=108/Δt.Figure 3
Sensitivity of the cylinder falling velocity toλ values.
(a)
(b)


## 3. Results and Discussion

### 3.1. Lescape’s Fluid-Solid Model Validation and Verification

The computations of forces, using the code Lescape (LEvel Set CArtesian PEnalization), on a static body have already been validated in a previous work [23]. In this section, computations of forces and moments for solid motions are validated using an oscillating cylinder and a flapping airfoil.


#### 3.1.1. Inline Oscillating Cylinder in a Fluid at Rest

The interaction of an oscillating circular cylinder with a fluid at rest is a problem well-documented in the literature [13, 14, 24–27]. The two key parameters when dealing with such flows are the Reynolds number Re=UmaxD/ν and the Keulegan-Carpenter number, KC=Umax/fD, where Umax is the maximum velocity of the cylinder, D is the diameter of the cylinder, and f is the characteristic frequency of the oscillation. The set-up for this computation corresponds to the LDA experimental and numerical simulations performed by [13].The motion of the cylinder is described by a simple harmonic oscillation,(17)xt=-H0sin⁡2πft,where H0 is the amplitude of the oscillation. The parameters were set to Re=100 and KC=5. The computational domain is 50D and 30D in the x and y directions, with the cylinder initially located at the center of the domain. The grid contains 120×120 cells in the cylinder. Simulations are started in a static fluid field, and the cylinder moves according to (17) until periodic vortex shedding is established (see Figure 4 for an overview of the vortex shedding). Quantitative comparisons with the experimental and numerical data from [13] are considered in Figure 5, which shows the computed velocity profiles at four x locations and three different phase angles. The numerical results from [13] were obtained by solving the Navier-Stokes equations for the velocity and pressure with a 2D finite volume method imposing an oscillatory flow around a cylinder at rest. The maximum difference between the numerical results is below 5% for both u and v velocity, which indicates a good agreement, taking into account the differences in the mathematical models and the numerical methods. In the paper by [13], the comparison between the numerical and the experimental data is considered to be very satisfactory, even if the experimental error is not given. It is expected that the velocity distributions should be symmetric for u and antisymmetric for v. The dimensional drag force divided by the cylinder diameter over one period is also in good agreement with literature results (see Figure 6).Figure 4
Vorticity contours.
(a)
(b)
(c)
(d)Figure 5
Comparison with results taken of velocity profiles at four cross sections from [13] with constant x locations; phase angles: (a-b) 180°, (c-d) 210°, and (e-f) 330°.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Comparison of the drag force obtained using the present CFD code, Lescape, with [14] and experimental results taken from [13].


#### 3.1.2. Validation of Force Calculations through an Oscillating Airfoil

An oscillating NACA0015 airfoil simultaneously experiencing pitching θ(t) and heaving h(t) motions is modelled. The pitching axis is located along the airfoil chord at position (xp,yp)=(1/3,0). The airfoil motion, described by [15], is defined as follows:(18)ht=H0sin⁡2πft+φ,θt=θ0sin⁡2πft.θ0 is the pitching amplitude, and the phase difference φ is set to 90°. The heaving velocity is then given by (19)Vyt=H02πfcos⁡2πft+φ.Based on the imposed motion and the upstream flow conditions, the airfoil experiences an effective angle of attackα(t) and an effective upstream velocity Veff(t) defined by(20)αt=arctan-VytU∞-θt,Vefft=U∞2+Vy2t,where the freestream velocity far upstream of the oscillating airfoil is U∞=1.To validate our simulations, a regime corresponding to the parametersRe=U∞c/ν=1100, H0/c=1, f=0.14, xp/c=1/3, and θ0=76.33° was computed. A view of the motion is sketched in Figure 7.Figure 7
Oscillating airfoil, sketch of the airfoil motion forH0/c=1, f=0.14, xp/c=1/3, and θ0=76.33°.The computational domain is[-14:16]×[-15:15]; the boundary conditions are set to Dirichlet on the left, top, and bottom, and a Neumann boundary condition is set on the right. The numerical results, obtained using dx=dy=5×10-3 and dt=1×10-3, are then compared to the force predictions presented by [15, 16].The force coefficientsCX, CY, and the moment coefficients Cm, are presented in Figure 8. The drag, lift, and moment coefficients are defined, respectively, as(21)CX=Fx1/2ρU∞2L,CY=Fy1/2ρU∞2L,Cm=T1/2ρU∞2L2,where the x-axis is aligned with U∞ and the y-axis is perpendicular to the far field velocity vector. Considering the Lescape’s drag coefficient (Figure 8(a)), the time evolution of the drag coefficient coincides with the results of Kinsey et al. and Campobasso et al. The numerical error for CX is below 3% for the maximum and minimum amplitudes. Lescape slightly overpredicts (by around 10%) CX during the decreasing phase of the drag. The evolution of the lift (Figure 8(b)) coincides with literature results. Lescape CY is in phase advance compared to Kinsey et al. and Campobasso et al., and the maximum time lag is below 5% just before the mid cycle T/2. The error in amplitude is globally below 5% during the cycle, with the maximum error occurring at the end of the cycle, and being about 10%. The pitching moment (Figure 8(c)) is the one with the least agreement with the literature results, with a maximal error of 20% at t/T=0.3 and t/T=0.8, the remaining error being below 10%. It must be noticed that, for pitching moments, Lescape shows a tendency to underestimate maximum (positive values) and overestimate minimum (negative values). The evolution in time of Cm follows the literature results; it increases when literature results increase and decreases in accordance with published results (meaning that, in the situation of fluid-solid interaction, the Lescape code will predict the rotation angle and its sign to an acceptable extent). For Cm, a small time lag at a maximum of 3% is observed. CX, CY, and Cm numerical results are satisfactory in the perspective of trajectory computations.Figure 8
Oscillating airfoil, comparison with the results obtained by [15, 16]. (a) Drag coefficient, (b) lift coefficient, and (c) pitching moment.
(a)
(b)
(c)


### 3.2. Plate Trajectories Database

For the parametric study, a sharp corner plate is used, as in Figure2. During these simulations, the rotation axis of the rectangular plate coincides with its barycenter. The effects of K, Fr, In, and Re on flat plate trajectories are studied in order to identify the important parameters at the start of the trajectory computation. A full factorial design with four factors and two levels [28] is used to build the test matrix.The simulations are performed usingdx=dy=1.×10-3 and a domain of [-3,8]×[-3,5], with the plate’s center initially located at (0,0), and with starting angles of ±π/3. For each factor, minimal and maximal values are selected around a base value: Re=7000, K=30, Fr=57.7, and In=8.44×10-2. Except for the Reynolds number, these are the typical values for shed ice pieces around an aircraft in a deicing situation [29]. For In, two values of the aspect ratio h/l are used, namely, 0.05 and 0.15. This is only a 1% variation, but this is a significant thickness change. Table 1 gives the 16 run factors. The 16 runs are executed for two plate starting angles and two release times. Thus, the database consists of 64 trajectories.Table 1
Parameters for numerical experiments.
Run number
Fr
Re
K
I
n
(
1
)
51.7
8750
37.5
8.52
×
1
0
-
2
(
2
)
51.7
8750
37.5
8.35
×
1
0
-
2
(
3
)
51.7
8750
22.5
8.52
×
1
0
-
2
(
4
)
51.7
8750
22.5
8.35
×
1
0
-
2
(
5
)
51.7
5250
37.5
8.52
×
1
0
-
2
(
6
)
51.7
5250
37.5
8.35
×
1
0
-
2
(
7
)
51.7
5250
22.5
8.52
×
1
0
-
2
(
8
)
51.7
5250
22.5
8.35
×
1
0
-
2
(
9
)
66.6
8750
37.5
8.52
×
1
0
-
2
(
10
)
66.6
8750
37.5
8.35
×
1
0
-
2
(
11
)
66.6
8750
22.5
8.52
×
1
0
-
2
(
12
)
66.6
8750
22.5
8.35
×
1
0
-
2
(
13
)
66.6
5250
37.5
8.52
×
1
0
-
2
(
14
)
66.6
5250
37.5
8.35
×
1
0
-
2
(
15
)
66.6
5250
22.5
8.52
×
1
0
-
2
(
16
)
66.6
5250
22.5
8.35
×
1
0
-
2The starting angle creates a separated flow around the plate. Vortices are shed periodically at low frequencies. The flow is first solved around a static plate, which goes on until the periodic flow is well established. The plate is released into the flow either at timet=20.0 s or at time t=21.1 s. The dimensionless position Y′ of the plate when it crosses the X′=1.125×10-3 line will be studied.The computed trajectories are illustrated in Figures9 and 10. To help visualize the individual trajectories, the results are separated into four. For a positive plate starting angle, these results are presented in Figure 9, with a further separation between the two shedding times. Similarly, the results for a negative plate starting angle are shown in Figure 10. The positions X and Y are plotted against each other. Since the initial angle of the plate is not null, vortex shedding is created behind the plate. The influence of vortex shedding is not negligible for a particular trajectory. For different starting times, the plate positions at the end of the trajectory calculations are different, as comparison between (a) and (b) shows. However, the results at different launching time will be grouped together to study the effects of other parameters. It can be seen that, in general, the plate accelerates faster in the X direction than in the Y direction. The gravity effects curve the trajectory. Here, the plate motion is mostly driven by aerodynamic forces and the starting angle has a dominant effect.Figure 9
Dimensional trajectories representative of the test matrix for a positive starting angle of the plate; Red Fr number51.7; Blue Fr number 66.6; square Re number 8750; circle Re number 5250; filled symbolsK37.5; empty symbolsK22.5; line inertia parameter 8.52×10-2; dashed line inertia parameter 8.35×10-2.
(a)
(b)Figure 10
Dimensional trajectories representative of the test matrix for a negative starting angle of the plate; Red Fr number51.7; Blue Fr number 66.6; square Re number 8750; circle Re number 5250; filled symbolsK37.5; empty symbolsK22.5; line inertia parameter 8.52×10-2; dashed line inertia parameter 8.35×10-2.
(a)
(b)Figure11 presents an example of the trajectories obtained for a plate with an initial positive angle. The dimensionless positions of the center of gravity for the 16 run parameters with an initial positive angle of attacks and a release at t=20.0 s is shown. The computed average trajectory of the 16 runs is also plotted, together with the upper and lower limits for a confidence interval. There is a 90% probability that the mean position of the plate gravity center at a given X′ position lies between the upper and lower limits. The confidence interval is computed using the sample variance and a Student probability distribution [30]. As can be seen, many trajectories are clearly outside the confidence interval. The mean trajectory gives a representative idea of the plate motion, but all the trajectory must be plotted to take into account the stochastic nature of the motion.Figure 11
Mean trajectory and individual trajectories for the16 runs with an initial positive angle of attack for the plate.The plate geometry, corresponding to run(1), is drawn at a constant time interval in Figure 12. For this run, the plate is rotating only slightly. The horizontal acceleration of the plate is visible as the spacing between its locations increases with the X′ position. The mean trajectory of the center of gravity is also plotted with the corresponding confidence interval. For run (1), the average trajectory is representative of the plate motion.Figure 12
Mean trajectory and plate orientation for run(1) for a plate with an initial negative angle of attack.The trajectory footprint probability maps used for the 3D trajectories analysis by [3] become histograms for 2D trajectories. At X′=1.125, a line is drawn in the Y direction and split into intervals. The heights of the histogram bar are computed by dividing the count of plate passage into an interval ΔY′ by 32, which represents the total number of runs [28]. The histograms associated with the two initial plate angular orientations are shown in Figure 13. In this figure, the Y′ position is on the horizontal axis, and two histograms are plotted, one for each plate orientation. The length of each bar gives the probability of passage of the plate gravity center in an interval of ΔY′=0.4×10-4. The initial plate orientation has a clear effect on the histograms, with the positive angle giving positive Y′ values, and the negative angle giving negative Y′ values. This is expected since an initial positive angle creates a positive lift coefficient, and a negative angle creates a negative lift coefficient. The two histograms are placed almost symmetrically around the Y′=0 line, indicating that gravity is not the dominant force of the problem at the beginning of the motion.Figure 13
Histogram of the probability of passage atY′ for X′=1.125×10-3.The initial plate orientation is the parameter that has the most effect on the trajectories. To study the effects of the other parameters, the database is split into two sets of32 runs: set one for the starting positive angle, and set two for the starting negative angle. Box plots are used to compare parameter effects [28]. Figure 14 shows the relation between the box plot and the probability histogram for the results of the flat plate with a starting negative angle. The central mark is the median of the 32 runs. The edges of the box are Q1, the 25th percentile, and Q3, the 75th percentile. The whiskers extend to Q1-1.5(Q3-Q1) and Q3+1.5(Q3-Q1), or the most extreme data points not considered outliers. The points for which Y′ are smaller than Q1-1.5(Q3-Q1) or greater than Q3+1.5(Q3-Q1) are considered outliers. One outlier is visible in Figure 14.Figure 14
The box plot and the corresponding histogram for the32 runs of the flat plate with a negative initial angle.For the trajectories computed with a starting angle ofπ/3 (positive orientation) (Figure 15), the box plots show that the Fr parameter has a significant effect on the median value. Also, the median value for K- is lower than the 25th percentile for the K+ distribution, which confirms that gravity forces are more important for lower K values. For the two other dimensionless parameters, the median does not change significantly between minimum and maximum values. For example, for Re+, the median is between the edges of the box for Re- distribution, while that for Re- is between the edges of the box for Re+ distribution. The result dispersion is also changed by Fr values and increases with increasing Fr values. The whiskers extend from 3×10-4 to 6.7×10-4 for Fr+ and from 4.3×10-4 to 6.5×10-4 for Fr-. According to the box plot, lower Fr values will make the plate fly higher and will reduce the dispersion. As Fr is only present in (12), the result dispersion is caused by the plate rotation. We should recall that Fr is the dimensionless length of the plate; Ln=1/Fr2. It means that longer dimensionless length plates fly more easily and in a more predictable fashion.Figure 15
Boxplot of theY′ position at X′=1.125×10-3 for a plate with a starting positive angle orientation.The box plots for a plate at a starting angle of-π/3 (negative orientation) are shown in Figure 16. The Fr number once again has the most significant effects on the plate location Y′ at X′=1.125×10-3. The standard deviation increases noticeably with an increase of the Fr number. Decreasing the K decreases the median Y′ position of the plate, as was the case for the plate with a starting positive angle. Decreasing Re and In has almost no effects on the final center of gravity position.Figure 16
Boxplot of theY′ position at X′=1.125×10-3 for a plate with a starting negative angle orientation.This database parametric study focuses on the start of the trajectory. To give an order of magnitude, in typical ice trajectory experiments [29], the ice piece travels a dimensionless horizontal distance of 4.×10-3, three times the horizontal distance in this database. The major limitation of this study is the 2D computational domain that precludes any comparisons against 3D experimental results. The computed 2D results should be equivalent to 3D results with a plate of infinite span and different from the results obtained with a plate of short span extent. Short span plates are typical of most engineering applications and experiments. Another limitation is the small number of runs used for the statistical analysis. Although 64 runs are used to build the database, the effect of the plate starting angle on the final position is strong, and runs must be separated into two databases. Trajectory depends on launch time and multiple launch times should be used to populate a larger database. Nevertheless, the database size is comparable to the database size used by [7] for parametric study.The numerical errors on aerodynamic coefficients observed during the oscillating airfoil validation could also have an impact on the database results. The aerodynamic moment coefficient oscillates typically between −0.6 and 0.6 during one trajectory calculation. The number of oscillation cycles is usually above three. Assuming the airfoil validation errors apply for the moving plate, the moment coefficient maximum values are underestimated by20% and the minimum values overestimated by 20%. For one cycle, the numerical errors have probably neutral effects. Also, it is worth mentioning that, for Lagrangian methods, often used for trajectory calculations, 10% of errors on static aerodynamic coefficients are usually assumed, together with an unknown uncertainty on the damping function used to model the dynamic caused by rotational velocity [5]. An additional limitation specific to aeronautic is due to the values of the Re, which are one level of magnitude lower than those for aerodynamic applications. This may not be a critical limitation as the static experimental lift, drag, and momentum correlation for a rectangular plate do not show any dependence on Re [31]. However, Lagrangian trajectory calculations use corrections to the lift, drag, and momentum that depend on Reynolds. Actual CFD results show some Re number effects.


## 3.1. Lescape’s Fluid-Solid Model Validation and Verification

The computations of forces, using the code Lescape (LEvel Set CArtesian PEnalization), on a static body have already been validated in a previous work [23]. In this section, computations of forces and moments for solid motions are validated using an oscillating cylinder and a flapping airfoil.


### 3.1.1. Inline Oscillating Cylinder in a Fluid at Rest

The interaction of an oscillating circular cylinder with a fluid at rest is a problem well-documented in the literature [13, 14, 24–27]. The two key parameters when dealing with such flows are the Reynolds number Re=UmaxD/ν and the Keulegan-Carpenter number, KC=Umax/fD, where Umax is the maximum velocity of the cylinder, D is the diameter of the cylinder, and f is the characteristic frequency of the oscillation. The set-up for this computation corresponds to the LDA experimental and numerical simulations performed by [13].The motion of the cylinder is described by a simple harmonic oscillation,(17)xt=-H0sin⁡2πft,where H0 is the amplitude of the oscillation. The parameters were set to Re=100 and KC=5. The computational domain is 50D and 30D in the x and y directions, with the cylinder initially located at the center of the domain. The grid contains 120×120 cells in the cylinder. Simulations are started in a static fluid field, and the cylinder moves according to (17) until periodic vortex shedding is established (see Figure 4 for an overview of the vortex shedding). Quantitative comparisons with the experimental and numerical data from [13] are considered in Figure 5, which shows the computed velocity profiles at four x locations and three different phase angles. The numerical results from [13] were obtained by solving the Navier-Stokes equations for the velocity and pressure with a 2D finite volume method imposing an oscillatory flow around a cylinder at rest. The maximum difference between the numerical results is below 5% for both u and v velocity, which indicates a good agreement, taking into account the differences in the mathematical models and the numerical methods. In the paper by [13], the comparison between the numerical and the experimental data is considered to be very satisfactory, even if the experimental error is not given. It is expected that the velocity distributions should be symmetric for u and antisymmetric for v. The dimensional drag force divided by the cylinder diameter over one period is also in good agreement with literature results (see Figure 6).Figure 4
Vorticity contours.
(a)
(b)
(c)
(d)Figure 5
Comparison with results taken of velocity profiles at four cross sections from [13] with constant x locations; phase angles: (a-b) 180°, (c-d) 210°, and (e-f) 330°.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Comparison of the drag force obtained using the present CFD code, Lescape, with [14] and experimental results taken from [13].


### 3.1.2. Validation of Force Calculations through an Oscillating Airfoil

An oscillating NACA0015 airfoil simultaneously experiencing pitching θ(t) and heaving h(t) motions is modelled. The pitching axis is located along the airfoil chord at position (xp,yp)=(1/3,0). The airfoil motion, described by [15], is defined as follows:(18)ht=H0sin⁡2πft+φ,θt=θ0sin⁡2πft.θ0 is the pitching amplitude, and the phase difference φ is set to 90°. The heaving velocity is then given by (19)Vyt=H02πfcos⁡2πft+φ.Based on the imposed motion and the upstream flow conditions, the airfoil experiences an effective angle of attackα(t) and an effective upstream velocity Veff(t) defined by(20)αt=arctan-VytU∞-θt,Vefft=U∞2+Vy2t,where the freestream velocity far upstream of the oscillating airfoil is U∞=1.To validate our simulations, a regime corresponding to the parametersRe=U∞c/ν=1100, H0/c=1, f=0.14, xp/c=1/3, and θ0=76.33° was computed. A view of the motion is sketched in Figure 7.Figure 7
Oscillating airfoil, sketch of the airfoil motion forH0/c=1, f=0.14, xp/c=1/3, and θ0=76.33°.The computational domain is[-14:16]×[-15:15]; the boundary conditions are set to Dirichlet on the left, top, and bottom, and a Neumann boundary condition is set on the right. The numerical results, obtained using dx=dy=5×10-3 and dt=1×10-3, are then compared to the force predictions presented by [15, 16].The force coefficientsCX, CY, and the moment coefficients Cm, are presented in Figure 8. The drag, lift, and moment coefficients are defined, respectively, as(21)CX=Fx1/2ρU∞2L,CY=Fy1/2ρU∞2L,Cm=T1/2ρU∞2L2,where the x-axis is aligned with U∞ and the y-axis is perpendicular to the far field velocity vector. Considering the Lescape’s drag coefficient (Figure 8(a)), the time evolution of the drag coefficient coincides with the results of Kinsey et al. and Campobasso et al. The numerical error for CX is below 3% for the maximum and minimum amplitudes. Lescape slightly overpredicts (by around 10%) CX during the decreasing phase of the drag. The evolution of the lift (Figure 8(b)) coincides with literature results. Lescape CY is in phase advance compared to Kinsey et al. and Campobasso et al., and the maximum time lag is below 5% just before the mid cycle T/2. The error in amplitude is globally below 5% during the cycle, with the maximum error occurring at the end of the cycle, and being about 10%. The pitching moment (Figure 8(c)) is the one with the least agreement with the literature results, with a maximal error of 20% at t/T=0.3 and t/T=0.8, the remaining error being below 10%. It must be noticed that, for pitching moments, Lescape shows a tendency to underestimate maximum (positive values) and overestimate minimum (negative values). The evolution in time of Cm follows the literature results; it increases when literature results increase and decreases in accordance with published results (meaning that, in the situation of fluid-solid interaction, the Lescape code will predict the rotation angle and its sign to an acceptable extent). For Cm, a small time lag at a maximum of 3% is observed. CX, CY, and Cm numerical results are satisfactory in the perspective of trajectory computations.Figure 8
Oscillating airfoil, comparison with the results obtained by [15, 16]. (a) Drag coefficient, (b) lift coefficient, and (c) pitching moment.
(a)
(b)
(c)


## 3.1.1. Inline Oscillating Cylinder in a Fluid at Rest

The interaction of an oscillating circular cylinder with a fluid at rest is a problem well-documented in the literature [13, 14, 24–27]. The two key parameters when dealing with such flows are the Reynolds number Re=UmaxD/ν and the Keulegan-Carpenter number, KC=Umax/fD, where Umax is the maximum velocity of the cylinder, D is the diameter of the cylinder, and f is the characteristic frequency of the oscillation. The set-up for this computation corresponds to the LDA experimental and numerical simulations performed by [13].The motion of the cylinder is described by a simple harmonic oscillation,(17)xt=-H0sin⁡2πft,where H0 is the amplitude of the oscillation. The parameters were set to Re=100 and KC=5. The computational domain is 50D and 30D in the x and y directions, with the cylinder initially located at the center of the domain. The grid contains 120×120 cells in the cylinder. Simulations are started in a static fluid field, and the cylinder moves according to (17) until periodic vortex shedding is established (see Figure 4 for an overview of the vortex shedding). Quantitative comparisons with the experimental and numerical data from [13] are considered in Figure 5, which shows the computed velocity profiles at four x locations and three different phase angles. The numerical results from [13] were obtained by solving the Navier-Stokes equations for the velocity and pressure with a 2D finite volume method imposing an oscillatory flow around a cylinder at rest. The maximum difference between the numerical results is below 5% for both u and v velocity, which indicates a good agreement, taking into account the differences in the mathematical models and the numerical methods. In the paper by [13], the comparison between the numerical and the experimental data is considered to be very satisfactory, even if the experimental error is not given. It is expected that the velocity distributions should be symmetric for u and antisymmetric for v. The dimensional drag force divided by the cylinder diameter over one period is also in good agreement with literature results (see Figure 6).Figure 4
Vorticity contours.
(a)
(b)
(c)
(d)Figure 5
Comparison with results taken of velocity profiles at four cross sections from [13] with constant x locations; phase angles: (a-b) 180°, (c-d) 210°, and (e-f) 330°.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Comparison of the drag force obtained using the present CFD code, Lescape, with [14] and experimental results taken from [13].


## 3.1.2. Validation of Force Calculations through an Oscillating Airfoil

An oscillating NACA0015 airfoil simultaneously experiencing pitching θ(t) and heaving h(t) motions is modelled. The pitching axis is located along the airfoil chord at position (xp,yp)=(1/3,0). The airfoil motion, described by [15], is defined as follows:(18)ht=H0sin⁡2πft+φ,θt=θ0sin⁡2πft.θ0 is the pitching amplitude, and the phase difference φ is set to 90°. The heaving velocity is then given by (19)Vyt=H02πfcos⁡2πft+φ.Based on the imposed motion and the upstream flow conditions, the airfoil experiences an effective angle of attackα(t) and an effective upstream velocity Veff(t) defined by(20)αt=arctan-VytU∞-θt,Vefft=U∞2+Vy2t,where the freestream velocity far upstream of the oscillating airfoil is U∞=1.To validate our simulations, a regime corresponding to the parametersRe=U∞c/ν=1100, H0/c=1, f=0.14, xp/c=1/3, and θ0=76.33° was computed. A view of the motion is sketched in Figure 7.Figure 7
Oscillating airfoil, sketch of the airfoil motion forH0/c=1, f=0.14, xp/c=1/3, and θ0=76.33°.The computational domain is[-14:16]×[-15:15]; the boundary conditions are set to Dirichlet on the left, top, and bottom, and a Neumann boundary condition is set on the right. The numerical results, obtained using dx=dy=5×10-3 and dt=1×10-3, are then compared to the force predictions presented by [15, 16].The force coefficientsCX, CY, and the moment coefficients Cm, are presented in Figure 8. The drag, lift, and moment coefficients are defined, respectively, as(21)CX=Fx1/2ρU∞2L,CY=Fy1/2ρU∞2L,Cm=T1/2ρU∞2L2,where the x-axis is aligned with U∞ and the y-axis is perpendicular to the far field velocity vector. Considering the Lescape’s drag coefficient (Figure 8(a)), the time evolution of the drag coefficient coincides with the results of Kinsey et al. and Campobasso et al. The numerical error for CX is below 3% for the maximum and minimum amplitudes. Lescape slightly overpredicts (by around 10%) CX during the decreasing phase of the drag. The evolution of the lift (Figure 8(b)) coincides with literature results. Lescape CY is in phase advance compared to Kinsey et al. and Campobasso et al., and the maximum time lag is below 5% just before the mid cycle T/2. The error in amplitude is globally below 5% during the cycle, with the maximum error occurring at the end of the cycle, and being about 10%. The pitching moment (Figure 8(c)) is the one with the least agreement with the literature results, with a maximal error of 20% at t/T=0.3 and t/T=0.8, the remaining error being below 10%. It must be noticed that, for pitching moments, Lescape shows a tendency to underestimate maximum (positive values) and overestimate minimum (negative values). The evolution in time of Cm follows the literature results; it increases when literature results increase and decreases in accordance with published results (meaning that, in the situation of fluid-solid interaction, the Lescape code will predict the rotation angle and its sign to an acceptable extent). For Cm, a small time lag at a maximum of 3% is observed. CX, CY, and Cm numerical results are satisfactory in the perspective of trajectory computations.Figure 8
Oscillating airfoil, comparison with the results obtained by [15, 16]. (a) Drag coefficient, (b) lift coefficient, and (c) pitching moment.
(a)
(b)
(c)


## 3.2. Plate Trajectories Database

For the parametric study, a sharp corner plate is used, as in Figure2. During these simulations, the rotation axis of the rectangular plate coincides with its barycenter. The effects of K, Fr, In, and Re on flat plate trajectories are studied in order to identify the important parameters at the start of the trajectory computation. A full factorial design with four factors and two levels [28] is used to build the test matrix.The simulations are performed usingdx=dy=1.×10-3 and a domain of [-3,8]×[-3,5], with the plate’s center initially located at (0,0), and with starting angles of ±π/3. For each factor, minimal and maximal values are selected around a base value: Re=7000, K=30, Fr=57.7, and In=8.44×10-2. Except for the Reynolds number, these are the typical values for shed ice pieces around an aircraft in a deicing situation [29]. For In, two values of the aspect ratio h/l are used, namely, 0.05 and 0.15. This is only a 1% variation, but this is a significant thickness change. Table 1 gives the 16 run factors. The 16 runs are executed for two plate starting angles and two release times. Thus, the database consists of 64 trajectories.Table 1
Parameters for numerical experiments.
Run number
Fr
Re
K
I
n
(
1
)
51.7
8750
37.5
8.52
×
1
0
-
2
(
2
)
51.7
8750
37.5
8.35
×
1
0
-
2
(
3
)
51.7
8750
22.5
8.52
×
1
0
-
2
(
4
)
51.7
8750
22.5
8.35
×
1
0
-
2
(
5
)
51.7
5250
37.5
8.52
×
1
0
-
2
(
6
)
51.7
5250
37.5
8.35
×
1
0
-
2
(
7
)
51.7
5250
22.5
8.52
×
1
0
-
2
(
8
)
51.7
5250
22.5
8.35
×
1
0
-
2
(
9
)
66.6
8750
37.5
8.52
×
1
0
-
2
(
10
)
66.6
8750
37.5
8.35
×
1
0
-
2
(
11
)
66.6
8750
22.5
8.52
×
1
0
-
2
(
12
)
66.6
8750
22.5
8.35
×
1
0
-
2
(
13
)
66.6
5250
37.5
8.52
×
1
0
-
2
(
14
)
66.6
5250
37.5
8.35
×
1
0
-
2
(
15
)
66.6
5250
22.5
8.52
×
1
0
-
2
(
16
)
66.6
5250
22.5
8.35
×
1
0
-
2The starting angle creates a separated flow around the plate. Vortices are shed periodically at low frequencies. The flow is first solved around a static plate, which goes on until the periodic flow is well established. The plate is released into the flow either at timet=20.0 s or at time t=21.1 s. The dimensionless position Y′ of the plate when it crosses the X′=1.125×10-3 line will be studied.The computed trajectories are illustrated in Figures9 and 10. To help visualize the individual trajectories, the results are separated into four. For a positive plate starting angle, these results are presented in Figure 9, with a further separation between the two shedding times. Similarly, the results for a negative plate starting angle are shown in Figure 10. The positions X and Y are plotted against each other. Since the initial angle of the plate is not null, vortex shedding is created behind the plate. The influence of vortex shedding is not negligible for a particular trajectory. For different starting times, the plate positions at the end of the trajectory calculations are different, as comparison between (a) and (b) shows. However, the results at different launching time will be grouped together to study the effects of other parameters. It can be seen that, in general, the plate accelerates faster in the X direction than in the Y direction. The gravity effects curve the trajectory. Here, the plate motion is mostly driven by aerodynamic forces and the starting angle has a dominant effect.Figure 9
Dimensional trajectories representative of the test matrix for a positive starting angle of the plate; Red Fr number51.7; Blue Fr number 66.6; square Re number 8750; circle Re number 5250; filled symbolsK37.5; empty symbolsK22.5; line inertia parameter 8.52×10-2; dashed line inertia parameter 8.35×10-2.
(a)
(b)Figure 10
Dimensional trajectories representative of the test matrix for a negative starting angle of the plate; Red Fr number51.7; Blue Fr number 66.6; square Re number 8750; circle Re number 5250; filled symbolsK37.5; empty symbolsK22.5; line inertia parameter 8.52×10-2; dashed line inertia parameter 8.35×10-2.
(a)
(b)Figure11 presents an example of the trajectories obtained for a plate with an initial positive angle. The dimensionless positions of the center of gravity for the 16 run parameters with an initial positive angle of attacks and a release at t=20.0 s is shown. The computed average trajectory of the 16 runs is also plotted, together with the upper and lower limits for a confidence interval. There is a 90% probability that the mean position of the plate gravity center at a given X′ position lies between the upper and lower limits. The confidence interval is computed using the sample variance and a Student probability distribution [30]. As can be seen, many trajectories are clearly outside the confidence interval. The mean trajectory gives a representative idea of the plate motion, but all the trajectory must be plotted to take into account the stochastic nature of the motion.Figure 11
Mean trajectory and individual trajectories for the16 runs with an initial positive angle of attack for the plate.The plate geometry, corresponding to run(1), is drawn at a constant time interval in Figure 12. For this run, the plate is rotating only slightly. The horizontal acceleration of the plate is visible as the spacing between its locations increases with the X′ position. The mean trajectory of the center of gravity is also plotted with the corresponding confidence interval. For run (1), the average trajectory is representative of the plate motion.Figure 12
Mean trajectory and plate orientation for run(1) for a plate with an initial negative angle of attack.The trajectory footprint probability maps used for the 3D trajectories analysis by [3] become histograms for 2D trajectories. At X′=1.125, a line is drawn in the Y direction and split into intervals. The heights of the histogram bar are computed by dividing the count of plate passage into an interval ΔY′ by 32, which represents the total number of runs [28]. The histograms associated with the two initial plate angular orientations are shown in Figure 13. In this figure, the Y′ position is on the horizontal axis, and two histograms are plotted, one for each plate orientation. The length of each bar gives the probability of passage of the plate gravity center in an interval of ΔY′=0.4×10-4. The initial plate orientation has a clear effect on the histograms, with the positive angle giving positive Y′ values, and the negative angle giving negative Y′ values. This is expected since an initial positive angle creates a positive lift coefficient, and a negative angle creates a negative lift coefficient. The two histograms are placed almost symmetrically around the Y′=0 line, indicating that gravity is not the dominant force of the problem at the beginning of the motion.Figure 13
Histogram of the probability of passage atY′ for X′=1.125×10-3.The initial plate orientation is the parameter that has the most effect on the trajectories. To study the effects of the other parameters, the database is split into two sets of32 runs: set one for the starting positive angle, and set two for the starting negative angle. Box plots are used to compare parameter effects [28]. Figure 14 shows the relation between the box plot and the probability histogram for the results of the flat plate with a starting negative angle. The central mark is the median of the 32 runs. The edges of the box are Q1, the 25th percentile, and Q3, the 75th percentile. The whiskers extend to Q1-1.5(Q3-Q1) and Q3+1.5(Q3-Q1), or the most extreme data points not considered outliers. The points for which Y′ are smaller than Q1-1.5(Q3-Q1) or greater than Q3+1.5(Q3-Q1) are considered outliers. One outlier is visible in Figure 14.Figure 14
The box plot and the corresponding histogram for the32 runs of the flat plate with a negative initial angle.For the trajectories computed with a starting angle ofπ/3 (positive orientation) (Figure 15), the box plots show that the Fr parameter has a significant effect on the median value. Also, the median value for K- is lower than the 25th percentile for the K+ distribution, which confirms that gravity forces are more important for lower K values. For the two other dimensionless parameters, the median does not change significantly between minimum and maximum values. For example, for Re+, the median is between the edges of the box for Re- distribution, while that for Re- is between the edges of the box for Re+ distribution. The result dispersion is also changed by Fr values and increases with increasing Fr values. The whiskers extend from 3×10-4 to 6.7×10-4 for Fr+ and from 4.3×10-4 to 6.5×10-4 for Fr-. According to the box plot, lower Fr values will make the plate fly higher and will reduce the dispersion. As Fr is only present in (12), the result dispersion is caused by the plate rotation. We should recall that Fr is the dimensionless length of the plate; Ln=1/Fr2. It means that longer dimensionless length plates fly more easily and in a more predictable fashion.Figure 15
Boxplot of theY′ position at X′=1.125×10-3 for a plate with a starting positive angle orientation.The box plots for a plate at a starting angle of-π/3 (negative orientation) are shown in Figure 16. The Fr number once again has the most significant effects on the plate location Y′ at X′=1.125×10-3. The standard deviation increases noticeably with an increase of the Fr number. Decreasing the K decreases the median Y′ position of the plate, as was the case for the plate with a starting positive angle. Decreasing Re and In has almost no effects on the final center of gravity position.Figure 16
Boxplot of theY′ position at X′=1.125×10-3 for a plate with a starting negative angle orientation.This database parametric study focuses on the start of the trajectory. To give an order of magnitude, in typical ice trajectory experiments [29], the ice piece travels a dimensionless horizontal distance of 4.×10-3, three times the horizontal distance in this database. The major limitation of this study is the 2D computational domain that precludes any comparisons against 3D experimental results. The computed 2D results should be equivalent to 3D results with a plate of infinite span and different from the results obtained with a plate of short span extent. Short span plates are typical of most engineering applications and experiments. Another limitation is the small number of runs used for the statistical analysis. Although 64 runs are used to build the database, the effect of the plate starting angle on the final position is strong, and runs must be separated into two databases. Trajectory depends on launch time and multiple launch times should be used to populate a larger database. Nevertheless, the database size is comparable to the database size used by [7] for parametric study.The numerical errors on aerodynamic coefficients observed during the oscillating airfoil validation could also have an impact on the database results. The aerodynamic moment coefficient oscillates typically between −0.6 and 0.6 during one trajectory calculation. The number of oscillation cycles is usually above three. Assuming the airfoil validation errors apply for the moving plate, the moment coefficient maximum values are underestimated by20% and the minimum values overestimated by 20%. For one cycle, the numerical errors have probably neutral effects. Also, it is worth mentioning that, for Lagrangian methods, often used for trajectory calculations, 10% of errors on static aerodynamic coefficients are usually assumed, together with an unknown uncertainty on the damping function used to model the dynamic caused by rotational velocity [5]. An additional limitation specific to aeronautic is due to the values of the Re, which are one level of magnitude lower than those for aerodynamic applications. This may not be a critical limitation as the static experimental lift, drag, and momentum correlation for a rectangular plate do not show any dependence on Re [31]. However, Lagrangian trajectory calculations use corrections to the lift, drag, and momentum that depend on Reynolds. Actual CFD results show some Re number effects.


## 4. Conclusions

A database of64 trajectories for a 2D plate released in a moderate Reynolds number flow has been built using CFD. Apart from the plate starting angles and launch times, four parameters influence the trajectories: the Reynolds number, the Froude number, the Tachikawa number, and the dimensionless mass moment of inertia parameter. Both the Froude and the Tachikawa numbers have significant effects on the trajectories. It appears that increasing the Froude number increases the dispersion of the trajectories in a vertical plan. Lowering the Tachikawa numbers decreases the median value of the dimensionless vertical position. Launch times change significantly the trajectory, as indicated by scatter plots of the dimensionless vertical position as a function of horizontal position, due to the unsteady separated flow around the static plate.The CFD tool used to build the database solves the penalized Navier-Stokes equations. The numerical method for the flow is based on a VIC scheme. The aerodynamic forces and moments are computed by integrating the penalty term and are used to compute the velocity from the solid motion equation. The velocity is imposed on the solid area of the flow, identified by a level set function. The forces and moments are validated against literature results for an inline oscillating cylinder and a flapping airfoil.The mean trajectory with a confidence interval gives a representative idea of the plate motion as long as the two starting angles are processed separately. The histograms of passage probabilities at a vertical position show the strong effect of the starting plate angle. The box plot representation allows a comparison of the effects of the four parameters on a single graph.In the future, it could be interesting to carry out the same study with a higher Reynolds number more representative of the aerodynamic flow in deicing situations. Also, the same kind of study could be conducted in 3D to see if the same conclusions stand. Progress in high performance computing is however sorely needed if large enough databases are to be built.



---

*Source: 9638430-2018-02-01.xml*