# Modeling fluid structure interaction over a flexible fin attached to a NACA0012 airfoil

TABLE OF CONTENTS ACKNOWLEDGMENTS ii LIST OF FIGURES vii CHAPTER I. INTRODUCTION 1 1.1 Motivation 1 1.2 Basic design 2 1.3 Fluid structure interactions 4 1.4 Numerical representation of FSI problem .. 9 1.5 Literature review of immersed boundary methods 11 1.6 Manuscript organization 13 II. THE NAVIER STOKES FINITE DIFFERENCE FLUID DYNAMIC SOLVER 17 2.1 Transforming Cartesian to curvilinear coordinates 18 2.2 Spatial discretization of continuity and momentum equations 27 2.3 Reynolds stress modeling 3 2 2.4 Temporal discretization of continuity and momentum equations 3 9 III. THE SUBDIVISION FINITE ELEMENT STRUCTURAL DYNAMIC SOLVER 41 3.1 Basic formulation 41 3.2 Finite element discretization 47 iv

Table of Contents—continued CHAPTER IV. IMMERSED BOUNDARY METHODS 52 4.1 Non-boundary conforming methods 52 4.2 Treatment of immersed boundary (boun dary reconstruction) 57 V. VALIDATION OF THE FLUID DYNAMIC AND STRUCTURAL DYNAMIC SOLVERS 66 5.1 Turbulent flow over a NACA0012 airfoil (body-fitted coordinates) 67 5.2 Unsteady deflection of a simply sup ported square plate 75 VI. VALIDATION OF THE FINITE DIFFERENCE IMMERSED BOUNDARY NAVIER-STOKES SOLVER 78 6.1 Flow over a circular cylinder 80 6.2 Laminar flow over a NACA0012 airfoil 85 6.3 Unsteady flow over a circular cylinder at Reynolds number of 200 90 6.4 Turbulent flow over a NACA0012 airfoil at Reynolds number of 170000 96 6.5 Laminar flow over an infinitesimally thin flat plate 101 6.6 Flow over a flapping flat plate using immersed boundary method 105 VII. FLUID STRUCTURE INTERACTION OVER A PASSIVE FLEXIBLE FLAT PLATE 121 7.1 Problem description 122 7.2 Computation domain 123 7.3 Coupling CFD and CSD solvers 126 7.4 Results and discussions 129 v

Table of Contents—continued CHAPTER VIII. MODELING FLUID STRUCTURE INTERACTION OVER A FLEXIBLE FIN ATTACHED TO THE UPPER SURFACE OF A NACA0012 AIRFOIL 138 8.1 Experimental set-up 13 9 8.2 Computational set-up 14 0 8.3 Domain and grid details 141 8.4 Results and discussion 145 IX. CONCLUSIONS AND FUTURE DIRECTIONS 157 BIBLIOGRAPHY 162 VI

LIST OF FIGURES 1.1. Flexible fin attached to a NACA0012 airfoil .... 3 1.2. One-way coupling vs. two-way coupling 5 1.3. Numerical representation of the FSI problem .... 10 3.1. Shell geometry in the reference (left) and deformed (right) configuration 42 4.1. Classification of non-boundary conforming me thods 53 4.2. Schematic interpretation of one-direction li near interpolation for solid boundary with finite thickness 59 4.3. Linear interpolation technique for forcing at one grid point 60 4.4. Two-sided one-direction interpolation tech nique 63 4.5. Two-direction one-sided interpolation tech nique 64 5.1. C-type computational grid with nodes of 200*95 68 5.2. Close up of the body-fitted C-grid with a NACA0012 airfoil 69 5.3. Grid independence study by comparing surface pressure distributions 71 5.4. Pressure contours for three different angles of attack . . 73 5.5. Lift coefficient curve for different angels of attack for a NACA0012 airfoil . 74 5.6. Drag coefficient computed at different angles of attack 75 5.7. Time variation of central deflection of a simply supported plate 7 7 vii

Table of Contents—continued 6.1. Circular cylinder domain for 181*5*181 mesh .... 82 6.2. Surface pressure distribution over a circular cylinder at Re=4 0 83 6.3. Pressure contours over a circular cylinder at Re=40 84 6.4. Velocity contours over a circular cylinder at Re=40 85 6.5. Airfoil domain for 281*3*281 mesh 87 6.6. Pressure distributions along the airfoil 88 6.7. Pressure contours over a NACA0012 airfoil us ing immersed boundary method at Re=500 89 6.8. Velocity contours over a NACA0012 airfoil us ing immersed boundary method at Re=500 90 6.9. Contours of pressure at three different times around a circular cylinder at Re=200 using immersed boundary technique (t=30, 40 and 50 sec) 92 6.10. Velocity contours at three different time steps of 30, 40 and 50 sec 94 6.11. Lift coefficient variation with time 94 6.12. Cartesian airfoil domain for 281*3*281 mesh .... 97 6.13. Pressure contours over a turbulent NACA0012 airfoil at Re=170,000 using immersed boundary technique 99 6.14. Velocity contours over a turbulent NACA0012 airfoil at Re=170,000 using immersed boundary technique 100 6.15. Grid independence test on a NACA0012 computed using immersed boundary method 100 6.16. 2D grid with the infinitesimally thin flat plate 102 viii

Table of Contents—continued 6.17. Grid independence test for a infinitesimally thin plate computed using immersed boundary- method 103 6.18. Y-vorticty contours after t=20 104 6.19. Kinematic model showing the object position at different times 107 6.20. Kinematic model showing the object position at different times 109 6.21. Computational grid employed by both codes 110 6.22. Lift histogram comparisons (Solid line- Immersed boundary method, Dotted line- FLUENT) 112 6.23. Surface pressure distributions on the plate at the maximum amplitude position 112 6.24. Surface Pressure distributions on the plate at the base level (zero amplitude position) .... 113 6.25. Stream lines with respect the plate position at maximum, baseline and minimum amplitude positions 114 6.26. Time varying x-velocity contours with respect the plate position at maximum, baseline and minimum amplitude positions 116 6.27. Convergence of LX norm of error for the veloc ity field for flow over a circular cylinder .... 119 6.28. Convergence of L„ norm of error for the veloc ity field for flow over a NACA0012 airfoil 120 7.1. Schematic representation of the Mylar flap ping flat plate hinged at the leading edge 122 7.2. Computation fluid dynamics grid with the plate at baseline position before the flow is impulsively started 124 ix

Table of Contents—continued 7.3. Unstructured triangular mesh used in the structural code for the flexible flapping plate hinged at x=0 at the starting position . . . 126 7.4. Non-iterative scheme over all time 127 7.5. Iterative schemes over all time 128 7.6. Non-iterative schemes over each time step 129 7.7. Surface pressure Variation at t=1.3 130 7.8. Tail amplitude variations with respect to time 131 7.9. Lift histogram of the flexible flat plate 132 7.10. a) Pressure, bending stress and y-vorticity magnitude after 0.1 sec 133 7.10. b) Pressure, bending stress and y-vorticity magnitude after 0.2 sec 134 7.10. c) Pressure, bending stress and y-vorticity magnitude after 0.3 sec 135 7.10. d) Pressure, bending stress and y-vorticity magnitude after 0.4 sec 136 7.10. e) Pressure, bending stress and y-vorticity magnitude after 0.5 sec 137 8.1. Experimental setup of the fin attached to a NACA0012 airfoil (Courtesy of Dr. Liu and Dr. Montefort) 140 8.2. Schematic representation of the Mylar flapping flat plate attached to the upper surface of NACA0012 airfoil 141 8.3. Body fitted C-Grid generated using Meshpiolt ... 142 8.4. Zoomed up body-fitted grid showing the flexi ble fin orientation before the start of simu lation 143 x

Table of Contents—continued 8.5. Unstructured triangular mesh for the flexible flapping plate hinged at x=0 at the start of the simulation 144 8.6. Running Pressure, viscous and total drag coefficients compared with the experiments 14 6 8.7. Lift coefficient history with the fin 147 8.8. Tail end displacements with respect to time .... 147 8.9. a-1 Pressure forces on the fin after 1.8 sec ... 148 8.9. a-2 Stream lines with U-velocity contours af ter 1.8 sec 149 8.9. b-1 Pressure forces on the fin after 1.87 sec .. 149 8.9. b-2 Stream lines with U-velocity contours af ter 1.87 sec 150 8.9. c-1 Pressure force on the fin after 1.96 sec ... 150 8.9. c-2 Stream lines with U-velocity contours af ter 1.96 sec 151 8.9. d-1 Pressure force on the fin after 2.06 sec ... 151 8.9. d-2 Stream lines with U-velocity contours af ter 2 . 06 sec 152 8.9. e-1 Pressure force on the fin after 2.06 sec ... 152 8.9. e-2 Stream lines with U-velocity contours af ter 2 .14 sec 153 8.10. Drag coefficient as a function of angle of attack for fin and baseline 154 8.11. Pressure drag as a function of angle of at tack for baseline and fin (computations) 155 8.12. Viscous drag as a function as a function of angle of attack for baseline and fin (compu tations) 156 xi

CHAPTER I INTRODUCTION 1.1 Motivation Natural flyers like birds and insects are speculated to utilize their wing flexibility, particularly thin flexible fins, for more efficient flight and effective flow control in different flight regimes. The flexible fins seem to play an important role in flapping flights where highly unsteady aerodynamics is nonlinearly coupled with the deforming wing. Bird flight inspired many early aviation pioneers like Lilienthal and Wright brothers who used flexible thin wings for flight control. In Lippisch's [1] first (and last) successful man-powered ornithopter test in 1929, the dramatic effect of quasi- flexible trailing edges (made up of bamboo pieces attached to the rigid wing near the tips) on improving flapping propulsion over that of a rigid flapping wing was observed. However, as remarkable advances were made in fixed-wing aircraft, the potential benefit of wing flexibility had been largely ignored, partially because flexibility has usually been considered a dangerous factor and the associated unsteady aerodynamics is too complicated to handle. Recently, the potential advantage 1

of wing flexibility has been re-discovered, and relevant research has been supported by NASA's Morphing Program (McGowan 2001) [2]. Separation flow control is of massive importance to the performance of air, land, sea vehicles, turbo machinery and diffusers. Generally, it is desired to postpone flow separation so that form drag is reduced, stall is delayed, lift is enhanced, and pressure recovery is improved. Therefore, considerable research efforts have been made over years for separation control (or stall control) by using various techniques like synthetic jets (Glezer & Amitay [3], Mittal et al. [4]), vortex generators (Gad-el-Hak [5]), passive and active blowing (Gad-el-Hak [5]), local suction (Atik et al. [6]), flapping wing (Jones et al. [7]), and oscillating camber (Munday & Jacob [8] ) . In this thesis the concept of separation control (or stall control) for a post-stall NACA0012 airfoil is done by using a flexible fin to passively manipulate the interactions of the organized vortical structures in the separation region. 1.2 Basic design The basic design is illustrated below in the figure 1.1. A thin flexible fin is attached to the upper surface of a NACA0012 airfoil. The oscillation of the membrane,

induced by the separated wake from the post-stall airfoil interacts passively with the flow field to alter the global aerodynamic properties of the NACA0012 airfoil. The oscillations allied with the shape deformations change the overall pressure distribution on the fin, which in turn affect the fin dynamics. Thus this mutual effect of inertial forces and elastic forces can be considered through fluid-structure interactions (FSI) . Hence in this thesis a computational fluid dynamic solver is combined with a computational structural dynamics solver in order to model these fluid structure interactions around the thin flexible fin attached to the upper surface of a NACA0012 airfoil that passively manipulates the flow field in fully separated flows. Figure 1.1: Flexible fin attached to a NACA0012 airfoil

1.3 Fluid structure interactions 4 There are two major numerical techniques to compute the solution of fluid structure interaction problems. They can be classified as monolithic methods and partition/segregated methods. In monolithic methods the complete system of non linear equations for the fluid and structure are coupled and integrated into one system and solved at their common interface [9, 10]. This procedure leads to a single matrix containing all equations and couplings [11] . This matrix might be large and ill- conditioned [11] and there is a chance for numerical difficulties in convergence of the solution. This could be a major problem when dealing with large geometries. On the other hand monolithic methods are considered to be more robust of the two numerical techniques. The second technique, segregated method, is the widely used method by commercial software packages where different software and different meshes are employed by the fluid and structural problems. In this method both the fluid as well as structural field are defined separately and solved and the interface conditions from the structure and the fluid are applied as boundary conditions at different times. These methods are very popular because the individual codes can be modified accordingly depending on the complexity of the problem. The

5 methodology presented in this thesis comes under the second category. Coupling the CFD and CSD solvers From a physical view point the fluid structure interaction problem is a combination of two problems, flow field and the structural field. These problems use different numerical procedures to compute the solution on different domains and meshes. In addition to using different meshes there is a need for data transfer (pressure from fluid to structure and displacement, velocity from structure to fluid) across the interface. The transfer of data which are the boundary conditions is a very important feature of fluid structure interactions. There are basically two different ways to transfer the data. Fluid Field (CFD solution) Pressure Force Structural Field (FEA solution) Fluid Field (CFD solution) Pressure Force Displacement & velocity- Structural Field (FEA solution) Figure 1.2: One-way coupling vs. two-way coupling

The easiest is one-way coupling. In this procedure the forces are transferred only one way from the fluid (CFD) solution to the structure (e.g. transferring static pressure loads from the fluid solution onto a structural model). The underlying assumption here is that the deformation of the solid is so small that it doesn't affect the overall fluid flow solution. The second coupling method is the widely used two-way coupling procedure for most problems involving large displacements. This is required when the fluid forces cause a significant oscillation of the structure. The results are mapped from the first solution to the second and back from the second to the first. Either the structural solution or the fluid solution takes the lead and at significant intervals the solution is mapped to and fro from fluid/structural to structural/fluid. For example when the total pressure force causes the structure to deflect, the solution is mapped from fluid to structure and in return the displacement (or position) and velocity of the structure is transferred back to the fluid domain. Immersed boundary method In the present thesis the term fluid structure interaction is considered as interaction of forces

(pressure) and the corresponding movement of fin (momentum interaction) rather than thermal interaction. Hence coupling the computational fluid dynamics (CFD) solver with the computational structural dynamics (CSD) solver provides an effective tool for calculating the fin dynamics. There are a number of factors that need to be handled in order to couple the CFD solver with the CSD solver. Other than handling the different spatial as well as temporal characteristics of each solver difficulty arises because of the moving interface (i.e. fin) present in the domain. These dynamically moving boundary problems are amongst the most demanding problems in contemporary computational fluid dynamics. The major complexity arises from the fact that generally all the fluid dynamic domains are described in Eulerian frame of reference. This method is suitable and works well when the boundary location doesn't change with respect to time. This becomes a problem when the boundary location changes with time. There are different techniques that have been proposed to account for this time dependent movement of the boundary, such as the overset grid method, dynamic meshing and coordinate transformations which can be applied to body conformal grids. In order to account for the motion of the boundary these grids need to be regenerated at every time step and also the old solution

needs to be projected onto the new grid. For problems involving large deformations and or large motions these grid regeneration methods are not only complex and time consuming but can have adverse effects on the simplicity, accuracy and efficiency of the solver. A study done by Liou & Pantula [13] where a dynamically moving flat plate was simulated using the commercial codes Fluent and Ansys CFX supports the above argument. Hence there is a need for developing a cost efficient numerical procedure that can deal with large boundary motions. An alternate to the boundary conforming methods which do not require the regeneration of grid at every time step is the widely used non conforming boundary fitted technique called immersed boundary method (IBM). In this method the boundary location need not be dependent on the mesh layout. The basic idea of the immersed boundary technique lies on the definition of the solid boundaries which may be static or dynamic. The immersed boundary technique mimics a solid body by means of suitably defined body forces applied to the discretized set of the momentum equations. These body forces impose a kinematic condition such that the velocity at each node point is coupled to the interpolated fluid velocity. The body force-field f is imposed so that a desired velocity distribution V can be specified on an immersed boundary [14] . This means

that we just add the body force f to the Navier-stokes equations and solve for u from the equation — + V(uu) = -S/P + vV(Vu) + f, V.W = 0 .X dt The main advantage of this approach is that f can be prescribed on a regular mesh so that the accuracy and efficiency of the solution procedure on simple grids are maintained. Another advantage of these formulations is the simplification of grid generation, especially in the case of moving boundaries where the need for regeneration or deformation of the grid is eliminated. 1.4 Numerical representation of FSI problem In figure 1.3 the numerical representation of procedure that is adopted and used in the present thesis is represented. The first step in the numerical procedure is to calculate the fluid solution on a eulerian grid. Then the total pressure force acting is exported as a condition for the structural field to obtain the velocity and the displacement. Then the location of the body is traced in a lagrangian fashion and appropriate virtual forces at interface locations are formulated and smoothly transferred onto the eulerian grid nodes using the immersed boundary technique. Then the solution is

10 advanced to next time step. The most convenient way for transferring the data between the CFD mesh and the FEA mesh would be if the nodes are concurrent. The pressure forces as well as the velocities and displacements need to be interpolated to the nearest node points if they are not coincident. CFD Mesh (Steady w.r.t time) Fluid Field (CFD solution) Export pressure forces from the nodes Structural Field (FEA solution) Immersed Boundary Method (Interpolate velocity and position onto the fixed CFD mesh) Export di spl acement & vel oci t y at each node poi nt FEA Mesh (Same node pos i t i ons as CFD mesh) Figure 1.3: Numerical represent at i on of the FSI problem The low Reynolds number K- s model of Launder and Sharma [15] i s al so incorporated i nt o the CFD code t o compute the t urbul ent scal es in the flow. One major fact or t hat needs to be addressed here when using the

immersed boundary technique is that the thickness of the fin is very small. The thickness of fin in the experiments is considered to be 0.0001 m but in the numerical procedure we consider the fin to be infinitesimally thin (i.e. thickness=0). Different interpolation techniques have been formulated to take care of this problem. These techniques are discussed in detail in the chapter IV of the thesis. 1.5 Literature review of immersed boundary methods Fluid structure interaction modeling has a wide range of industrial applications ranging from blood vessels, functioning heart value, elastic arteries to flexible tubes to parachutes to tents bridges to flapping flag to swimming fish. The first non-conformal boundary method to conduct fluid structure interactions was proposed by Peskin [16] where a fluid structure interaction over a cardiovascular circulation was studied assuming a low Reynolds number (Re) flow. In this calculation the boundary was modeled as a set of elements linked by springs. The body forces were easily computed using Hooke's law. Numerical difficulties arose when this method was applied to solid/rigid boundaries because the assumption of elements being elastic becomes unacceptable. Goldstein [17] applied it to solid

boundaries by introducing a feedback forcing approach that asymptotically enforces the desired boundary conditions on the solid. These come under the category of continuous forcing functions according to immersed boundary methods review by Mittal and Iaccarino [18]. The forcing is incorporated into continuous equations even before discretization. This feedback mechanism of Goldstein combined with the spectral method was used to simulate two-dimensional flow around a circular cylinder, as well as three-dimensional plane [17] and ribbed- turbulent channel flow [18] . These results were in good agreement with the reference data. The major drawback of this procedure is in order to calculate the feedback forcing two empirical constants related to flow frequencies were introduced. These two free constants need to be tuned according to the frequency of the flow. The equations became stiffer when the magnitude of these constants was high. This induced spurious oscillations and numerical instability which restricted the computational time step size. Saiki and Biringen [20] used the same forcing to compute the flow around fixed and rotating circular cylinders using a fourth-order central finite-difference scheme. Their results showed that the use of finite difference scheme eliminated the occurrence of spurious oscillations of flow at the

boundary. Even though this approach was successful at very low Reynolds number flows they cannot be used for high/moderate Reynolds number flows. Recently, Mohd-Yusof [21] and Fadlun [14] proposed a direct forcing embedded boundary formulation and showed that discrete time forcing is much more accurate and efficient compared to the feedback forcing. There are no empirical methods in this approach and the derivation of the forcing is explicit making the derivation of f flow independent. The methodology of Fadlun [14] was to introduce the forcing at the first grid point external to the immersed boundary using a velocity that, in a linear approximation, this point would have if the boundary had a desired velocity. Fadlun [14] reconstructed the solution at the fluid zones closest to the solid zone where as Kim [22], Majumdar [23] reconstructed the solution at ghost cells, which are solid zones closest to the fluid zone. Both the direct and feedback forcing procedures are discussed in detail in the later chapters. 1.6 Manuscript organization The second chapter of the manuscript/thesis begins by describing the Navier-Stokes fluid dynamic solver. The governing equations, spatial and temporal discretizations of the momentum and continuity, pressure velocity

couplings are discussed. Chapter three describes the subdivision finite element structural dynamics solver. The kinematics of deformation, finite element discretization and unstructured mesh generation techniques are presented in this chapter. Chapter four describes the immersed boundary method in detail. The application to moving boundaries as well as coupling between the CFD and the CSD solvers is discussed for the couple of FSI simulations conducted in this thesis. The velocity forcing applied on the fin is also discussed. In Chapter five the CFD solver and the CSD solver are validated. The turbulent flow over a NACA0 012 airfoil is simulated at Re=170,000 and compared with the published results. Different grid independence tests were done to validate the accuracy of the CFD solver. Then the CSD solver is validated by comparing unsteady simply supported beam solution with analytical results In Chapter six, different validations for the immersed boundary technique are presented. The major aim is to model the fluid structure interaction of a flexible fin attached to the upper surface of a post stall NACA0012 airfoil. The fin is considered infinitesimally thin. The Reynolds number of the flow is fixed at 63000.

15 In order for the immersed boundary technique to be applicable to the present design the following validations are imperative. • Steady state immersed boundary (laminar flow) a) Laminar flow over a NACA0012 airfoil at Reynolds number of 500. b) Laminar flow over a circular cylinder at Reynolds number of 40. • Unsteady immersed boundary (laminar flow) a) Unsteady laminar flow over a circular cylinder at Reynolds number of 200. • Unsteady flow over a infinitesimally thin flat plate. a) Unsteady flow over a flat plate at AOA of 30. • Steady state turbulent immersed boundary a) Turbulent flow over a NACA0012 airfoil on a rectangular grid using immersed boundary technique. • Moving boundary solution using immersed boundary method. a) Unsteady laminar flow around a flexible flat plate validated against a boundary fitted data of Pantula and Liou [13]. In Chapter seven combined fluid structure interaction modeling around a passively flapping flat plate at angle of attack of five is described.

In Chapter eight the fluid structure interaction around a flexible fin attached to a NACA0012 airfoil is described. Firstly the experimental set up is described followed by the numerical procedure used to model the fluid structure interactions around the fin. Finally the conclusions and future directions are presented in chapter nine.

17 CHAPTER II THE NAVIER STOKES FINITE DIFFERENCE FLUID DYNAMIC SOLVER In this chapter the details of the fluid dynamic part of the coupled fluid structure interaction solver is discussed. The basic solver in curvilinear coordinates is discussed in detail. The chapter starts with transforming the coordinate system from Cartesian to curvilinear coordinates followed by spatial and temporal discretizations. The numerical method is based on partial transformation approach, where the appropriate forms of the incompressible governing Navier-Stokes equations expressed in curvilinear coordinates, with velocity components expressed in Cartesian coordinates. The basic advantage of utilizing general curvilinear coordinates comes from the fact that the numerical fluxes can easily be estimated for non-orthogonal grids. The Navier-Stokes equations 2.1 to 2.4 are discretized in space, on a non- staggered mesh, using second order finite difference scheme for the pressure gradient and viscous terms, and second order upwind finite differencing for the convective terms. The upwind differencing of the convective terms eliminates the need for adding

18 artificial dissipation terms, to the right hand side of the momentum equations, to stabilize the numerical algorithm. This is due to the fact that a fixed amount of dissipation is inherent in the upwind differencing. The finite difference schemes employed for the pressure gradient, which pressure located on mesh points, and viscous terms, which calculating velocity components on half mesh points, are two-point central finite differencing, and for the convective terms are three- point one-sided finite differencing. The pressure- velocity equation is solved using the alternate- direction- implicit (ADI) approximate factorization method. A four-stage Runge-Kutta method is also used to advance the discrete equations in time. 2.1 Transforming Cartesian to curvilinear coordinates The Navier-Stokes equations in Cartesian coordinates can be written as dt dxyy x) dyyy yJ dzKH z) eux TT dux TT eux TT eux —- +ux —- +uv —- +uz —x - dt dx By dz d d d T H T H r ~\ xx ^ yx ^ z ox oy oz 2.1 dp dx +P8X