Large-eddy simulations of an autorotating square flat plate

Accepted Manuscript

Large-eddy simulations of an autorotating square flat plate Patricia X. Coronado Domenge, Carlos A. Velez, Tuhin Das PII: DOI: Reference:

S0307-904X(16)30054-3 10.1016/j.apm.2016.01.058 APM 11018

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

27 February 2015 11 January 2016 28 January 2016

Please cite this article as: Patricia X. Coronado Domenge, Carlos A. Velez, Tuhin Das, Large-eddy simulations of an autorotating square flat plate, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.01.058

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A model to simulate the complex FSI of an autorotating plate is proposed.

CR IP T

• LES and hybrid LES models ability to predict complex FSI found in autorotation.

• LES resolved the small flow structures that RANS is unable to compute.

• Hybrid models save computational time compared to pure LES providing good results.

AC

CE

PT

ED

M

AN US

• Computational results agree well with experimental data.

1

ACCEPTED MANUSCRIPT

Large-eddy simulations of an autorotating square flat plate

CR IP T

Patricia X. Coronado Domengea,1,∗, Carlos A. Veleza,1 , Tuhin Dasa,2 a Department of Mechanical and Aerospace Engineering University of Central Florida, Orlando, FL 32816, USA

Abstract

AN US

Large-eddy simulation (LES) turbulence models are underdeveloped in the area of fluid-structure interaction (FSI), specifically in autorotation applications. To gain a better understanding of FSI simulations, under the influence of strong turbulent interactions, several LES simulations were conducted to study the autorotation of a thin plate and compared to experimental measurements and RANS simulations found in literature. The plate is allowed to spin freely about

M

its center of mass and is located near the center of a wind tunnel with a free stream velocity of 5 m/s. Wall effects from the wind tunnel enclosure are ne-

ED

glected. In this work, a coupled Computational Fluid Dynamics (CFD) - Rigid Body Dynamics (RBD) model is proposed employing the delayed-detachededdy simulation (DDES) and the Smagorinsky turbulence models to resolve

PT

the subgrid-scale stress (SGS). The qualitative prediction of vortex structures and the qualitative computation of pressure coefficients are in good agreement with experimental results. When compared to RANS, the results from the LES

CE

models provide better predictions of the pressure coefficient. Moreover, LES accurately captures the transient behavior of the plate and close correspondence

AC

is found between the predicted and measured moment coefficient. Keywords: ∗ Corresponding

author Email addresses: [email protected] (Patricia X. Coronado Domenge), [email protected] (Carlos A. Velez), [email protected] (Tuhin Das) 1 PhD Candidate 2 Assistant Professor

Preprint submitted to Applied Mathematical Modelling

February 8, 2016

ACCEPTED MANUSCRIPT

autorotation, fluid-structure interaction, large-eddy simulation, delayed-detached-eddy simulation

CR IP T

1. Introduction

The term autorotation was first introduced by Riabouchinsky in 1935 [1] and was later defined as the continued rotation of an object lacking an external

power source due to a stream of air [2]. Accurate prediction of the fluid and rigid body dynamics of autorotation is known to be very challenging due to the

AN US

complex unsteady flow dynamics, involving fluid-structure interaction (FSI),

turbulent flow and flow separation, and most importantly the strong coupling of the two-way interaction between the fluid and solid. As a result, accurate unsteady and non-dissipative turbulence modeling is critical when resolving the aerodynamic non-linearity of autorotation. Many industrial applications, such as wing flutter and bridge oscillations, experience strong FSI effects which can

M

cause catastrophic failure, especially when materials susceptible to fatigue are involved [3]. Further research is required to develop current state of the art FSI methods, in specific the turbulent and rigid body motion interaction, before op-

applications.

ED

timization, coupled with CFD, can be used in the design of successive industrial

PT

At present, Reynolds average Navier-Stokes (RANS) equations based models are widely used for low run times when simulating turbulence, although they lack precision in the accurate prediction of flow separation given its time aver-

CE

aging of fluctuations and modeling of the Reynolds stress tensor. As a result, RANS models make it hard to preserve the vortex characteristics involved leav-

AC

ing a flow field void of incoherent structures. In contrast, large-eddy simulation (LES) models have proven to be a good alternative to RANS models as they preserve the characteristics of the vortex as shown by Eisenback and Friedrich [4], and most recently the Smagorinsky model with Van Driest damping has been implemented by Im et al. [5] achieving good agreement with experimental observations. LES has also been successful when used in particle tracking [6],

3

ACCEPTED MANUSCRIPT

compressible flow [7], combustion [8, 9, 10] and now in autorotation. Recently, LES has been successfully implemented in the study of turbomachinery as published by Li et al. [11], Tyacke et al. [12] and Watson et al.

CR IP T

[13]. However, as shown by Breuer et al. [14] and Feymark [15], there is still the need for further LES investigation in the area of FSI, especially in autorotation

without a prescribed set rotation. A great deal of the work in FSI is conducted using hybrid models, combining both RANS and LES methods, as presented by

Wang [16], Wang and Zha [17], and Shinde et al. [18]. The work on this paper

aims to provide new contributions to LES in the developing field of autorotation,

AN US

a subfield of FSI.

Not only are LES models promising, but they perform better than RANS models, and overcome several of the RANS model disadvantages. In LES, the governing equations are spatially filtered on the scale of the numerical grid. The large energy containing scales are resolved numerically, and the small scale eddies, which are generally more homogeneous and universal, are modeled. The

M

large eddies are strongly affected by the flow field geometry boundaries, therefore the direct computation of the large eddies by LES is more accurate than

ED

modeling the large eddies by RANS. This in turn helps to better simulate flow separation, an important factor in autorotation. Separated vortices create large pressure gradients on the surface of origin, which is the driving force of this

PT

autorotation, i.e., in laminar flow the autorotation is not sustained. The effect of the unresolved small scales of motion in LES is typically modeled by a

CE

subgrid-scale (SGS) model[19, 20, 21, 22, 23] or by the inherent dissipation in the numerical schemes[24, 25, 26, 27, 28, 29]. Because the statistics of the small scale turbulence are more isotropic and universal, a general physical model for

AC

small scale eddies is more plausible. For certain applications and complex flows that require solving for the wall

boundary layer, the CPU resource needed by LES is close to that of the Direct Numerical Simulation (DNS). As a result, pure LES might not be rigorously implemented for another 3 decades in engineering applications [30] and several hybrid RANS and LES models have been developed to overcome the intensive 4

ACCEPTED MANUSCRIPT

CPU requirements for LES, with runtimes between those of RANS and LES. One of these models was introduced by Spalart et al. in 1997 called detachededdy simulation (DES) [30] which divides a flow domain into a LES region far

CR IP T

away from a solid wall and a RANS region near a solid wall. Previous work for turbulence simulations for airfoils, cylinders and forbodies using DES have

shown encouraging results as seen in work done by Travin et al. [31], Spalart [32], Hansen and Forsythe [33], Viswanathan et al. [34], and Subbareddy and Candler [35].

The original DES model suffers from the downside that the transition from

AN US

RANS to LES may not be grid independent and as a result Spalart suggested a modification to his original model in 2006 [36], called delayed-detached-eddy simulation (DDES). In order for the transition to be independent from grid spacing, Spalart used a blending function to limit the DES length scale similar to the one used by Menter and Kuntz [37] for the Shear Stress Transport (SST) model. The DDES model has shown excellent agreement with experimental

M

data as well as a significant improvement from DES in work done by Wang and Zha [38], Coronado Domenge et al. [39] and Im et al. [5].

ED

For turbulence modeling, additional adjustments for the DDES model, improved DDES (iDDES), are then outlined by Travin et al. [40] concerning the definition of the LES length scale and the Wall-Modelled Large-Eddy Simula-

PT

tion (WMLES) helping to resolve the turbulence definition near solid walls. For a more comprehensive review of the LES and hybrid models mentioned in this

CE

paper the reader is directed to the review paper by Argyropoulos and Markatos [41].

In this paper, the flow around an autorotating flat plate is simulated, using

AC

several LES models, which was analyzed experimentally by Martinez-Vazquez et al. [42] and previously simulated using RANS models by Hargreaves et al. [43].

5

ACCEPTED MANUSCRIPT

2. Computational Method and Models 2.1. Computational Method 2.1.1. Smagorinsky

CR IP T

The LES Smagorinsky model performs spatial filtering of the velocity fluc-

tuations to decompose them into large scales, which are numerically resolved, and small scales, which are modeled. The Smagorinsky-Lilly model with Van

Driest damping [19] models the SGS by employing an eddy viscosity approach. In this approach it is hypothesized that a turbulent eddy viscosity exists at the small scales and that the stresses are in equilibrium at the interface between the

AN US

large and small scales. The eddy viscosity (˜ µt ) is defined as: µ ˜t = ρ¯Cs2 l2

p

2Sij Sij

(1)

where Cs is the Smagorinsky constant, Sij is the rate-of strain tensor, and l is the model length scale given by:

p

M

l = (∆)1/3

(1 − exp(−y + /26))3

(2)

ED

where ∆ is the volume of the cell.

Typical variations of the Smagorinsky model involve modifications to the SGS length scale expression (l) and the filter method used. Although in this

PT

approach the cube root of the cell volume is used to filter the scales, Gaussian filtering can be performed based on the strain rate or velocity fluctuations

CE

in flow. Additionally, other LES models seek to improve upon the definition of the model constant Cs . In the Smagorinsky-Lilly model this value is held constant and typically equal to 0.1. In the dynamic Smagorinsky method [22],

AC

later studied by [44] for complex geometries, the model constant is computed dynamically and allowed to change in time and space. In this work, a dynamic choric Smagorinsky model is employed with the following definition for the Smagorinsky model constant:

Cs =

6

(3)

ACCEPTED MANUSCRIPT

where the angled brackets represent averaging over the whole domain, valid in homogeneous turbulence, K is defined as:

and m as:

1 (ug ˜i u ˜j ), i uj − u 2

g 2 ). ˜ 2 − ||D|| m = ∆2 (4||D||

(4)

CR IP T

K=

(5)

˜ is modeled as the deviatoric component of the symmetric gradient of the D

AN US

˜ = dev(symm(∇U )). The dynamic Smagorinsky constant velocity field, i.e., D

ranges from 0 to 1. The constant grows near the boundary region towards unity and is zero in regions where velocity gradients are small [45]. This model assumes homogeneous turbulence at the SGS scales, although non-homogeneous models do exist in the literature [46].

M

2.1.2. DDES

The DDES formulation introduced by Spalart et al. [36], based on the Spalart-Allmaras (S-A) one equation turbulence model [47], suggests some mod-

ED

ifications to his previous DES model [30], given that in wide boundary layers and shallow separation regions the DES simulation can present erroneous behavior. This may occur when the thickness of the boundary layer is greater

PT

than the grid spacing parallel to the wall, making the transition from RANS to LES earlier. With the new modified DDES, the RANS model is retained longer

CE

for thick boundary layers independent of the grid spacing. The DES model is adjusted as follows. The SGS formulation from the S-A definition is modified by redefining the

AC

non-dimensional parameter r to:

rd =

νt + ν (Ui,j Ui,j )0.5 k 2 d2

(6)

where Ui,j are the velocity gradients, and the subscript d refers to delayed for DDES. This parameter is modified in this form so it can be applied to any 7

ACCEPTED MANUSCRIPT

eddy-viscous model. rd is applied in the following function:

(7)

CR IP T

fd = 1 − tanh([8rd ]3 )

The coefficients 8 and 3 are acquired from DDES flat plate boundary layer

tests[36] by matching the solution to the RANS values. Using Eq. 7, the DES distance to the nearest wall d˜ can be modified and be redefined for DDES as follows:

AN US

d˜ = d − fd max(0, d − CDES ∆)

(8)

where d˜ is filtered and ∆ is the largest spacing of the grid cell in all the directions. This modification in d˜ reduces the grey transition area between RANS and LES. The qualitative change of the new d˜ is very significant, depending now on the eddy-viscosity field. The DDES model can now refuse the transition to

M

LES if not ready, when the function fd , using the value of rd , indicates that the point still lies within the boundary layer. The opposite also occurs; when there

ED

is massive separation indicated by fd , the change from RANS to LES takes place in the simulation.

PT

2.2. Case description

2.2.1. Experimental details The computational model presented in this paper uses as its base model

CE

and for validations purposes the experiments conducted by Martinez-Vazquez et al. [42] in the wind tunnel at the University of Auckland. The plate was

AC

mounted on a turntable in an open test section with a range of wind speeds of Uw = 5, 7.5 and 10 m/s for the autorotation experiments. The 2.7 kg

square plate was made of polystyrene with 1 m length and 0.0254 m thickness, equipped with twenty-four pressure transducers arranged as seen in Fig. 1. The reader is referred to [42] for further details on the experimental setup. The data collected by the pressure transducers was used to compute the force coefficients

8

ACCEPTED MANUSCRIPT

given by CN = FN /(0.5ρUw2 A), where A = 1m2 . From this data set, drag, lift and moment coefficients were inferred as defined by CD = FNx /(0.5ρUw2 A), CL = FNy /(0.5ρUw2 A) and CM = T /(0.5ρUw2 l3 ), respectively, where D is the

CR IP T

drag force, L is the lift force, T is the torque and l is the characteristic length which in this case is 1 m. In this paper, the results for a wind speed of 5 m/s

are used for the model validation, as the experimental results of the pressure

M

AN US

coefficients for the 7.5 and 10 m/s cases are not readily available in [42].

ED

Figure 1: Sensor distribution on plate surface for Auckland experimental setup [42]

It is important to note that the experimental results show an increase in the dynamic pressure in the vicinity of the plate compared to the entrance of the

PT

tunnel and this is taken into consideration when calculating the aerodynamic force coefficients. The measured increase and employed in the CFD simulation

CE

is of 22% for a Uw = 5 m/s. The 22% increased in dynamic pressure measured by [42] has been employed in the CFD simulations, changing the free stream

AC

velocity in the vicinity of the plate from 5 m/s to 6.1 m/s.

2.2.2. CFD model Several LES models are used to study the autorotation of a square flat plate.

The mesh is completely structured and consists of 4.6 million hexahedral cells, with y + values between 10 and 100. Based on the ranges of y + , wall functions are employed for the turbulence viscosity (νt ) to model the shear and sub-layer 9

ACCEPTED MANUSCRIPT

profiles of the boundary layer. The computational domain is made up of two zones, a cylinder containing the plate and the outside rectangular outer domain, connected by an Arbitrary Mesh Interface or AMI [48] as seen in Fig. 2 and

CR IP T

Fig. 3. The outer rectangular test section has dimensions of 15c × 7c × 3.5c

where c = 1m. The plate’s center is located 1.2c from the bottom wall and 5c

PT

ED

M

AN US

from the inlet.

CE

Figure 2: Computational domain and boundary conditions setup for LES

As indicated by Spalart et al. [36], DDES is mesh independent (grid density)

and therefore a y + values between 10 and 100 should not be of concern. However,

AC

LES simulations, such as the Smagorinsky-Lilly model, need a better defined boundary layer and therefore preferably a y + < 1, without the use of a wall

function to model the boundary layer. Fig. 3(a) and 3(b) show the cross-section of the mesh around the flat plate for the two different y + values. The mesh with the lower y + consists of around 100,000 more cells. Fig. 4 shows the comparison of the lift and drag coefficients for the Smagorin10

CR IP T

ACCEPTED MANUSCRIPT

(b) y+ < 1

AN US

(a) y+ > 10

M

(c) Plate inside rotating cylinder

Figure 3: Fully structured computational domain cross-section of mesh around flat plate for

ED

(a) y+ > 10 and (b) y+ < 1 for inner rotating mesh

sky simulations with the low and high y + values. The plots show there is very

PT

little difference in the results and grid independence, therefore, because of computational time, the mesh with the larger y + will be used for the rest of the

CE

simulations and results presented in this paper. For the computation, the inlet boundary condition was set uniform in the

x-direction with an imposed turbulence of 5%. This value was chosen as it is

AC

typical for low turbulence wind tunnels. The turbulence is modeled as random fluctuations about a mean velocity (freestream) and it showed little effect on the end results given that the distance between the inlet and the plate was considerable (d = 5m). A Re = 3.34 × 105 for Uw = 5m/s is used. Note,

as previously mentioned, that the inlet velocity used in the experiment and

11

ACCEPTED MANUSCRIPT

smag y+ > 10 smag y+ < 1

4

smag y+ > 10 smag y+ < 1

6

4

Cl

Cd

2

2

0 −2 0

45

90

135

180

225

270

315

−2

360

0

45

90

CR IP T

0

135

180

225

270

315

360

Angle(◦ )

Angle(◦ )

Figure 4: Comparison of lift and drag coefficients for a full rotation for different y + values for LES

AN US

simulation is different than the dynamic pressure directly in front of the plate. In order to replicate the experiment as close as possible, the front, back and a single top section were set as fixed static pressure outlets at atmospheric pressure. All the walls have a no-slip boundary condition, except for the plate surfaces which are set as moving wall velocity to take into account the motion computed by the FSI of the rotating mesh so there is not a flux across the plate.

M

A coupled CFD and rigid body dynamics (RBD) model is carried out in OpenFOAM [49] using incompressible DDES, iDDES, Smagorinsky with Van-

ED

Driest damping and dynamic Smagorsinsky models to simulate the small scale turbulence. The RBD model implemented into OpenFOAM is adopted from the symplectic splitting method described by Dullweber et al. [50]. In every

PT

time step, this model calculates the plate angular velocity in response to the instantaneous fluid forces computed which in turns define the mesh motion,

CE

demonstrating the two way strong coupling between the solid body and fluid flow around it. The rigid body motion consists of a series of consecutive rotations to update the orientation matrices as well as the angular momentum. For

AC

the detailed equations, the reader is referred to [50]. Given that the plate is only free to rotate on the z-axis, the RBD model is reduced to a single degree of freedom. To provide a better match between the simulation and experimental condi-

tions, a bearing friction submodel was implemented as suggested by Hargreaves

12

ACCEPTED MANUSCRIPT

et al. [43]. This model adds a bearing friction torque, Tf ric , to the aerodynamic torques already acting on the plate during the simulation. Tf ric is defined as:

(9)

CR IP T

p Tf ric = (0.5µr d) (mg − L)2 + D2

where µr is the rolling friction coefficient of the bearings, and d is the bore

diameter of the bearing block. These were assumed from typical roller bearing

parameters to be 0.003 and 0.0254 m respectively, as these values where not given for the Auckland experiment. It was observed that the contribution to the total torque, Ttotal , from Tf ric was less than 1%.

AN US

The PISO [51] algorithm is used to solve the Navier-Stokes equation in time with the finite volume method. The 3-D incompressible version of OpenFOAM [26] is employed to solve the coupled pressure-velocity equations from the discretized momentum equation. A second order time marching Euler scheme is

used for the temporal derivatives, while a 3rd order cubic scheme is used for the

M

gradient, divergence and laplacian operations used in the finite volume calculations. No relaxation of the field variables or plate acceleration is used; instead, the Courant number is limited to 0.5, as suggested for LES with FSI problems

ED

to avoid any numerical instabilities [52]. All simulations were run in parallel on 102 Intel Xeon 64-bit processors with a time step of 0.0005 s. When running for a total of 120 hours, DDES ran for 5.9

PT

computational seconds, while Smagorinsky ran for 3.8 computational seconds. This shows the limitations of pure LES models compared to hybrid models as it

CE

takes approximately 50% more clock time to simulate the same computational time.

AC

3. Results and Discussion 3.1. Model Validation To validate the CFD-RBD model presented in this paper, a comparison is

shown between the results of several LES models, the Auckland experiment data [42] and the RANS simulations performed by Hargreaves et al. [43]. 13

ACCEPTED MANUSCRIPT

ddes iddes smag dynsmag

−5

Experiment Hargreaves2014 0

45

90

0

135 180 225 270 315 360

Experiment Hargreaves2014 0

45

90

Angle(◦ )

AN US Cp

5

−5

Experiment Hargreaves2014 45

90

135 180 225 270 315 360 Angle(◦ )

(c) Sensor 7

ED

−5

CE AC

45

90

90

135 180 225 270 315 360 Angle(◦ )

(d) Sensor 8

ddes iddes smag dynsmag

Cp

0

−5

Experiment Hargreaves2014

0

45

5

PT

0

Experiment Hargreaves2014

0

ddes iddes smag dynsmag

5

Experiment Hargreaves2014 0

135 180 225 270 315 360

ddes iddes smag dynsmag

0

M

Cp

0

0

Cp

(b) Sensor 6 ddes iddes smag dynsmag

−5

135 180 225 270 315 360 Angle(◦ )

(a) Sensor 1

5

CR IP T

0

−5

ddes iddes smag dynsmag

5

Cp

Cp

5

45

90

135 180 225 270 315 360 Angle(◦ )

Angle(◦ )

(e) Sensor 11

(f) Sensor 13

Figure 5: CFD and experimental pressure coefficients at various sensor locations

Fig. 5 shows the pressure coefficients for half a rotation at various sensor

locations. The plate rotates about its z-axis in a clockwise motion. The experimental data presented is the average over a large number of cycles for a 120 14

ACCEPTED MANUSCRIPT

sec time frame, while all the simulation data is from one single cycle. This is the case because of the large clock time needed to run CFD simulations. Running the full time frame in CFD would require approximately 2880 hours of

CR IP T

clock time. Qualitatively the computed results agree well with the experimental data. At sensors 6 and 11, located near the edge of the plate, the LES results are significantly better compared to the RANS results. The maximum deviation from simulation to experiment is reduced by about 50% at sensor 6 and sensor 11, and about 10% at sensor 7 and 8. The RANS results predict a much

larger pressure peak which is resolved by LES, which can be attributed to LES

AN US

being able to correctly represent the vortex core pressures as it has the ability to resolve a large portion of the small flow structures that RANS is unable to do.

Numerical results for the moment coefficient show close resemblance to experimental values as seen in Fig. 6. From Fig. 6 it can be seen, from the peaks at each rotation, that every other Cm is over-predicted compared to the exper-

M

imental values. These plots show that Cm consists of several harmonics as seen in Fig. 7. The frequency-amplitude signal shown in Fig 7 was generated using a

ED

discrete Fourier Transform as outlined in [43]. LES results show adequate prediction of the frequencies response, but over-predicts the dominant frequency by approximately 0.1 Hz. This dominant frequency is the same as the vortex

PT

shedding frequency which occurs twice in every rotational cycle. 0.3

0.2

0.2

CE

0.1

0

Cm

Cm

0 −0.1

−0.2

−0.2 ddes iddes Experiment

AC

−0.4

0

1

2

3

4

5

6

7

8

9

ddes iddes Experiment

−0.3 −0.4 0.2

10

0.4

0.6

Time (s)

0.8

1

1.2

1.4

1.6

1.8

2

Time (s)

Figure 6: Moment coefficients

Comparisons for the lift and drag coefficients between the simulations and

15

ACCEPTED MANUSCRIPT

0.12 ddes iddes smag Experiment

0.08 0.06 0.04 0.02 0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Frequency (Hz)

CR IP T

Amplitude

0.1

2.2

2.4

Figure 7: Frequency domain representation of computed moment coefficient

the experimental results can be seen in Fig. 8. The CFD results present an

AN US

over-estimation for both aerodynamic coefficients, however, this same issue was found in [53] where RANS simulations were performed, and in the same way the differences cannot currently be fully explained. Alternatively, mass eccentricity, a possible inadequate bearing friction model and inaccuracies in the CFD-RBD simulation due to numerical limitation where the plate experiences stalled conditions may further affect the prediction of lift and drag forces. Also,

M

it is important to note that this over-prediction can be attributed to the time averaging done for the experimental results, as it can be seen in Fig. 6 every

simulation. 4

ED

other peak presents lower values compared to the ones given by the numerical

iddes smag dynSmag

4 Cd

Cl

PT

2

0

−2

iddes smag dynSmag

6

2

0 Experiment

CE

Experiment

0

45

90

135

180

225

270

315

−2

360

0

45

90

135

180

225

270

315

360

Angle(◦ )

Angle(◦ )

AC

Figure 8: Lift and drag coefficients for a full rotation for DDES, iDDES and experimental data

Fig. 8 shows that variations in the LES model predictions are greater for the

lift coefficient, where as all three models predict almost identical behavior for the drag coefficient. This is expected as pressure forces are less dominant in the 16

ACCEPTED MANUSCRIPT

lift direction, when compared to drag. Consequently, the lift force is influenced more by the viscous forces and turbulent surface interactions which are highly dependent on the LES model in use.

CR IP T

In general, the CFD-RBD model presented in this paper presents results that very closely resemble those from the Auckland experiment. It is also seen that

the four different turbulence models used in this paper give relatively similar results. Therefore, it is concluded that hybrid methods such as DDES and

iDDES can be used in FSI autorotation instead of pure LES methods, saving

AN US

computational time. 3.2. Results

Fig. 9 shows the instantaneous static pressure contours at the back surface of the plate at different angles of rotation. The contour levels for the Smagorinsky (top three plots) and the iDDES (middle three plots) present the same contour levels as the RANS solutions (bottom three plots) presented by Hargreaves et al.

M

[43]. Although qualitatively the contours for all three results presented seem to be very similar, both Smagorinsky and iDDES results seem to show a slightly

ED

better resolution of the vortices near the corners and side edges of the plate than the RANS solution, especially at an angle of attack of 120o . In Fig. 11 and Fig. 10 the formation of the leading and trailing edge vortices

PT

in the wake of the autorotating plate can be seen. At the initial time (T = 0.5sec) the plate is near vertical, which creates strong tip vortices at both the leading and trailing edges making the plate unstable. After 1 second the plate

CE

rotates an additional 90◦ allowing for vortex shedding to occur and recirculate downstream. At the T = 1sec, the tip vortex, that has shed from the bottom of

AC

the plate, collides with the bottom surface of the plate. As the vortex collides with the plate large pressure gradients are rapidly formed and then dissipated. These vortex-structure interactions are undoubtedly chaotic and are the root of instability and inconsistent behavior in autorotation. Notice that even at the last time (T = 3.5sec) shed vortices have not traveled far downstream and still interact directly with the plates downstream pressure field. These vortex17

(b) α = 20o

(d) α = 10o

(e) α = 20o

(c) α = 120o

AN US

(a) α = 6o

CR IP T

ACCEPTED MANUSCRIPT

ED

M

(f) α = 120o

(g) α = 0o

(h) α = 30o

(i) α = 120o

PT

Figure 9: Instantaneous pressure contours on the rear face of the plate at (a) α = 6o , (b) α = 20o , and (c) α = 120o for Smagorinsky; at (d) α = 10o , (e) α = 20o , and (f) α = 120o for iDDES; and at (g) α = 0o , (h) α = 30o , and (i) α = 120o from Hargreaves et al. RANS

CE

solutions [43]

blade interactions [54] strongly couple turbulence closure models to the dynamic

AC

motion of the rigid body. This two-way interaction is the largest source of error in FSI, since the interacting turbulent eddies are never fully resolved.

18

ACCEPTED MANUSCRIPT

(b) T=1.0 sec, α = 186o

(c) T=1.5 sec, α = 294o

(d) T=2.0 sec, α = 372o

(e) T=2.5 sec, α = 484o

(f) T=3.0 sec, α = 558o

AN US

CR IP T

(a) T=0.5 sec, α = 110o

(g) T=3.5 sec, α = 667o

of rotation for Smagorinsky

ED

4. Conclusions

M

Figure 10: Instantaneous velocity contours on XY Plane at different times and different angles

A CFD-RBD model was proposed to simulate the complex fluid-structure interaction of an autorotating flat plate, subject to a 5 m/s free-stream air ve-

PT

locity. Two Smagorinsky LES models are compared to two Hybrid LES models to predict the autorotating plate motion. The computational results agree well

CE

with the experimental data from the Auckland experiment [42], especially the qualitative prediction of vortex structures, as well as the qualitative computation of pressure coefficients at the plate’s surface. This work also shows better

AC

agreement compared to RANS results previously found in literature [43], as LES is able to resolve the small flow structures that RANS is unable to compute. Hybrid models, such as DDES and iDDES, were shown to save computational time compared to pure LES and have been found to be an acceptable alternative in the simulation of a flat plate autorotation. The next step is to apply similar

19

ACCEPTED MANUSCRIPT

(b) T=1.0 sec, α = 191o

(c) T=1.5 sec, α = 300o

(d) T=2.0 sec, α = 388o

(e) T=2.5 sec, α = 502o

(f) T=3.0 sec, α = 554o

(g) T=3.5 sec, α = 664o

(h) T=4.0 sec, α = 742o

AN US

M

(i) T=4.5 sec, α = 858o

(k) T=5.5 sec, α = 1035o

ED

(j) T=5.0 sec, α = 919o

CR IP T

(a) T=0.5 sec, α = 113o

Figure 11: Instantaneous velocity contours on XY Plane at different times and different angles

PT

of rotation for iDDES

techniques to the autorotation of a helicopter rotor. The results of this work

CE

show great promise in the ability of LES and hybrid LES models to predict the complex FSI found in autorotation, which warrants further research.

AC

5. Acknowledgments The authors acknowledge the University of Central Florida Stokes Advanced

Research Computing Center for providing computational resources and support that have contributed to results reported herein. Carlos Velez would like to thank McKnight Fellowship by Florida Education Fund. Research at UCF was

20

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

supported by NSF Grant # 1250280.

21

ACCEPTED MANUSCRIPT

References [1] D. P. Riabouchinsky, Thirty years of theoretical and experimental research

282–348.

CR IP T

in fluid mechanics, Journal of the Royal Aeronautical Society 39 (1935)

[2] E. H. Smith, Autorotating wings: an experimental investigation, Journal of Fluid Mechanics 50 (3) (1971) 513–534.

[3] K. Billah, R. Scanlan, Resonacen, Tacoma Narrows Bridge Failure and

Undergraduate Physics Textbooks, American Journal of Physics 59 (2)

AN US

(1991) 118–124.

[4] S. Eisenback, R. Friedrich, Large-eddy simulation of flow separation on an airfoil at a high angle of attack and Re=105 using cartersian grids, Theoretical and Computational Fluid Dynamics 22 (2008) 213–225.

M

[5] H.-S. Im, G.-C. Zha, B. Dano, Large eddy simulation of coflow jet airfoil at high angle of attack, Journal of Fluids Engineering 136 (2) (2014) 1–10.

ED

[6] G. Duan, B. Chen, Large eddy simulation by particle method coupled with sub-particle-scale model and application to mixing layer flow, Applied Mathematical Modelling DOI: 10.1016/j.apm.2014.10.058.

PT

[7] D. A. Lysenko, I. S. Ertesv˚ ag, R. K. E., Large-eddy simulation of the flow over a circular cylinder at Reynolds number 3900 using the OpenFOAM

CE

toolbox, Flow, Turbulence and Combustion 89 (4) (2012) 491–518. [8] J. Janicka, A. Sadiki, Large eddy simulation of turbulent combustion sys-

AC

tems, Proceeding of the Combustion Institute 30 (1) (2005) 537–547.

[9] K. Li, L. Zhou, C. Chan, H. Wang, Large-eddy simulation of ethanol spray combustion using SOM combustion model and its experimental validation, Applied Mathematical Modelling 39 (1) (2015) 36–49.

22

ACCEPTED MANUSCRIPT

[10] C. Velez, S. Martin, A. Jemcov, S. Vasu, LES simulation of an enclosed turbulent reacting methane jet with the tabulated premixed CMC method, Journal of Engineering for Gas Turbines and Power (2015) GTP–15–1357.

CR IP T

[11] C. Li, S. Zhu, Y.-L. Xu, Y. Xiao, 2.5D large eddy simulation of vertical

axis wind turbine in consideration of high angle of attack flow, Renewable Energy 51 (2013) 317–330.

[12] J. Tyacke, P. Tucker, R. Jefferson-Loveday, N. R. Vadlamani, R. Watson, I. Naqavi, X. Yang, Large eddy simulation for turbines: Methodolo-

AN US

gies, Cost and Future Outlooks, Journal of Turbomachinery 136 (6) (2013) 061009.

[13] R. Watson, P. Tucker, Z. Wang, X. Yuan, Towards robust unstructured turbomachinery large eddy simulation, Computers and Fluids 118 (2015) 245–254.

M

[14] M. Breuer, G. D. Nayer, M. Munsch, T. Gallinger, R. Wuchner, Fluidstructure interaction using a partitioned semi-implicit predictor-corrector scheme for the application of large-eddy simulation, Journal of Fluids and

ED

Structures 29 (2012) 107–130. [15] A. Feymark, A large eddy simulation based fluid-structure interaction

PT

methodology with application in hydroelasticity, Ph.D. thesis, Chalmers University of Technology, 2013.

CE

[16] B. Wang, Detached-eddy simulation of flow non-linearity of fluid-structural interactions using high order schemes and parallel computation, Ph.D. the-

AC

sis, University of Miami, 2008.

[17] B.-Y. Wang, G.-C. Zha, Dettached eddy simulation of transonic airfoil limited cycle oscillation with high order WENO scheme, in: 47th AIAA Aerospace Sciences Meeting and Exhibit, 2009.

[18] V. Shinde, T. Marcel, Y. Hoarau, T. Deloze, G. Harrah, F. Baj, J. Cardolaccia, J. Magnaud, E. Longatte, M. Braza, Numerical simulation of the 23

ACCEPTED MANUSCRIPT

fluid-structure interaction in a tube array under cross flow at moderate and high Reynolds number, Journal of Fluids and Structures 47 (2014) 99–113.

tion, Monthly Weather Rev. 91 (1963) 99–164.

CR IP T

[19] J. S. Smagorinsky, General circulation experiments with the primitive equa-

[20] D. Lilly, The presentation of small-scale turbulent in numerical simulation experiments, IBM Scientific Computing Symp. on Environmental Sciences (1967) 195.

[21] J. Deardorff, A numerical study of three dimensional turbulent channel flow

AN US

at large reynolds numbers, Journal of Fluid Mechanics 41 (1970) 453–480.

[22] M. Germano, U. Piomelli, P. Moin, W. H. Cabot, A dynamic subgrid-scale eddy viscosity model, Physics of Fluids A 3 (1991) 1760–1765. [23] U. Piomeli, W. Cabot, P. Moin, S. Lee, Subgrid scale backscatter in tur-

M

bulent and transitional flows, Physics of Fluids A 3 (1991) 1766–1771. [24] T. Kawamura, K. Kuwahara, Computation of high reynolds number flow

1984.

ED

around a circular cylinder with surface roughness, AIAA Paper-84-0340,

[25] J. P. Boris, On large eddy simulation using subgrid turbulence models, in:

PT

J. Lumley (Ed.), Whither Turbulence? Turbulence at the Crossroads, vol. 357 of Lecture Notes in Physics, Springer-Verlag, 344–353, 1990.

CE

[26] J. Boris, F. Grinstein, E. Oran, R. Kolbe, New insights into large eddy simulation, Fluid Dynamics Research 10 (4-6) (1992) 199–228.

AC

[27] F. Grinstein, Dynamics of coherent structures and transition to turbulence in free square jet, AIAA Paper-96-0781, 1978.

[28] J. Boris, More for LES: A brief historical perspective of MILES, Implicit Large Eddy Simulation Cambridge University Press (2007) 9–38.

24

ACCEPTED MANUSCRIPT

[29] Y.-Q. Shen, G.-C. Zha, Comparison of high order schemes for large eddy simulation of circular cylinder flow, in: 47th AIAA Aerospace Sciences Meeting and Exhibit, 2009.

CR IP T

[30] P. Spalart, W.-H. Jou, M. Strelets, S. Allmaras, Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach, Advances in

DNS/LES, 1st AFOSR Int. Conf. on DNS/LES, Greyden Press, Columbus, H., 1997.

[31] A. Travin, M. Shur, M. Strelets, P. Spalart, Detached-eddy simulations past

AN US

a circular cylinder, Flow Turbulence and Combustion 63 (1999) 293–313.

[32] P. R. Spalart, Young-person’s guide to detached-eddy simulation grids, NASA/CR-2001-211032, 2001.

[33] R. P. Hansen, J. Forsythe, Large and detached eddy simulation of a circular cylinder using unstructured grids, AIAA-2003-0775, 2003.

M

[34] A. Viswanathan, K. Klismith, J. Forsythe, K. D. Squires, Detached-eddy simulation around a forebody at high angle of attack, AIAA-2003-0263,

ED

2003.

[35] P. Subbareddy, G. V. Candler, Numerical investigations of supersonic base flows using DES, AIAA Paper 2005-0886, 2005.

PT

[36] P. Spalart, S. Deck, M. Shur, K. Squires, A new version of detached-eddy simulation, resistant to ambiguous grid densities, Theoretical and Compu-

CE

tational Fluid Dynamics 20 (2006) 181–195.

[37] F. Menter, M. Kuntz, The Aerodynamics of Heavy Vehicles: Trucks, Buses

AC

and Trains, chap. Adaptation of eddy-viscosity turbulence models to unsteady separated flow behind vehicles, Lecture notes in applied and computational mechanics, Springer-Verlag, 339–352, 2004.

[38] B.-Y. Wang, G.-C. Zha, Detached eddy simulations of a circular cylinder using a low diffusion E-CUSP and high-order WENO scheme, in: AIAA 38th Fluid Dynamics Conference, 2008. 25

ACCEPTED MANUSCRIPT

[39] P. X. Coronado Domenge, B. Wang, G.-C. Zha, Delayed-detached-eddy simulation of shock wave/turbulent boundary layer interaction, in: 48th AIAA Aerospace Sciences Meeting, AIAA 2010-109, 2010.

CR IP T

[40] A. Travin, M. Shur, P. Spalart, M. Strelets, Improvement of delayed detached-eddy simulation for LES with wall modelling, in: European Conference on Computational Fluid Dynamics, 2006.

[41] C. Argyropoulos, N. Markatos, Recent advances on the numerical modelling of turbulent flows, Applied Mathematical Modelling 39 (2) (2015) 693–732.

AN US

[42] P. Martinez-Vazquez, C. Baker, M. Sterling, A. Quinn, P. Richards, The flight of wind borne debris: An experimental, analytical and numerical investigation. Part II (experimental work), in: The Seventh Asia-Pacific Conference on Wind Engineering, 2009.

[43] D. Hargreaves, B. Kakimpa, J. Owen, The computationl fluid dynamics

M

modelling of the autorotation of square, flat plates, Journal of Fluids and Structures 46 (2014) 111–133.

ED

[44] A. Scotti, C. Meneveau, M. Fatica, Dynamic Smagorinsky model on anisotropic grids, Physics of Fluids 9 (6) (1997) 1856–1858. [45] A. Tejada-Martinez, K. Jansen, A dynamic Smagorinsky model with dy-

PT

namic determination of the filter width ratio, Tech. Rep., Rensselaer Polytechnic Institute, 2004.

CE

[46] E. L´evˆeque, F. Toschi, L. Shao, J.-P. Bertoglio, Shear-improved Smagorinsky model for large-eddy simulation of wall-bounded turbulent flows, Jour-

AC

nal of Fluid Mechanics 570 (2007) 491–502.

[47] P. Spalart, S. Allmaras, A One-equation Turbulence Model for Aerodynamic Flows, AIAA-92-0439, 1992.

[48] P. Farrell, J. Maddison, Conservative interpolation between volume meshes by local Galerkin projection, Computer Methods in Applied Mechanics and Engineering 200 (1-4) (2011) 89–100. 26

ACCEPTED MANUSCRIPT

[49] H. G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Computers in Physics 12 (6) (1998) 620–631.

CR IP T

[50] A. Dullweber, B. Leimkuhler, R. McLachlan, Symplectic splitting methods for rigid body molecular dynamics, Journal of Chemical Physics 107 (15) (1997) 5840–5851.

[51] H. Jasak, A. Jemcov, A. Tukovic, OpenFOAM: A C++ library for complex

physics simulations, in: International Workshop on Coupled Methods in

AN US

Numerical Dynamics IUC, Dubrovnik, Croatia, 2007.

[52] C. Brebbia, G. Rodriguez (Eds.), Fluid Structure Interaction VII, WIT Press, 2013.

[53] B. Kakimpa, D. Hargreaves, J. Owen, P. Martinez-Vazquez, C. Baker, M. Sterling, A. Quinn, CFD modelling of free-flight and auto-rotation of

M

plate type debris, Wind and Structures 13 (2) (2010) 169–189. [54] P. Coronado, C. Velez, M. Ilie, Rotorcraft blade-vortex street interactions;

ED

critical aerodynamic aspects, in: 51st AIAA Aerospace Sciences Meeting,

AC

CE

PT

AIAA 2013–0805, 2013.

27