Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0001085/00001
## Material Information- Title:
- Membrane and adaptively-shaped wings for micro air vehicles
- Creator:
- Lian, Yongsheng (
*Author, Primary*) - Publication Date:
- 2003
- Copyright Date:
- 2003
## Subjects- Subjects / Keywords:
- Aerodynamics ( jstor )
Aircraft wings ( jstor ) Airfoil camber ( jstor ) Coordinate systems ( jstor ) Low reynolds number ( jstor ) Rigid wings ( jstor ) Shape optimization ( jstor ) Turbulence models ( jstor ) Velocity ( jstor ) Wing flaps ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Lian, Yongsheng. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 9/9/1999
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

MEMBRANE AND ADAPTIVELY-SHAPED WINGS FOR MICRO AIR VEHICLES By YONGSHENG LIAN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2003 Copyright 2003 by Yongsheng Lian To my parents and to my wife Lihui, for their love and support ACKNOWLEDGMENTS First, I would like to express my sincere gratitude to my advisor, Dr. Wei Shyy. He provided a good study environment and offered plenty of opportunities for me to explore and pursue my research interest. I would like to thank him for his enduring enthusiasm in educating me as a researcher and a person. I hope that I emulate his fine example of a successful international student, an intelligent and industrious scholar, a responsible and knowledgeable advisor, and an optimistic and humorous person. I owe thanks to several colleagues for their willingness to share their experiences and knowledge with me. I have benefited much from their collaboration in various areas: Dr. Peter G. Ifju, micro air vehicle; Dr. Andrew Kurdila, structural dynamics; and Dr. Raphael T. Haftka, optimization. Dr. Siddharth Thakur generously shared his time and experience, which made my work much easier. I enjoi, l. the good relationships with my professors. Dr. Renwei Mei gave me many helpful -,.-.-- -I'r ii and improved my knowledge of advanced fluid mechanics. Dr. Mark Sheplak spent much time gathering class materials for me. Both the University of Florida and the Mechanical and Aerospace Engineering Department provided a pleasant environment to study and live. I am proud to have been here; I hope that they will be proud to have me here. My family members in ('Ch1o i have given me tremendous support for my study abroad. Their trust and love have supported me over the years and will support me in the future. My wife, Lihui Bai, has been with me and supported me all the time. Her patience, understanding, and encouragement are the invaluable wealth of my life. But no acknowledgement could possibly state all that I owe to her. TABLE OF CONTENTS page ACKNOWLEDGMENTS ................... ...... iv ABSTRACT .. .. .. .. ... .. .. .. .. .. .. .. ... .. .. vii CHAPTER 1 INTRODUCTION ........................... 1 1.1 MAV and Membrane Wing Aerodynamics ........ .. .. 1 1.2 Motivations and (C! iII. 'ges ......... ........ ... 4 1.3 Outline ............ ...... ........ ..... 5 2 THEORETICAL AND COMPUTATIONAL MODELS FOR FLUID F L O W S . . . . . . . . . 8 2.1 Governing Equations and Fluid Solver ............... 9 2.1.1 Governing Equations for Fluid Flow ............. 9 2.1.2 Fluid Flow Solvers ........... ........... 12 2.1.3 Geometric Conservation Law ....... ......... 12 2.2 Turbulence Models . . ..... ....... 13 2.2.1 Reynolds Equations and k-E Turbulence Model ...... 13 2.2.2 Boundary Conditions for Turbulence Models . ... 15 2.3 Laminar-turbulent Transition ............. .. .. .. 18 3 MOVING GRID TECHNIQUE .................. ... .19 3.1 Perturbation Method .................. ..... .. 20 3.2 Master/Slave Concepts .................. ... .. 21 3.3 Numerical Tests .................. ......... .. 26 4 ACTIVE FLOW CONTROL WITH AN ADAPTIVELY-SHAPED WING 27 4.1 Introduction ............. . . ...... 27 4.2 Experimental and Computational Set-up . . ...... 28 4.3 Numerical Study .................. ...... .. .. 29 4.3.1 Grid Sensitivity Ain -i ........... .. ... 29 4.3.2 Assessment of the Turbulence Model . . 30 4.3.3 Effects of Gaps between Flaps ............ .30 4.3.4 Effects of Flapping Amplitude ............ .32 4.4 Conclusions ............... ........... .. 34 5 MEMBRANE STRUCTURAL SOLVER ........ .......... 36 5.1 Introduction ................... ....... 36 5.2 Structural Solver ................... .... 38 5.2.1 Mooney-Rivlin Model ................ .. .. 38 5.2.2 Principle of Virtual Work .................. .. 39 5.2.3 Triangular Membrane Element . . ..... 40 5.2.4 Green-Lagrange Strain Tensor ............... .. 42 5.2.5 Internal Forces .................. .... 43 5.2.6 External Forces .................. .... .. 45 5.2.7 Dynamic Equations for the Membrane . . .... 46 5.2.8 Solution of the Dynamic Equations . . ..... 47 5.3 Numerical Test of Membrane Model ................ .. 48 6 MEMBRANE WING AERODYNAMICS FOR MAVS . . .. 53 6.1 Introduction .................. ......... .. 53 6.2 Coupling between the Fluid and Structural Solvers . ... 55 6.3 MAV Wing Aerodynamics and Structural Response . ... 58 6.3.1 Computational Background ............. .. 59 6.3.2 Rigid Wing Aerodynamics ............. .. .. 64 6.3.3 Membrane Wing Dynamics ............. .. 76 7 WING SHAPE OPTIMIZATION .................. ..... 88 7.1 Introduction ..... .......... ......... .. 88 7.2 Optimization and Computational Framework . . ... 89 7.3 Results and Discussion ............... ...... 93 7.3.1 Sensitivity Analysis . . ....... .... 93 7.3.2 Performance of the Optimized Rigid Wing . ... 95 7.3.3 Flexible Wing Optimization ............... .. 97 7.4 Summary and Conclusions ................ 102 8 CONCLUSIONS AND FUTURE WORK ..... . . 104 APPENDIX A APPENDIX STRAIN-DISPLACEMENT MATRIX . . .... 107 REFERENCES ................... ... ........ 109 BIOGRAPHICAL SKETCH .................. ......... .. 118 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MEMBRANE AND ADAPTIVELY-SHAPED WINGS FOR MICRO AIR VEHICLES By Yongsheng Lian August 2003 C('! I: Wei Shyy Major Department: Mechanical and Aerospace Engineering Micro air vehicles (ilAVs), with wing span of 15 cm or less and flight speed around 10 m/s, have many applications in both civilian and military areas. The Reynolds number based on the given parameters is around 104, which often yields insufficient lift-to-drag ratio. Furthermore, one expects the unsteady effect to be noticeable for such flight vehicles. The flexible wing has been demonstrated to exhibit favorable characteristics such as passive adaptation to the flight environment and d. 1 .1 stall. The present study focuses on developing computational and modeling capabilities to better understand the MAV aerodynamics. Both flexible wings, utilizing mem- brane materials, and adaptively-shaped wings, utilizing piezo-actuated flaps, have been studied. In the adaptively-shaped wing study, we use piezo-actuated flaps to actively control the flow. We assess the impacts of the flap geometry, flapping am- plitude, and turbulence modeling on the flow structure with a parallel experimental effort. The membrane wing uses a passive control mechanism to delay the stall angle and to provide a smoother flight platform. Our study focuses on the mutual interac- tions between the membrane wing and its surrounding viscous flow. We compare the lift-to-drag ratio and the flow structure between the flexible wing and the correspond- ing rigid wing. We also investigate the aerodynamic characteristics associated with the low Reynolds number and low aspect ratio wing. To assist our study, we propose an automatic and efficient moving grid technique to facilitate the fluid and structure interaction computations; we also present a dynamic membrane model to study the intrinsic large deformation of the flexible membrane wing. Solutions obtained from the three-dimensional Navier-Stokes equations are presented to highlight the salient features of the wing aerodynamics. Besides the aerodynamic study, we also perform shape optimization to improve the membrane wing performance. Since direct optimization of a membrane wing is too time consuming to be practical, we optimize a surrogate rigid wing model based on an integrated optimization algorithm, which consists of a N ,',i i-Stokes solver, an automatic grid generation tool, and a gradient-based optimizer. Then, we assess the membrane wing performance based on the outcome from the surrogate model. Our numerical results confirm that the membrane wing exhibits consistent improvement in the lift-to-drag ratio with the surrogate model. CHAPTER 1 INTRODUCTION 1.1 MAV and Membrane Wing Aerodynamics Micro air vehicles ( \!AVs) with a maximal dimension of 15 cm or less and a flight speed around 10 m/s are of interest to both military and civilian applications. Equipped with a video camera or a sensor, these vehicles can perform surveillance, reconnaissance, targeting, and biochemical sensing at a remote or otherwise hazardous location [83]. With the rapid progress made in structural and material technologies, miniaturization of power plants, communication, visualization, and control devices, several groups have developed successful MAVs [24, 34]. As its size is reduced, MAV gains favorable scaling characteristics, including low inertia and reduced stalling speed [83]. However, in terms of aerodynamics, control, range, and maneuverability, many outstanding issues are still present. First, an MAV flies in a low Reynolds number regime (104r105) due to its low flight speed and small dimension. Such a flight environment is often accompanied by laminar boundary 1.-r separation, transition, and low lift-to-drag ratio. Second, an MAV wing typically has a low aspect ratio, which causes strong vortical flow structures and increases the induced drag. Third, an MAV is susceptible to rolling instabilities, which becomes even more serious because of the existence of tip vortices. Fourth, an MAV is sensitive to wing fluctuations, which can be comparable to MAV's flight speed and makes both the instantaneous flight Reynolds number and angle of attack vary substantially [83]. In the development of MAVs, two approaches exist in the MAV wing design: flapping wings and fixed wings. In the flapping wing area, pioneering work was done by Weis-Fogh [112], and Lighthill [57]. Recent experimental and numerical studies in this aspect were given by ('!ii i-i,, i and C'!i ,I.:i varthy [8], DeLaurier [11], Smith [90], Dickenson et al. [13], Ellington [17], Jones and Platzer [43, 44, 45], Katz [50], Kamakoti et al. [48], Liu and Kawachi [59], Vest and Katz [106, 107], Walker [109], and Wang [110]. Shyy et al. [83] gave an extensive review of flapping and fixed wing characteristics. Templin [96] presented the flapping wing spectrum of flight animals. Our study focuses on a fixed wing. Specifically, we will study a flexible membrane wing for MAV applications. Figure 1-1 shows a 15 cm MAV with a flexible wing designed by Ifju and coworkers [34]. The MAV has a flight speed between 20 and 40 km/h. Powered by an electric motor, it can fly for about 15 minutes. With a build-in video camera and a transmitter, the vehicle has a total mass around 40 g. The MAV has a rigid fuselage and leading edge spar which are made of laminate material. On each side of the wing there are three unidirectional carbon fiber battens. The spar and battens are covered with a latex membrane material. Overall, the wing has a span of 15.2 cm, a root chord of 13.2 cm, a mean chord of 10.5 cm, and a wing area of 160 cm2, which result in a low aspect ratio of 1.4. The wing has a variable camber decreasing from 7.5'. at the root to 2' at the tip. In this configuration, the maximum camber is located at about '7'. chord at the root. SA Figure 1-1: MAV with membrane wing. Flexible wings have different aerodynamic performance from rigid wings. One advantage of flexible membrane wings is that they facilitate passive shape adaptation and result in a d. 1,. .1 stall [83, 111]. Fig. 1-2, adopted from Waszak et al. [111], compares the lift curves versus angles of attack for rigid and membrane wings with configurations similar to those shown in Figure 1-1. Under modest angles of attack, the rigid and membrane wings demonstrate a similar lift curve slope of 2.9, with the rigid wing having a slightly higher lift coefficient. However, the membrane wings have a much higher stall angles of attack than the rigid wing, which is a key element in enhancing the stability and agility of MAVs. The rigid stalls at 24" while the flexible wings have stall angles between 30" and 49", which are similar to that of much lower aspect ratio rigid wings (AR=0.5 to 1.0). But, the very low aspect wings exhibit lower lift curve slopes of 1.3 to 1.7 instead of 2.9. Hence, flexible wings appear to combine the desirable performance of rigid wings with higher and lower aspect ratios by exhibiting stall behavior similar to that of rigid wings with aspect ratio of 0.5 to 1.0 and the lift generating capability similar to rigid wings with aspect ratio of 2.0 [111]. Another advantage of flexible membrane wings is that they can adapt to the wind gust and provide smoother flight platforms, which are verified by wind tunnel experiments [84]. 0 -- Rigid (Graphite) 3.0 -- --- 6-Batten (Monofilm) -- -- 2-Batten (Latex) - 2.2 2'I ++ C 1.5 L 0.75 - 0 .0 L -10 0 10 20 30 40 50 Angle of Attack (deg) Figure 1-2: Lift coefficient versus angle of attack. (From Waszak et al. [111]) 1.2 Motivations and ('!C iII. iges Though experiments have provided information on the membrane wing aero- dynamics, to fully understand the membrane aerodynamics and to appreciate the membrane wing's passive control mechanism, we need to perform a detailed numer- ical study. A numerical study of the membrane wing brings two challenges: the low Reynolds number and the low aspect ratio condition, and the interaction between the membrane wing and its surrounding viscous flows. First, the low Reynolds number condition presents tremendous challenges on MAVs study. In our study, the chord Reynolds number of the MAV is around 9x 104. In the range of Reynolds number of 104W 106, complex flow phenomena often occur. The laminar boundary 1lv- r separation, transition, and reattachment usually coexist on the upper wing surface. An accurate predication of the transition from laminar to turbulent state is important to evaluate the wing performance. However, a transi- tion model for general application is not yet available. Furthermore, low aspect ratio wing is usually accompanied by complicated vortical flow structures. The vortical flow structures introduce unsteady phenomena. Valuable information regarding low Reynolds number and low aspect ratio wings was reviewed by Lissaman [61] and Tani [95], as well as addressed in several proceedings [65, 66]. Gad-el-Hak [20] discussed approaches to prevent/alleviate flow separation on MAVs with microelectronicsme- chanical systems (\!I \!S). Lian and Shyy [53] performed a detailed numerical study of the flow separation on the MAV wing and its impacts on aerodynamic performance. Second, a membrane wing exhibits self-initiated vibrations even under steady state freestream condition. Waszak et al. [111] conducted wind tunnel experiments and found that the membrane wings exhibited vibration at 0(102) Hz in a steady freestream. Similar observation was reported numerically by Lian and Shyy [53]. Such vibrations and the associated shape deformation change the pressure distribution on the membrane wing, which, in turn, affects the membrane dynamics. This recipro- cal action leads to a fluid and structure interaction problem. Aeroelastic analysis method coupling a computational fluid dynamics (CFD) solver with a computational structural dynamics (CSD) solver provides an effective tool for the membrane wing study. Although both CFD and CSD have achieved great success individually, a cou- pled computation is still a challenging task [80]. To date, the study of a membrane wing and fluid flow interactions is limited. Jackson and C('! i-I1 ,, [38] adopted a three- dimensional potential flow-based solver and a membrane wing model to analyze the aeroelastic behavior of marine sails. Smith and Shyy [92] and Shyy et al. [88] pre- sented a computational approach to model the interaction between a two-dimensional flexible membrane wing and surrounding viscous flows. A combined numerical and experimental ,a i1, ~i-; of a flow over a two-dimensional flexible sail was conducted by Lorillu and Hureau [62]. In the coupled fluid and membrane wing interaction study, we need a fluid solver to accurately capture the transient phenomenon of the fluid flow surrounding the membrane wing; we need a structural solver to predict the in- herent nonlinear behaviors of the membrane material. To couple these two solvers, we also need to handle different time characteristics associated with each solver. Be- sides these, we need a moving grid technique to facilitate the fluid and structure interaction problem. An interface technique is also required to exchange information between these two solvers. 1.3 Outline This work study both the membrane wing and adaptively-shaped wing. First, we will develop a computational capability to study the fluid and structure interaction. Second, based on this capability, we will study the piezo-actuated flapping wing the membrane wing. Third, with the gained knowledge, we will optimize the membrane wing performance with modern optimization technique. The outline of the current work is as the following: In ('!C ipter 2 we present the fluid flow solver and physical models of the fluid dynamics. The Navier-Stokes equations for incompressible flow on curvilinear coordi- nates are employ, l along with a two-equation k-e turbulence model. In this chapter, the fluid solver is considered under the moving boundary framework. A brief in- troduction is given to the geometry conservation law [100] that p1 .,i- an important role in moving boundary problems. Even though the present work does not employ laminar-turbulent transition model, a brief review is offered. In ('!i lpter 3 we devote to moving grid technology, which is used by both the coupled fluid-structure study and the optimization. A robust, efficient and effective moving grid technique is proposed, which uses transfinite interpolation-based pertur- bation method in single block grid [56, 77], and employs the master/slave concept [29] for a multiblock grid. In ('! Ilpter 4 we investigate the active flow control with adaptively-shaped flaps for a low-Reynolds number wing. Results of the interaction between piezo-actuated flaps and surrounding flow are presented. In C'!i lpter 5 we focus on the structural solver, more specifically, a nonlinear dynamic membrane model. Detailed derivation of the membrane model, with an emphasis on the nonlinearity as well as the time integration will be presented. The membrane model will be assessed using well defined test problems. In C'!i lpter 6 we propose a coupled fluid and structure interaction algorithm. The coupling is realized through an interface technique and time synchronization. The numerical implementation of geometry conservation law is given and tested. The low Reynolds and low aspect ratio wing aerodynamics, including the laminar boundary l -.-r separation, tip vortices, and the unsteady phenomenon are investigated. Com- putational results will be presented to elucidate the membrane wing aerodynamics. 7 In ('!i ipter 7 we optimize of a flexible membrane wing with a gradient-based method. Since a direct optimization of a membrane wing requires substantial com- putational resource and time, a rigid wing is used as a surrogate. The membrane wing performance is then evaluated based on the surrogate model. In brief, the present work is motivated by, first, the practical need to understand the aerodynamics of membrane and adaptively-shaped wings, second, our interests in micro air vehicles. Original work include a moving grid algorithm for multi-block grid, a dynamic membrane model, numerical study of interaction between fluid and a flexible wing, and shape optimization of a flexible wing. CHAPTER 2 THEORETICAL AND COMPUTATIONAL MODELS FOR FLUID FLOWS Our study investigates the aerodynamics of a membrane wing as well as adap- tively shaped wings under low Reynolds number conditions. During flight, the mem- brane wing undergoes shape change under the external forces; meanwhile, its shape variation affects the pattern and structure of its surrounding fluid flow. The mem- brane wing vibration was observed both experimentally [111] and numerically [53, 55]. In the adaptively shaped wing study, the piezo-actuated flaps vibrated with different frequencies and amplitudes. The flap vibrations cause unsteady phenomena like vor- tex shedding [30, 56]. It is desirable to have a computational capability to accurately capture the transient behavior of fluid flow and structural dynamics. Even thought steady flow computations have achieved significant success in their development and applications, little attention has been given to the systematic analy- sis of unsteady computations [79, 80]. Several factors are responsible for this situation. First, steady CFD is still the mainstream aerodynamic analysis, and the techniques used to enhance the performance of steady CFD computations may not be suitable to unsteady computations. For example, in steady computations, a variety of tech- niques have been proposed to eliminate the low frequency transients to accelerate the convergence to steady state; however, low frequency phenomena are intrinsic to many fluid and structure problems, and such techniques will be inappropriate to the unsteady analysis [80]. Second, the applicability of the turbulence model to unsteady flows is not well understood [79]. However, numerous approaches have been proposed for unsteady flow computa- tions. Pulliam [74] incorporated subiteration to improve time accuracy of conven- tional implicit schemes. Jameson [39] and Rumsey et al. [79] developed a dual time stepping technique in the context of multigrid methodology to enhance the efficiency and accuracy of traditional factored schemes. Issa [35], and Oliveira and Issa [71] pro- posed the Pressure-Implicit with Splitting of Operators (PISO) method to improve the convergence. For relevant information of these and some of the most frequently adopted approaches, see Shyy and Mittal [85]. The SIMPLE method [72, 82] and the PISO method [35, 71, 98] are two popular methods to solve the three-dimensional Navier-Stokes equations for incompressible flows. Both methods belong to the pressure-based approach [82] which devises an artificially derived pressure (correction) equation by manipulating the mass continuity and momentum equations. Three-dimensional N ,i.-, r-Stokes equations for incompressible flow in curvilinear coordinates are adopted in our study. Two methods are emplo-, '1 to solve the gov- erning equations, one is the SIMPLEC method [102] that is belong to the SIMPLE family, and the other is the PISO method [35]. In this chapter we briefly introduce these fluid solvers under the moving boundary framework. For detailed descriptions of these methods, we refer to Refs. [35, 72, 82, 98, 102]. To solve the moving bound- ary problems, we also emphasize the importance of the boundary conditions and the geometric conservation law [100]. An introduction of turbulence models, with an emphasis on the k-e turbulence model is presented. Though an accurate transition model for general application is beyond the scope of our study, we offer a brief review of the transition models at the end of this chapter for completeness purpose. 2.1 Governing Equations and Fluid Solver 2.1.1 Governing Equations for Fluid Flow The three-dimensional N ,1-i. r-Stokes equations for incompressible fluid flows in curvilinear coordinates [82] in strong conservative form can be written as the follow- ing, aU aV aw + + (2.1) ^+ a9+^ (2.1) a(Jpu) a (pUu) a (pVu) a (pWu) t+ + a + a(Jpv) 9 (pUv) 9 (pVv) a (pWv) at + 9+ + 9 a(Jpw) a(pUw) a(pVw) a(pWw) at + 9 + a + 9( a [ (j11iu + qi1- + q13u()] + P [(q21 + q__.. : + 923U()] + [a (qU31+ 32U+ "+ ) - [(fiP) + (f21P) + (f31p) + P (q21V + q._ + q23Vc)] + P, (31q +(I + q33V()] - [(f12P)+ (f22P) )+ (f32P) + G2a(, ).J, (2.2) P (qdjl(w + l-12 + q13w()] + [ (q21W+q_, + q23UW)] OTl J + P[(q3iw + 4 _: ,, + q33') [(fl1P) + (f2P)+ (f32P)] + G3(~ )J, (2.4) + [( + 41 _,+ Y1: )1 + [ (f21p + q(f'p) + i(fq)2 + G3((, ).J, (2.4) where u, v, and w are the velocity components in the x-, y-, and z-directions, re- spectively, p is the viscosity, G1, G2, and G3 are the body force components per unit volume in the x-, y-, and z-directions. The coefficients appearing in the equations are expressed as follows: f2 + f2 + f 12 211 122 13, 912 f221 + f222 + f223, 13 f321 + f322 + f323, 923 fll21 + f12f22 + f13f23, fll31 + f12f32 + f13f33, f31f21 + f32f22 + f13f23. The expressions for fj, i 1, 3, j=1, 3, can be found in [82]. The Jacobian transfor- mational matrix is defined as follows: (2.6) The determinant of the Jacobian matrix is expressed as J = x0yxz + xyzx, + xyjzz xyZzq xy,'z xyzq,. (2.7) In the governing equations, U is the flux through the control surface whose normal direction is ; V and W are the fluxes through the Tr and ( control surfaces. When the governing equations are considered under a moving grid framework, the grid velocities should be considered in the flux computations. In curvilinear coordinates, they are expressed as the following: U = fll(u V = f21(U W = f31(U i) + fl2(v I) + f22(V -) + f32(V ) + f13(W ) + f23( S+ f33(w (2.8) where, for example, x is the grid velocity component in the x-direction that can be approximated by the first order Euler method, x Xo At (2.9) where At is the time step and the superscript 0 refers to the previous time level. Similarly, the values of the other two components i and z can be estimated. (2.5) 9(x, y,z) J 9= ,C 2.1.2 Fluid Flow Solvers The SIMPLEC method [102] is a variant of the SIMPLE method originally pro- posed by Patankar [72] in Cartesian coordinates, and its extension to curvilinear coordinates was discussed by Shyy [82] to accommodate the complex geometry com- putation. The SIMPLEC method uses coordinated under-relaxations for the momen- tum and pressure corrections to improve the convergence that is intrinsically slow in the original SIMPLE method. The SIMPLEC method was also used to solve mov- ing boundary problems [53, 88]. Contrary to the SIMPLE family method, the PISO method is a non-iterative approach to handle the pressure-velocity coupling which splits the process of solution into a series of predictor and corrector steps. The PISO method is generally more efficient for transient flow computations [35, 98]. To solve the Navier-Stokes equations, we require proper boundary conditions for the computations. At the interface between the fluid and the solid structure, the following boundary conditions are applied, Xf = x, on F, (2.10) Oxf ax8 t at on F. (2.11) 9t 9t ' Eq. (2.10) expresses the compatibility between the wet fluid surface Xf and the struc- ture surface x, at the interface F; Eq. (2.11) requires that the fluid velocity at the interface F matches the surface velocity of the structure. These boundary conditions are provided by computing the displacement of the solid structure. 2.1.3 Geometric Conservation Law When body-fitted curvilinear coordinates are used in the computation, a trans- formation maps a physical flow region (x, y, z) onto a uniform computational space ( T1, () where conservation laws are carried out. To facilitate the coordinate trans- formation, the Jacobian matrix J is introduced. The Jacobian matrix is defined as in Eq. (2.6). The determinant of J represents a volume element in the transformed coordinates. In a moving grid problem, the computational grid needs to be updated with time; meanwhile, the Jacobian J needs to be updated simultaneously. Specific procedures are required to compute the effective value of J; otherwise, errors arise due to the inconsistency in evaluating geometric quantities and the adding of non-physical quantities. As pointed by Thomas and Lombard [100], in the updating process, the following geometric conservation law (GCL) needed to be satisfied, d Jdidqd(= (V W,)ddqd(, (2.12) where V is the volume bounded by a close surface S, W8 denotes the local velocity of the boundary surface S. Thomas and Lombard [100] proposed to evaluate J while maintaining the geometric conservation law by setting p=1, V=0 from the continuity equation, i.e., Jt + (JW) + (Jt)k, + (J(c = 0. (2.13) Implications of the geometric conservation law with fluid-structure interaction prob- lems were discussed by Shyy et al. [88]. 2.2 Turbulence Models 2.2.1 Reynolds Equations and k-E Turbulence Model For simplicity, the Reynolds averaged momentum and continuity equations for incompressible flows are presented in Cartesian coordinates au, a(uJ U) aTj (2.4) at + x (2.14)xj 0, (2.15) where Ti is the stress tensor, repeated indices in any term indicate a summation over all three values of the index, and Ui is the Reynolds averaged mean velocity. The stress tensor in a turbulent flow may be written as T, = -P6,j + 2pSj ,i, i,. (2.16) 6ij is the Kronecker delta, which is equal to one if i= j and zero otherwise. P is the mean hydrodynamic pressure and p is the dynamic viscosity. The mean rate of strain Sij is defined by 1 9U, au. Si ~ + a ) (2.17) 2 Oxj Oxi ui and uj are the fluctuation velocities in the turbulence flow which contribute to the mean stress by T j -7, ,., (2.18) Tij is the Reynolds stress tensor which is symmetric. The diagonal components of Tij are normal stresses which make up the turbulent kinetic energy k k ui (2.19) 2 In many flows, these normal stresses contribute little to the transport of mean mo- mentum. The off-diagonal components of ij are shear stresses which pl ', a dominant role in the mean momentum transfer by turbulent motion. Eq. (2.14) contains nine unknown components of ij (of which six are independent of each other) besides P and the three components of Ui. Therefore, we need more equations for ij to close the system of equations. A turbulence closure model serves that purpose. Different versions of turbulence closure models have been proposed. The success of such models is largely depend on the applications of the wall functions that relate surface boundary conditions to fluid points away from the boundaries. Thereby we avoid the problem of direct modelling the viscosity influence. Here we use the wall function treatment by Launder and Spalding [52] as one example to illustrate the closure problem. The governing equations for standard k-E models can expressed in the form shown below, a(pk) a(pkUi) 9 ( it 9k t+ ) ax l ) + 2tSijSj pE, (2.20) at axi 9xi Jk xi 9(pE}) 9(pEU) a Pt E E E2 (p + +( 2CIptSijSij C2p. (2.21) at + xi ax"i aE xi k k The eddy viscosity pt is defined as k2 pt =pC-, (2.22) where C, is a dimensionless constant. Sij is the mean rate of strain. These equations have five adjustable parameters, CO, Jk, E, C.1, and C.2. Here, for brevity, we only give the coefficients for standard k-E model -ir_-'-. -1. I by Launder and Spalding [52]. These parameters are valid for a wide of range of turbulent flows. C, 0.09, rk = 1.00, a 1.30, C1 = 1.44, CE2 1.92. (2.23) To close the system of equations for turbulence computation, we need an additional Boussinesq relationship that relates the Reynolds stress to the mean rate of strain and turbulence kinetic energy, -, Uj = ( + O ) -pkJij 21 tSij pkbij. (2.24) axj axi 3 3 2.2.2 Boundary Conditions for Turbulence Models The elliptic characteristic of the k and E model equations gives rise of bound- ary conditions. Usually we assume zero gradient at the outlet boundary. At inlet boundary, if there is no measurement of k and E at disposal, one way is to make an approximation based on the freestream velocity U,, turbulence intensity Ti, and the characteristic length scale L. In our study we adopt the following formulas to define the freestream boundary conditions for k and E, k (UoTi)2, (2.25) 2 C3/4k3/2 c 4 2 (2.26) 0.07L The turbulence intensity Ti is usually set between the range of 0.02 and 0.05. Another way is to adjust the turbulence Reynolds number Ret=pUL/pft. For practice we usually set the turbulence Reynolds number in the range of 100 and 500. The solid wall boundary conditions are more complicated. Near the wall bound- aries, the local Reynolds numbers are low and the viscous effects are important. However, the k-E turbulence model that we mentioned above is based on the high- Reynolds-number assumption in which the turbulent stress forces are dominated. Moreover, close to the wall, extremely fine grid is usually required to capture the steep variations in flow profiles to solve the Navier-Stokes equations directly. The CPU time based on such a fine grid is almost beyond the ability of current facilities. In order to solve these difficulties, two approaches are commonly emploi- l1 one is the wall function treatment [52], the other is the low-Reynolds-number versions of the k-E [46, 73]. The wall function treatment is popular in practice because of its simplicity and capability for complex geometries. Under the high Reynolds number condition, flow can be largely divided into two 1-v. rs: the outer 1-v -r and the surface -1i.-r. In the outer l-?v-r the Reynolds stress is dominated; while in the surface l-?v-r both the Reynolds stress and the shear stress are important. In the outer l1.v-r, the fluid field is influenced by the retarding effect of wall through the value of the wall shear stress, but not by the viscosity itself. In the wall function treatment, velocity profiles have different descriptions in the out and surface 1-~v.i--. The velocity in the out -1~ .-r is governed by the velocity-defect law expressed as max U g( ), (2.27) U, 6 where 6 is the boundary l-i-r thickness, Umax-U is the velocity deficit, g is an unknown function with the order of 1. The velocity in the surface l-v-r is given by the law of the wall UT = f(y ). (2.28) Here y+ is a non-dimensional normal distance from the wall, and it is represented by y+ = '-,. (2.29) P- In the above expression y is the distance from the solid wall, and u, is the friction velocity at the wall defined by U, = (2.30) The shear stress at the wall, 7-, is given by dU P (2.31) The low-Reynolds-number k-E model requires a relatively fine grid resolution near the wall regions. For detailed descriptions we refer to Refs. [46, 73, 87]. To make appropriate use of the two versions of turbulence models, we strive to adjust the value of y+ of the first point next to the wall to satisfy certain requirements. For the low-Reynolds-number approach, y+ should be less than 2 for an accurate result, while y+ in the range between 30 and 500 [103] can produce accurate results for the wall function treatment. An equation for the flat plate boundary l -.-r theory is used in the grid generation process to get an estimation of proper cell size for y+, + = 0.172 []Re0O9. (2.32) It also should be remember, while the position of the first cell is critical, it must have enough cells inside the boundary l~-.-r to resolve the steep variations of flow variables. The success of the recent turbulence closure is restricted to the sufficiently large enough Reynolds number turbulence flows. Further refinement is needed before they are used with confidence to calculate near-wall and low Reynolds number flows. A more detailed review about turbulence models for near-wall and low Reynolds number flows is given by Patel et al. [73]. 2.3 Laminar-turbulent Transition A reliable computation of the boundary 1~.-r transition remains one of the most challenging demands in aerodynamics. Since the transitional separation bubbles and their losses p1 i, an important role for the pressure and drag distributions, an ac- curate representation of both laminar and turbulent separated flows is critical when drag prediction needs to be accurately predicted. However, the fundamental fluid dy- namics problem of transition to turbulence has resisted any simple phenomenological description. There is still no comprehensive theory of transition. The range of existing transition prediction methods extends from simple empiri- cal relationships through those based on parallel and linear stability theories, like the eN methods [15, 75], or linear or nonlinear parabolized stability equations (PSE) [31], or the low-dimensional models [76], to direct numerical simulation, like large-eddy simulation techniques. CHAPTER 3 MOVING GRID TECHNIQUE For problems involving interactions between fluid flows and moving structures, the geometry of the solid object is not known a priori. In the course of computation it is necessary to track the geometry, the field equations, and the interfacial condition to ensure that all requirements are met. To facilitate the solution of such moving bound- ary problems, the moving grid technique (or dynamic grid technique) is employ, .1 to adjust the computational grid dynamically along with the geometric updates. It is desirable but difficult to develop an automated regridding procedure to ensure that the new grid not only matches the geometric changes but also maintains satisfactory characteristics such as non-ovi. iI1 pplli, smooth, and not excessively skewed. The moving grid technique can also be used by shape optimization [54]. For moving grid problems, several alternative approaches have been proposed. Schuster et al. [81] used an algebraic shearing process in their static aeroelastic anal- ysis of a fighter aircraft. The basic idea is to assume that the displacement of the subject is redistributed along the grid line, which connects the aeroelastic surface to the outer boundary. This simple method works quite well for single block grid with modest deflection; however, when applied to multiblock grid arrangements, this approach usually needs extensive interventions. The spring analogy method was first proposed for unstructured grids [3], and later was extended to structured grids by Robinson et al. [78]. This method regards all grids as connected by springs, and each spring has the stiffness inversely proportional to the spacing between target vertex and its neighbors. Compared with other currently used grid regeneration methods, the spring analogy method needs more memory and CPU time. Direct transfinite interpolation by Eriksson [18] is a popular method because of its efficiency for struc- tured grid. Hartwich and Agrawal [29] also proposed the Master/Slave concept to expedite the grid regeneration process and minimize the user intervention. 3.1 Perturbation Method For small amplitude displacement on one face of the computational block, a simple but efficient one-dimensional perturbation method [56, 77] is can efficiently regenerate the grid. The one-dimensional perturbation method obeys the following formula, Xnew X od + old (new Xld), (3.1) where xi represents the location of an arbitrary interior grid point, x, represents the location of a grid point on the moving boundary, and S represents the normalized arc length along the radial mesh line measured from the outer domain, which is given by e lil + (+1 yi2 + (1+1 Z)2 Si = (3.2) Z7 (x xI)2 + (y1 yi + (z1+1 () 2 To use this one-dimensional perturbation method, one needs to know the displace- ments of the two ending points of a mesh line. The positions of the interim points are solely determined by the displacements of these ending points and their original positions as shown in Eq.(3.1). These displacements are considered as perturbation sources. Intuitively, the perturbation will spread throughout the whole domain while exerting more effects on nearer points. Based on Eq.(3.1), a more complicated 3-stage perturbation method, analogous to the transfinite interpolation (TFI) method [18], was proposed for three-dimensional problems treated by the single-block approach. Unlike the original TFI method, this perturbation method considers the original grid distribution, and preserves the original grid quality. Take a two-dimensional case as an example, in Fig. 3-1, a rectangle undergoes shape change. To regenerate the grid for the changed geometry, the perturbation method proceeds as a 2-stage process. The first stage is to match corner points, i.e., to move the corner point A to its new position A', and B to B', etc. (See Fig. 3 -b). Once the corner points match their new positions, interim edges are generated based on Eq. (3.1). To adjust the position of an internal point O, effects from the four corner points need to be accounted for in a coordinated manner. In this approach, P', Q', S', and T', which connect 0' to its boundary, will be first adjusted. The displacements of the four edge points in, e.g., the x-direction (AP, AQ, AS, AT) can be calculated with Eq. (3.1), and the movement of 0' is evaluated as a combination of these displacements. Specifically, one can adopt the following expression: A' abs(AI)AI + abs(AJ)AJ max(abs(AI) + abs(AJ), 10-6) where, Al = (1 Si)AP + SiAQ, AJ (1 Sj)AS + SjAT, Si and Sj are the normalized arc lengths along mesh lines PQ and ST respectively. Displacements in other directions can be computed in the similar way. After the first stage, the four corner points have matched their new locations; however, the interim edges do not necessarily match their final positions. The second stage is to match the interim edges to their final edge locations. The difference be- tween the interim and final edge locations are taken to be the source of perturbation. In this stage, since the perturbations in both directions are independent, the contri- butions from all individual directions are added [See Fig. 3-lc]. For three-dimensional problems, a 3-stage process perturbation method was proposed in Refs. [56, 77]. 3.2 Master/Slave Concepts The aforementioned perturbation method can be implemented efficiently for sin- gle block and small displacement problems. However, when a multi-block grid is involved, this method needs to be enhanced by other techniques. The Master/Slave concept acts as a useful approach to maintain grid quality and to prevent potential grid cross-over. I- SC S' -P .. Q------- -- (a) o/ Q XA B A' B' SS' S c ' /(C P I(b) S" *o Q' Q T f ( A B A' ' T' T' Figure 3 1: Two stage perturbation method: (a) ABCD changes to A'B'C'D'; (b) Sstage: :': corner points; (c) second stage: matching edge points. D S Although unstructured grid has been used to solve fluid dynamics problems by different researchers [3, 26], here, our attention is limited to structured grid only. For a multi-block structured grid, CFD software often requires point-matched grid block interfaces. A basic requirement for grid regeneration in moving boundary problems is to maintain point-matched interfaces. While this updating process is easily satisfied by those block surfaces that coincide with the body surface, the remaining block surfaces need to be handled with due care. To use the perturbation method, one needs to specify the new positions of the off-body surface points, or, at least, the vertices of such block surfaces. When an object changes its shape/position, master points, which are defined as the points located at the body surface, first move, and then affect distributions of the slave points, the off-body points. While the task of redistributing interior grids can be (d ill. i-ii] in practice, it is more difficult to move the vertices of each block if a point-to-point matched interface without overlap is required. If two abutting blocks have congruent interfaces, such as the interface of block 2 and block 3 in Fig. 3-2a, once the corner vertices are determined, the edge points and the interior points can be uniquely settled based on Eq. (3.1). However, when two neighboring grid blocks do not share identical interfaces, such as block 1 and block 2 shown in Fig. 3-2a, this procedure can cause discontinuity at the interface because the regriding procedures at the interface are based on different stencils for different blocks. Slaving vertices to their master points can avoid creating undesirable grid discontinuities. The Master/Slave coupling is based on the distance between an off-body point (Slave point) and a surface point (\! ,Iter point). The distance is given by: |r| = [(X, Xb)2 + (Y, b) + (, )2]0.5, (3.4) where subscript v represents an off-body vertex, and, b, a body surface point. For each off-body vertex, the nearest body surface point is defined as its master point. (a) (b) (0,1.48,0) BloBlock ck 3 Block 1 (2,0.9,1) (0,0 B lock 2 \(2,0.5,1) (0,0.5,1) (c) (d) It. ;'lii (e) Figure 3-2: Tests on the moving grid scheme: (a) initial geometry; (b) side view of the grid; (c) deformed shape; (d) new grid distribution; (e) grid distribution with high stiffness parameter. To simplify the connection between the master and slave points, one can convert the three-dimensional parameter space (i, j, k) into a one-dimensional data structure. Each grid point has an identification number 1 defined by I = boff(bn) + i + imax X (j 1) + imax X max X (k 1), (3.5) where boff(bn) is the identification number of the first point in block bn in the one- dimensional array, imax and jmax are the maximal dimensions in the i and j directions, respectively. In this way a slave point is associated with its master point by storing its master's identification number into its address. Once the relation is established, the next step is to determine how the slave point moves according to the influence from its master point. Intuitively, a nearer vertex will have more effect (displacement) than those far away. Also, to avoid a lower surface to cross over its opposite face when the movement is relatively large, one needs to scale the movement. A simple but effective formula is expressed as follows: Xs = Xs + 0(X Xr), (3.6) where subscripts m and s represent the master and slave points, respectively, the tilde (~) indicates the new position, and 0 is a decay function. In our computations, a Gaussian distribution function is used 0 = exp{-3min[500, dv/(E +dm)]}, (3.7) where 3 is a parameter which affects the stiffness. dv = (x, Xm)2 + (yv m)2 + (z, z)2, (3.8) dm = (xm X,)2 + (,m- Vm) + (Zm- __)2. (3.9) 3.3 Numerical Tests The moving grid technique consists of two elements, one is the perturbation method, and the other is the master/slave concept. The proposed method is tested on a three-dimensional, multi-block, and structured computational grid. The initial configuration is shown in Fig. 3-2a, block 1 shares the interface with blocks 2 and 3. The initial grid is shown in Fig. 3-2b, it has clustered distributions at the top and bottom surfaces, and it also has point-to-point matched interfaces. The bottom two faces are the moving boundaries that experience large deformation. The coordinates of each vertex are labelled. Comparing Fig. 3-2a and Fig. 3-2c we can see the displacement is as large as half of the total height. The new computational grid is computed with the method discussed in this chapter. The new grid is shown in Fig. 3-2d. That figure does not show any cross-over in the new grid after the large deformation. The new grid not only keeps the point-to-point matched interfaces, but maintains the clustered distributions at the top and bottom surfaces. The CPU time for regenerating the computational grid is less than 10' of the CPU time used by the fluid solver. The parameter 3 controls the stiffness, we find /3=1/64 works well for our prob- lems. A larger 3 causes the block to behave more like a rigid body (See Fig. 3-2e). More details about the master and slave concept can be found in Refs. [29, 56]. CHAPTER 4 ACTIVE FLOW CONTROL WITH AN ADAPTIVELY-SHAPED WING 4.1 Introduction The main objectives for fluid dynamists and aeronautical engineers throughout the history of modern flight have ahlv-b-, been to increase lift and reduce drag. Micro air vehicles are airborne vehicles with characteristic lengths under 6 inches traveling at speeds around 30 ft/s. They are designed to perform an array of different tasks such as reconnaissance, target designation or border surveillance. To enhance the performance of these small vehicles, improvements should be made in the respects of power plants, thrust generators, navigation systems, cameras, data transmitting, and more importantly, their aerodynamic performance. As far as the aerodynamics is concerned, MAVs have to operate in one of the worst flight regimes imaginable, one of them is the low-Reynolds numbers which is usually norotiously accompanied with laminar boundary separation, transition, and low lift-to-drag ratio. Furthermore the small dimensions of the MAVs greatly enhances the trade off difficulties mentioned above. The small size inhibits the designer's freedom to achieve the desired aerody- namic performance. The limitations in the power sources available, at present time, make this equation even more complicated. To solve the aerodynamic challenges of these problems, two in ii r fields have been evaluated. The first approach is trying to copy nature's own way; that of flapping wings [43, 45, 59, 83]. The second approach is using a dynamical wing There are some different implementations in this field. One widely studied method is the use of an wing with a flexible membrane upper surface [53, 55, 83, 92]. Another method uses active devices to control the flow on the upper surface of the wing. The latter method has gained more attention with the success microelectromechanical systems ( \!i \!S). In active flow control area, a piezo-actuated flap mounted on a low Reynolds number wing have been studied both experimentally [19] and computationally [30, 56]. Gad- el-Hak [20] has discussed the possibility of employing MEMS to delay or prevent separation on MAVs. To understand the low Reynolds number wing characteristics and improve the wing performance, we will perform a computational investigation for flows surround- ing a dynamically shaped wing, at a chord Reynolds number of 78,800, along with a parallel experimental effort. We adopt a piezo-actuated flap on the upper surface of a fixed wing for active control. The actuation frequency they focus on is 500 Hz. Their computational framework consists of a multi-block, moving grid technique to handle the geometric variations in time, using two-equation turbulence closures, and a pressure-based flow solver. The dynamic grid redistribution employs the transfinite interpolation scheme with a spring network approach. Comparing the experimental and computational results for pressure and velocity fields, effects of the detailed ge- ometric configuration of the flap, the flapping amplitude, turbulence modeling, and grid distributions on the flow structure will be assessed in this chapter. 4.2 Experimental and Computational Set-up The wing has moving flaps on the upper surface as shown in Fig. 4-1. The flap can be actuated at prescribed frequencies [19]. The wing studied in the experiment has a span of 1 ft and a chord of 5.4". In the experimental set-up, there is a gap of 0.016" wide between two flaps, each is 2" wide and 0.564" long. The Reynolds number, based on the freestream velocity U=27.5 ft/s, is 78,800. The flow is initiated as laminar; after the transition point, the turbulence model is invoked. The transition point is identified to be at the tip of the flap. In the computation, the turbulent Reynolds number based on the chord length, Ret=pUL/pt, at the inlet of the turbulent flow region is set to 500. The amplitude and shape of the flap are prescribed by analyzing the experimental data obtained with PIV [19] For the 500 Hz case, two flapping amplitudes are considered, namely 0.001" and 0.0073". region of interest piezo actuator (flaps) for PIV \ / \o / pressure taps, 0.25 in. apart A X Figure 4-1: Moving flap wing experimental setup in the wind-tunnel In this chapter, a number of two-dimensional and three-dimensional cases, in- cluding different grid densities, with and without the gap, will be studied. For 3-D cases, an effort has been made to consider the effect of the gap between flaps as shown in Fig. 4-2. The original Launder-Spalding k-E model [52] with the wall function, as well as a low-Reynolds number version k-E model [46] are compared. Linear Chords - Figure 4-2: Flap wing with refined three dimensional bending model. 4.3 Numerical Study 4.3.1 Grid Sensitivity Analysis A number of 2-D grids with different distributions and sizes have been tested. Then a reference 2-D grid is selected by ensuring that satisfactory solution quality is obtained on a most economic grid size. Then, the 3-D grid is established by extending the reference 2-D grid along the spanwise direction. The resulting grid, for cases with gap, has 33 nodes in the spanwise direction, and consists of 20 blocks with approximately 533,000 grid points. For cases without gap, the grid system has 20 nodes in the spanwise direction, and consists of 8 blocks with about 253,000 grid points. With the multi-block strategy it is convenient to employ the transition model, i.e., when a transition location is identified, new computational blocks are generated surrounding the transition point, so that the laminar flow equations are solved in the blocks before the transition location, and the Reynolds-averaged turbulent flow equations are solved in the blocks after the transition point. In our computations, the transition point is identified from the experiments [19], and it is at the tip of the flap, therefore, our computational blocks are set up based on this transition point. 4.3.2 Assessment of the Turbulence Model In the early work by He et al. [30], the wall function approach is used exclusively. For such a low Reynolds number, the low-Reynolds number version of the k-e model seems more suitable [56]. Fig. 4-3 compares the time-averaged pressure distributions at 500 Hz between two different turbulence models for the case without gap between flaps. Both experimental and computational data are collected at the middle of the flap and along the chordwise direction. As far as the computational time is concerned, it is observed that the low-Reynolds number model uses less CPU time than that of the wall function treatment. 4.3.3 Effects of Gaps between Flaps The inclusion of a gap obviously increases the complexity of the flow field in the spanwise direction. It can be seen from Fig. 4-4 that streamlines enter the gap and move toward the middle of the flap. A weak recirculatary flow with a vortex stretching from the wing surface to the flap is present on both sides of the gap, rotating in the clockwise direction on the right and counter-clockwise direction on the left of the gap when viewed from above. The gap also affects the overall pressure distribution on the wing. Fig. 4-5 compares the pressure distributions obtained with and without gap between flaps. The low-Reynolds number turbulence model is used in both cases. Time average is - Low Re Model SExperimental -- Wall function model 0.5 0.4- 0.3 - 1 0.2 0.1 1.5 2 2.5 3 3.5 4 Figure 4 3: Comparison between tine-averaged pressure distributions for cases without gap, at 500 11z, obtained with different turbulence models. (: : : for the entire cycle, .i.' ude: I01"). -\ ' Figure 4 4: Streamline for 3-D case with gap at 500 Hz, amplitude is 0.'-i1" made for the entire cycle for an amplitude of 0.001" and an actuation frequency of 500 Hz. Since the pressure taps are located at the middle of a flap, the effect of the gap may not be clearly demonstrated by the figures shown. Fig. 4-6 compares the velocity profiles immediately downstream from the flap trailing edge. The low- Reynolds number turbulence model is used again and the comparison is made at 0.0162" downstream from the trailing edge of the flap. Cases involving both stationary flaps and actuated flaps, at 500 Hz, are examined. Consistent with that already observed, the velocity profiles indicate that the flow enter into the gap from the upper surface of the flap, maintaining the same streamwise direction while moving laterally toward the region between the flap and the wing surface. A recirculating flow pattern is formed between the flap and the wing surface. 0.7 No gap between flap 0.65 Experimental .--- Gap between flap 0.6 0.55 0.5- O 0.45 0.4 0.35- 0.3- 0.25 - 0.2- 1 1.5 2 2.5 3 3.5 4 X(in) Figure 4-5: Comparison of time-averaged pressure distributions for cases with and without gap between flaps. 4.3.4 Effects of Flapping Amplitude Depending on the flapping amplitude, the reattachment point can move back and forth, and the recirculating flow can either remain attached to the tip of the flap, or be shed away. In our computation the curvature of the flap is set by collecting 10 5 10' 105 10- 10- 95 95 / 9 9- S85 / At the gap 85 S- Middle of the flap 0 - 75- 75- 7 /7 655- 65 6 5 6 At the gap Middle of the flap 55 55 2 0 2 4 6 8 10 12 -5 0 5 10 15 U velocity foot/s U velocity foot/s Figure 4-6: Comparison of velocity profiles at two different spanwise positions with gap implemented: (a) steady computation; (b) f=500 Hz. the experimental data and then approximated by a simplified model. For the 500 Hz case in our computation two amplitudes are considered, one is 0.001" and the other is 0.0073". For the smaller amplitude case, the motion of the trailing edge is modeled as the following, Fxz y(z, t) = (8.6 x 10-4 + 1.6 x 10-4 sin( )) sin(27f t), (4.1) 2 where z is the spanwise distance from the study point to the gap, f is frequency, and t is time. The larger amplitude case simply scales the above equation based on the same expression. Fig. 4-7 compares the time-averaged pressure coefficients between the two flap- ping amplitudes. For the higher flapping amplitude, the time-averaged reattachment point moves closer to the trailing edge of the flap, because a larger movement of the flap introduces more high momentum fluids into the boundary 1 r. In the experi- ment we observe vortex shedding as shown in Fig. 4-8a. In our computations with the smaller flapping amplitude, the vortex core is anchored at the flap tip, exhibiting two co-rotating cores as the flaps move upwards, as shown in Fig. 4-8b. On the other hand, with the higher amplitude, vortex shedding appears for cases with or without 0.5 0 Disp=l.0x 10-3 in SDisp=7.3 x 10-3 in 0.45- 0.4 0 0.35- 0.3 0.25 0.2 1 1.5 2 2.5 3 3.5 X position(in) Figure 4-7: Pressure coefficient at different actuation amplitude without gap between flaps. Actuation frequency: 500Hz. the gap. Fig. 4-9 shows selected instantaneous plots of the vortex structure during a single cycle of the flap movement. As expected, the recirculation is stronger when the flap moves upward and becomes weaker when the flap moves downward. Figure 4-8: Instantaneous trailing edge streamlines for the frequency of 500 Hz as flap is moving upwards for amplitude of 0.001". Left: experimental result; right: computational result. 4.4 Conclusions Base on our results we find that the low-Reynolds number formulation of the k-e turbulence model produces longer circulation zone. The gap between the flaps affects the location of the re-attachment point, and creates a three-dimensional flow pattern. Aided by the limited experimental information, the computational model can offer t Trailing edge of iap (a) Phase angle a = -/6. (c) Phase angle a 7i /6 (b) Phase angle a = r3. (d) Phase angle ao 8-r/6. (e) Phase angle a = 11T/6. Figure 4-9: Trailing edge streamlines with an ..:.. i'tudc of 0.C" ' (f) Phase angle a = 2-. . the f- : ..-- of 500 Hz in the third cycle insight into the fluid physics. Vortex dynamics changes substantially, in accordance with the flapping amplitude. Vortex shedding occurs when the amplitude is higher, while anchored vortex is observed with the lower amplitude. As to the physical mechanisms, while it is expected that the vortex shedding appears when the amplitude is high, it is not clear why and when this phenomenon ceases to exist when the amplitude is reduced. Similarly, while it is obvious that the gap between flaps affects the flow reattachment, it is hard to predict the trend quantitatively. CHAPTER 5 MEMBRANE STRUCTURAL SOLVER 5.1 Introduction The deformations of membrane structures are often large and thus inherently nonlinear. Qualitative descriptions of the thin membrane behaviors are often more complex than that of three-dimensional bodies. Until recently, membrane analyses were primarily relied on trial and error. Green and Adkins [25] are the first to derive a general non-linear theory for membranes. The membrane theory has two fundamental assumptions because of their thinness [41]. First, the membrane tension stiffness is much larger than its bending stiffness, therefore, the stress couple effects can often be neglected. Second, the ratio of the membrane thickness to its radius of the curvature is small, which effectively decouples the strain-displacement relation from the curvature tensor. Practically speaking, there exist two approaches to model the response of a membrane structure. The first approach is to idealize a plate or shell structure as a membrane away from its boundaries, in which the stress couples are neglected in the interior region. By formally equating the plate stiffness to zero, this membrane model can be described by the Foppl-von Kdrmdn theory [40, 41]. In the second approach, the structure cannot sustain stress couples anywhere. The second approach is adopted in the current work because it is appropriate for MAV computations. Two-dimensional membrane wing study with fixed leading and trailing edges were originally conducted by N. 1- [68], Thwaites [101], Voelz [108]. These works considered inextensible membranes surrounding by steady, two-dimensional, and ir- rotational flows. As a consequence of the inextensible assumption, when cambers and incidence angles are small, problems can be linearized and expressed compactly in a integral form which is referred as the "sail equation" by N. i.--in ,, [69]. The following equation completely defines the linearized theory of inextensible membrane wings in steady and inviscid flow fields. d2(y/a) c, d<2 < d(y/a) 1 J- 0 d( = (5.1) 2 o 27r( -) drx where y(x) defines the membrane profile as a function of the x coordinate, a is the incidence angle, and CT is the tension coefficient. Based on the work of Voelz [108], the sail equation has been successfully solved analytically and numerically. N, ,.-i in. [69] gave a comprehensive review of the earlier works related to membrane wing aerodynamics. The development of a three-dimensional membrane model introduces several complicating factors not encountered in two-dimensional analyses. First, unlike the - il equation", the tension in a three-dimensional membrane is no longer a single constant value but is at best described by a state of biaxial tension along the lines of principal stresses [38]. Second, geometric and material properties vary along the spanwise direction, and they have to be described accurately. Third, if one of the principal tensions vanishes, the membrane will be compressed and wrinkled, which further complicates the analysis. Considering the complexity of the three-dimensional problem, simplifications are made in the early works. Oden and Sato [70] presented a finite element method for the analysis of large displacements and finite strains in elastic membranes of general shape. In their work the static equilibrium of an inflated membrane was considered. A static analysis can be considered as a special case of a dynamic analysis in which the inertia effect is included. Works on the dynamic analysis of membranes before 1991 were reviewed by Jenkins and Leonard [41]; most recent works were reviewed by Jenkins [40]. Verron et al. [105] performed a combined numerical and experimental study of the dynamic inflation of a rubber-like membrane. Membrane wrinkle problems were tackled in the work of Ding et al. [14]; they performed a numerical study of partial wrinkled membranes. The stability issues of a fluid flow past a membrane were reported by Thaokar and Kumaran [99]. In this chapter a three-dimensional membrane model is proposed for the mem- brane dynamics. The membrane material considered obeys the hyperelastic Mooney- Rivlin model [64]. We use the Green-Lagrange strain tensor for the description of large strains. The dynamic response of the membrane is described by a system of second-order ordinary differential equations. A finite element method with triangu- lar elements for spatial discretization is utilized. For time integration, the implicit Wilson-O method is used. 5.2 Structural Solver 5.2.1 Mooney-Rivlin Model For an initially isotropic membrane, Green and Adkins [25] showed that there existed a strain energy function W which was expressed as the following form W =W(II, 12,1), (5.2) where 1i, 12, and 13 are the first, second, and third invariants of the Green deformation tensor C. If the material is incompressible, namely, 13 1, then the strain energy is a function of IA and I2 only. It can be written as: W -= ci( 3) + c2(2 3), (5.3) where ci and c2 are two material parameters. A material that obeys Eq. (5.3) is known as the Mooney-Rivlin material [64], and the model is called the Mooney-Rivlin model. The Mooney-Rivlin model is one of the most frequently employ, ,1 hyperelastic models because of its mathematical simplicity and relatively good accuracy for reasonably large strains (less than 150 percent) [64]. For an initially isotropic membrane, the general stress-strain relationship is writ- ten as 8W S = -pC- + 2 (5.4) where S is the second Piola-Kirchhoff stress tensor, p is the hydrostatic pressure. If the membrane is incompressible, then the stress-strain relation can be further simplified as the following S -pC-1 + 2 [(cl + c21)I c2C] (5.5) where I is 3-by-3 identity matrix. In the membrane theory the stress field is essentially treated as two-dimensional, and the stress normal to the deformed membrane surface is negligible in comparison with the stresses in the tangent plane [70]. Under this consideration, the deformation and stress matrices, C(t) and S(t), are given by CiI(t) C12(t) 0 S11(t) S12(t) 0 C(t) c 12() 22 () 0 S(t) 12(t) S22(t) 0 (5.6) 0 0 C33(t) 0 0 0 Therefore, the hydrostatic pressure is determined by considering the condition that S33-0, p =2(ci + ca21 C2A2)A (5.7) where A3 is defined by A=3 C33(t) hw h(5.8) in which h(t) and ho are the membrane thickness in the deformed and undeformed configurations, respectively. 5.2.2 Principle of Virtual Work Consider an elastic membrane, which in its initial state occupies a finite region in space. The position of the region can be specified by an orthogonal coordinate system XYZ which is referred to as global coordinates. In global coordinates, the general form of the principle of virtual work in the Lagrangian description is given by f DT (t)poD(t)dV + E: SdV- f o DT(t)P(t)dS = 0, VOU, (5.9) where Vo and iOVo are, respectively, the volume and the boundary surface of the undeformed membrane, po is the membrane density, D(t) is the displacement vector in global coordinates, 6D(t) is the virtual of vector D, the superscript T refers to the transverse of a vector or a tensor, D(t) is the acceleration vector, P(t) is the generalized external surface force, E is the Green-Lagrange strain tensor, S is the second Piola-Kirchhoff stress tensor, and 6E:S is the scalar product of two tensors. In the rectangular Cartesian coordinates, the scalar product is defined by E : S =6EijSji, (5.10) with summation on the repeated indices. 5.2.3 Triangular Membrane Element The principle of virtual work can be written in the following discretized form 61I (f DTpoUdV + 6E: SdV / DTPdS =0, V6D, (5.11) i \JVo0 V0 ,, where 61 is the total virtual work, and the subscript e represents elementary quan- tities, N, is the number of elements covering the entire membrane. In the present work the triangular element is employ, As shown in Fig. 5-1, a typical triangular element, e, is defined by nodes, i, j, k, and three straight-line boundaries. In the global coordinates the three nodes have coordinates (Xi, Y1, Z.), (Xj, Y, Zj), and (Xk, Yk, Zk), respectively. To facilitate the development of a computational procedure, it is convenient to consider the membrane element and its displacement in the local coordinates. As shown in Fig. 5-1, in local coordinates, the triangular element lies in the plane of k' y=x+u Y k x 1 z Y X Figure 5-1: Schematic global and local coordinates. ;, Nodes i, j, and k, have local coordinates of (xi, yi, 0), (xj, yj, 0), and (xk, yk, 0), respectively. For each point within the triangular element, its local coordinates x and global coordinates X have a relation expressed as follows: X X + Tx, (5.12) where Xi is the global coordinate of node i, and T is an orthogonal transformation matrix which satisfies TT= T-1. Similarly the displacements have the following relation D =Td, (5.13) where D is the value of a displacement in global coordinates, d is its value in local coordinates. If the element is sufficiently small, we can approximate the displacement at any point within the element with a linear combination of the displacements of nodes, i, j, and k, u = N3d, (5.14) where N3 is an interpolation matrix defined as 1- -I 0 0 0 0 r 0 0 N3(,r)= 0 1 0 0 00 0 (5.5) 0 0 1- 0 0 0 0 rI where 0 < ~< 1, 0 de di dj dk (5.16) in which di, dj, and dk are the nodal displacement vectors of nodes i, j, and k, respectively. Each node has three degrees of freedom, di= u (5.17) 5.2.4 Green-Lagrange Strain Tensor When a displacement is introduced the location of a point originally given by coordinates x now has new coordinates y x + u, (5.18) where u is the displacement as shown in Fig. 5-1. Following Green and Adkins [25], the Green-Lagrange strain tensor is defined by the relation 1 ( Oyk Yk (5.9) E =i -xi -xj 5i (5.19) Since the thickness of the membrane is small compared to the other characteristic dimensions, the following simplifications are allowed to be made [25, 70] _1 (Qu 8 9uP 9ui 9ui\ E t ( -, + U-- + N O 2 x, 3 dx dOx, dxo E3 = 0, (5.20) E33 2 (A 1), C, = 1, 2, 2 where A is the same as A3 in Eq. (5.8). To simplify the computation, we divide the Green-Lagrange strain tensor into two parts. One part is Eo, which takes into account small strains, u 1 iu av + )O ax 2 O y ax 1 au av av Eo (= t + )v 0 2 ay Ox ay 0 0 0 and the other is EL, which is the nonlinear part of the strain tensor, ou u Ovu v Owv w Ouw u Ovu v Ov w Ow Ox Oxox ox x Oxx Ox Oy Ox 9y Ox Oy SOuu + Ovv v Oww Ouw u u Ovv + Oww 0w EL = + + + +0 Ox ay Ox y Ox y OyOy ay y ay y y 0 0 A2 (5.21) (5.22) To find each component in the strain tensors Eo and EL, it is necessary to compute the displacement gradient uV(x,y). After some operations we obtain the expression of displacement gradient, Su Ou ui uj Xk(ui j) ui Uk Ox Oy Xj Xjyk Yk a-- v av V Zk(X ) 5. -2) UV (,) (5.23) ) x y Xj xjyk yk aw aw wj Xk(I' Wj) i ii' , Ox Oy xj xjyk Yk Based on the expressions in Eq. (5.23) we can compute the components Eo and EL. 5.2.5 Internal Forces When large deformations and/or nonlinear material properties are involved, the equivalent internal elastic resisting forces of the structure can be expressed as the following, F'nt j BT((d))aT(e)dV, (5.24) JVoe where a is the nonlinear stress vector, e is the strain tensor, and B is the strain- displacement matrix which is also a function of the displacement for large deforma- tion. Recall that in the principle of virtual work the term related to the internal forces from the discretized form of principle of virtual work Eq. (5.11) reads S 6E : SdV. (5.25) Eq. (5.6) tells us that the membrane stress is essentially two-dimensional, therefore, we know from Eq. (5.10) that it necessary to consider the reduced strain components instead of all the components when calculating the internal forces based on Eq. (5.24). The general membrane strain vector can be conveniently decomposed into two parts according to the strain decomposition in Eqs. (5.21) and (5.22), E = o + L, (5.26) where so is the reduced strain vector for the linear part of the strain tensor, Eo, and EL for the nonlinear part EL. Based on Eqs.(5.21) and (5.22) these two strain vectors can be expressed as the following, S u 'u Ou Ovv Ov aw w Once the virtual displacements are given at any point within the element, the strain variation at any point can be determined. The relationships between the virtual oxOx Ox Ox Ox Ox OX OX displacement and strain variation arev v < 60 = Fgy B0de, 6L { j Ly } BLd, (5.28) 9u +9Ov 9i2-u u +9v v +9Ow Ow 6x,,, 5ELay f< y Ox Ox 9y Ox 9y Ox 9y , where and B are the stradisplacements are given at any point wicesthin the expressions can be straifound variation Appendix A. Based on Eq. (5.27) the virtual work in each triangular element canelationships between the virtual displacement and strain variation are 6Eo 6FOy Bo6de, 6L 6ey BMde, (5.28) 6-11, 6FLxy where Bo and BL are the strain-displacement matrices, their expressions can be found in Appendix A. Based on Eq. (5.27) the virtual work in each triangular element can be computed as / 6E : SdV = 6DT3 f BTa dV = D TFnt, (5.29) where B = Bo + BL, 6eT = Le + bg, and Fnt is the internal force in global coordinates Fi't = TS BTa dV. (5.30) T3 is the a 9-by-9 transformation matrix defined as T 0 0 TOO T = 0 T 0 (5.31) S0 T in which T is the same transformation matrix as shown in Eq. (5.13). 5.2.6 External Forces As the membrane changes its shape and position under the pressure, not only the direction of the pressure changes but also the area on which it acts changes [70, 105]. To account for this situation, we assume that the size of each finite element is sufficiently small that the pressure is uniform over each element. To further simplify the computation, we again use local coordinates. As shown in Fig. 5-2, the local coordinates lies in the deformed plane instead of in the undeformed plane. In the local coordinates the pressure P is given by, P = Pn, (5.32) with n [ 0 0 1 ]T is the external normal direction of the element. The virtual work by the external force at each element is D DT PdS, k X Zi 1 z / Y X Figure 5-2: Schematic of local coordinates for external forces. the net external force at each node is obtained by simply representing it by three equal forces at each node SDTPdS = DT3oP 0 0 0 0 3 0 0 1 6D Ft. (5.33) J n' ,, L3 3 3 Hence, in global coordinates the node force on element e due to a uniform pressure over that element is F T3OP o0 03 O 0 1 0 0 1 (5.34) 5.2.7 Dynamic Equations for the Membrane The first term of the left side of the principle of virtual work equation (5.11) can be written as VI 6DTpoDdVo D= DMeD, (5.35) Voe where Me is the mass matrix of an element, 21 I I 1 Me = posoho I 21 I (5.36) 12 I I 21 in which I is a 3-by-3 identity matrix. Up to now we have obtained the expressions for the internal force, external force, and mass matrix on an element e. It is necessary to assemble all the discrete elements into a single unit. For any element e the principle of virtual work is 6H1e = JDTpoDdV + 6E: SdV 6DTPdS. (5.37) Vo IVo ,, Assembling all the forces in one element we have the following equation, Se = 6D> (MD + F'nt Fext (5.38) This equation is for a single element, if we sum all elements in global coordinates, then we obtain the system of governing equations for membrane under external loads, MD(t) + F'nt Fext, (5.39) where M is a positive definite mass matrix, which remains constant, Fint is the internal force, and Fext is the external load. 5.2.8 Solution of the Dynamic Equations Either an implicit or an explicit scheme can be adopted to follow the time history of the system of equations. The explicit scheme is computationally economic for each time step and requires less storage. Its efficiency can be further improved by special techniques, such as lumping technique of the mass matrix. However, its time step needs to be small to satisfy the stability requirement. As far as the implicit scheme is concerned, the matrix system is updated more than once per time step to account for the nonlinear effect, which requires more computations per time step and more storage. Generally speaking, implicit algorithms are effective for structural dynamics problems in which the response is controlled by a small number of low frequency modes, while explicit algorithms are efficient for wave propagation problems in which the contribution of intermediate and high frequency structural modes to the response is important [94]. There are numerous v--v to integrate the structural dynamics equations in time; among these we utilize the Wilson-O method [113]. The Wilson-O method has a controllable algorithmic numerical dissipation, it is a one-step method, and it yields second order accuracy in time. These characteristics make it easier to code and behave well in our fluid/structure interaction system. Since the Wilson-O method is dissipative, we employ a time step smaller than that required by the overall accuracy criterion to avoid excessive dissipation. This approach, with 0=1.4, works well for the present problems. 5.3 Numerical Test of Membrane Model The inflation of a spherical membrane under an imposed pressure step was an- alyzed experimentally by Alexander [1] and theoretically by Beaty [4]. Recently, a numerical study was performed by Verron et al. [105]. The corresponding semi- analytical solution was reported by Verron et al. [104]. These studies focus on an incompressible spherical hyperelastic membrane. The membrane has an initial radius of Ro, an undeformed thickness of ho, and a mass density of po, respectively. The membrane material is assumed to obey the Mooney-Rivlin model. There are two con- stitutive parameters in this model, namely cl and c2. To simplify our computation, we define a non-dimensional Mooney-Rivlin constant based on these two constitutive parameters a =C2 (5.40) The analytical governing equation for the spherical membrane is detailed by Verron et al. [104]. Here we only give the dimensionless form =p*2 + ( A)(1 +a 2), (5.41) A5 where A is the circumferential principle extension, which is defined as the ratio of the deformed to the undeformed radius. The corresponding non-dimensional time r, and the non-dimensional imposed pressure step p* are given as follows: t 2 c, (5.42) Ro V Po P* = o* 4 (5.43) 4clho where Po is the dimensional pressure step. The initial conditions of the membrane may vary; here, we focus on the case that the membrane is initially in its equilibrium state, i.e., with the following initial conditions: A(r 0) = 1 and A(r 0) = 0. (5.44) One can obtain the semi-analytical solution to Eqs. (5.40) and (5.43) using any stan- dard solver. The direct numerical solutions of the membrane equation are obtained with the finite element method discussed before. The computational mesh has 816 nodes and 1568 triangular elements on a semi-spherical surface. In our work we choose At5 x 10-4 which is 500 times larger than the explicit scheme used by Verron et al. [105] where they use an explicit scheme with a time step of 10-6 for the same problem. All computations are performed using the double-precision. In Fig. 5-3 we com- pare the finite element results with the semi-analytical solutions. Except for one case in Fig. 5-3b, all numerical solutions match well with the semi-analytical solutions. The loss of linearity of oscillations is well reproduced. The mismatched case shown in Fig. 5-3b corresponds to a=0.1. For this value of the material parameter, as pointed out by Verron et al [105], there is an unstable equilibrium point at p*W0.687. Near the unstable equilibrium point, the finite element results can only approach qualitatively the real behavior. This is the reason why the numerical results are highly sensitive to the assigned pressure step. As expected, the larger the Mooney-Rivlin constant a is, the higher pressure the membrane can sustain. Two profound issues should be considered in the use of the Wilson-O method. First, the Wilson-O method is unconditional stable for linear problems. The time step in nonlinear problems is often dictated by the errors introduced in the approximation treatment of the nonlinearities; therefore, due care should be taken to choose an approximate time step. Second, a small time step is desirable for the Wilson-O method due to its excessive dissipation. As shown in Fig. 5-4, a larger time step tends to cause more numerical damping. In Fig. 5-5 we show the error norm of the implicit scheme, the figure exhibits a slope around 1. Therefore, the overall accuracy for time integration is approximately first order. The present membrane model gives a good agreement with the analytical solu- tions for membrane dynamics under large deformation. It leaves one issue of funda- mental importance since it cannot handle wrinkle phenomenon which happens when the membrane is compressed. The membrane wing is sometimes compressed, espe- cially at the leading edge. Therefore, further investigation is necessary in this area. p =0.6 I0 5 0 5 15 20 p=1 2 p=07 10 15 20 Figure 5-3: Dynamic inflation of a Mooney spherical membrane: (a) a = 0; (b) a 0.1; (c) a = 0.25. Solid line: Semi-analytical solution; Circle: Numerical solution. A t=2.0e-4 1.375 $1 ;~ < 1.25 1.125 10 15 20 25 0 5 10 15 20 25 Figure 5-4: The effect of time step to the numerical damping for p*=0.5 and a=0. -4.7 0) t-4.9 / 0 LU. -5.1 / -3.8 -3.6 -3.4 Time Step (login) -3.2 Figure 5-5: Error norm versus time step for the implicit scheme for p*=0.5 and a=0. 1.375 z 1.25 1.125 10 0 5 A t=5.0e-4 CHAPTER 6 MEMBRANE WING AERODYNAMICS FOR MAVS 6.1 Introduction Membrane structures are found in many engineering practices. Traditionally, they are utilized in parachutes and sail boats to improve performance. A recent interesting application is a membrane wing-based Micro Air Vehicle (\!A V). MAVs with a maximum dimension of 15 cm or less and flight speed around 10 m/s provide expendable platforms for surveillance and data collection in situations where larger vehicles are not suitable or else dangerous to reach. However, MAVs are sensitive to the environmental disturbances, and it is hard to avoid considering their small dimensions, light weights, and low flight speeds. Our experience indicates that an aeroelastic membrane wing can better adapt to the atmospheric disturbances, provide smoother viewing platform, and make the vehicle easier to fly [83]. During flight, the surrounding flow exerts aerodynamic forces on the membrane wing and causes structural deformations of the membrane surface; meanwhile, the wing deformations affect the flow pattern and change the aerodynamic force distri- butions. The mutual interactions between the fluid flow and the structure cause a fluid and structure interaction problem. Fluid and structure interaction problems have attracted interest of many researches. Numerous studies have been conducted on these subjects. Numerical schemes for fluid and structure interaction problems may vary. However, in general, there are two key elements in a coupled fluid and structure interaction system. One is the fluid solver that evaluates the aerodynamic forces, and the other is the structural solver that computes the shape deformation. Beside these, a moving grid technique is also needed to regenerate the computational grid; an interface algorithm is utilized to exchange information between the fluid and structural solvers. The interaction between a fluid and a membrane usually causes substantial struc- tural deformation that is inherently nonlinear. The membrane behavior is in general complex and earlier studies of interactions between fluid flows and membrane struc- tures are largely based on empirical observations [34], two-dimensional flow based computations [36, 86, 92, 93], or simplified three-dimensional analyses [37, 38]. In our study we present a coupled fluid and structure interaction system to investigate the interplay between a membrane and incompressible viscous flow. In the coupled system, the fluid dynamics is governed by the three-dimensional Navier-Stokes equa- tions for incompressible flows, and the structural dynamics is governed by a nonlinear dynamic membrane model. The Navier-Stokes equations are solved using a pressure- based method on a structured and multi-block grid. The membrane material follows the hyperelastic Mooney-Rivlin model [64]. The membrane model is derived from the principle of virtual work. It is spatially discretized with a triangular finite element method and temporally integrated with the Wilson-O method [113]. To accommodate the membrane shape change, the flow solver is further augmented with a moving grid technique discussed in C'! lpter 3. The interaction system is coupled through the exchange of the aerodynamic forces and shape deformations between the fluid and structure solvers. Since the fluid solver and structural solvers do not share the identical grid at the interface, the thin plate spline (TPS) interpolation [16] is used to map the external load and the shape deformation to the structural and fluid solvers, respectively. To avoid phase lag between the fluid and structural solvers, subiterations are performed to synchronize the fluid and structural solvers. In the following we first introduce the interface algorithm between the fluid and structural solvers, this is followed by a step-by-step description of the fluid and structure interaction algorithm. This coupling is achieved through exchanging infor- mation and synchronization between the fluid and structural solvers. In the aerody- namic analysis part, we study both the rigid and membrane wings. In the rigid wing study, we focus on the characteristics of the low Reynolds number and low aspect ratio rigid wing, including boundary 1 -v-r separation and reattachment, tip vortices, and unsteady phenomenon at high angle of attack. In the membrane wing study, we study the self-initiated membrane vibration, its impacts on the membrane wing performance, and we also compare the membrane and rigid wing performances. 6.2 Coupling between the Fluid and Structural Solvers When a fluid flow passes a flexible structure, the structure changes its shape due to the external forces situated in flow; meanwhile, the shape change affects the fluid flow around the structure and the force distributions on the structure. This leads to a fluid and structure interaction problem. Usually, the fluid problem and the structural problem have different mathematical and numerical properties. De- pending on the problem complexity (linear problem or nonlinear problem, simple or complex geometry, small-scale or large scale problem), available software, and com- puting platform, numerous techniques have been developed to provide coupled CFD and CSD methods. One approach is to treat the flow and structural equations as one monolithic system of equations [5] and solve the system of equations simultaneously on a "monolithic" grid. However, since multiple time scales present in a fluid and structural system, a unified treatment may not be efficient to handle the computa- tional stiffness. Liu et al. [58] solved the flow and structural fields simultaneously by a fully implicit method. However, unlike the monolithic method, their CFD solver and CSD solver are solved on different grid systems in which The fluid and structural governing equations are regarded as one single system of time-dependent equations in a pseudo-time. To better handle different characteristics in the flow and structural domains, most of the researchers solved the fluid and structural equations separately and then cou- pled them through a synchronization process. Smith and Shyy [93], Kamakoti et al. [49], Gordnier and Fithen [22], Kalro and Tezduyar [47], Zwaan and Prananta [115], and Lian et al. [53, 55] coupled the fluid and structural solvers with a subiteration approach. The fluid solver and the structural solver function independently on their own computational grids with their own time steps, subiterations between these two solvers are employ. .1 at each time step to avoid the phase lag. Alternatively, Farhat et al. [26, 28] formulated the fluid and structure interaction problem as a three-field problem: the fluid, the structure, and the dynamic mesh that is represented by a pseudo-structural system. In their work, the arbitrary Lagrangian Eulerian (ALE) finite volume scheme is emploi-y, for the fluid dynamics equations, a finite element method is applied for the structural dynamics equations. Partitioned procedures and -I I,.-:- red algorithms are then adopted for the coupled fluid and structure interaction problems in the time domain. In fluid-structure interaction computations, a moving grid technique can be fruit- fully employ. 1 to automatically regenerate the new CFD grid according for the de- formed configuration. In addition, since the CFD solver and CSD solver do not necessarily share identical grid at the interface, interpolations are needed to allow the information exchange. Smith et al. [91] evaluated suitable methods to transfer information between CFD and CSD grids. In our study a thin plate interpolation method [16] is adopted for this purpose. The thin plate spline interpolation is a global interpolation method, it is invariant with respect to rotations and translations, and hence is suitable for the interpolation on moving surfaces. A one-dimensional version of the interpolation is N H(x) = Y aix X2log Ix Xi, (6.1) i= 1 where H(x) is the displacement distribution function over the domain, ai is the unde- termined coefficient, N is the number of monitored structural nodes on the interface, and xi's are their locations. It can be seen from Eq. (6.1) that the interpolated function H(x) has a continuous first order derivative. Once these coefficients ai are determined based on the structural grid locations, the displacement vector defined on the CSD grid D is mapped to the CFD grid via a matrix G 6X = GD. (6.2) where 6X is the corresponding displacement vector on the CFD grid. Similarly the surface forces can also be mapped from the CFD grid to the CSD grid. Upon updating the aerodynamic forces in the structural equations and providing the new boundary conditions to the fluid solver, the temporal development between the aerodynamic and structural equations can be coordinated. To conduct time- accurate computations for the fluid-membrane dynamics, iterations between the fluid and structural solvers should be done at each time instant. A step-by-step coupled algorithm is described below: 1. Generate the CFD and CSD grids. 2. Set n 0, evaluate pressure Po on the undeformed wing based on the initial values po, no, and po. 3. n=n+l, and t,=t,_l + At. 4. Solve the Navier-Stokes equations and transfer the pressure P,+i to the struc- ture based on the TPS interpolation. 5. Use the structural solver to evaluate the displacement vector D,+i. 6. Perform subiterations to synchronize the fluid and structural solvers, and set i 0. (a) i i+1; (b) Map the displacements from the CSD grid to the CSD grid with the TPS interpolation 6X+ = GD+ . (c) Update the CFD grid with the moving grid technique. (d) Update the Jacobian to satisfy the geometric conservation law. (e) Update the boundary conditions for the fluid solver. (f) Compute the pressure PT4+1 field based on the new CFD grid and boundary conditions. (g) Map the aerodynamic load from the CFD grid to the CSD grid with the TPS interpolation. (h) Calculate the structural displacement D'i. (i) C('! 1: for convergence for the subiteration process. If 'No', return to step a), otherwise, continue. 7. Return to step 3 to process next time step. It should be noted that a fully converged solution in the process is hard to achieve. The number of subiteration is based on the computational resource, problem nonlinearity, and the time step. In our study we find two or three subiterations are enough to alleviate the phase lag errors. 6.3 MAV Wing Aerodynamics and Structural Response As shown in Fig. 1-1 the membrane wing consists of a leading edge spar and three chordwise battens on each side of the half wing. The membrane material is bonded to the spar and battens. While the membrane is flexible and does not sustain bending moment, the carbon fiber is rigid which provides support for the membrane and can sustain bending moment. To fully model the membrane wing, ideally, one needs to model both the rigid battens and the flexible membrane. Two principal difficulties arise by doing that. First, the membrane does not bear bending moment, it has three degree of freedoms at each node; the batten, if modeled as a beam, has five degrees of freedom at each node. A special treatment is required to treat the interface point which belongs to different regimes and has different degrees of freedom. Second, the membrane is much thinner than the batten, the mass matrix associated with these two different materials makes the assembled mass matrix ill-conditioned. Given these factors, we treat the batten as a special membrane material with a larger density while the other properties are the same; the density ratio between the batten and the membrane is set to be three. 6.3.1 Computational Background In our study only the membrane wing is considered while ignoring the fuselage and propellers. A schematic geometry of the wing is shown in Fig. 6-1. The shaded area shown corresponds to the fuselage which is fixed. The physical problem under consideration is viscous flow over a wing in an unbounded domain. Based on the freestream velocity of 10 m/s, the root chord Reynolds number is 9x104. Computa- tions are first performed to assess the sensitivity of the computed results to the wing location and the types of boundary conditions imposed to the outer flow boundaries. Basically, the boundaries should be set sufficiently far away from the membrane so that they do not have much impact on the computed aerodynamic coefficients. Based on our tests a proper wing position is chosen and shown in Fig. 6-2. The upstream inlet boundary is 7c from the leading edge, here c is the chord length at the root. At surface ABCD a zero gradient boundary condition is imposed for the velocity compo- nents; all the other boundaries are assigned as a Dirichlet type boundary condition. Since the freestream velocity is parallel to the chordwise direction and no propeller is modelled, only half of the domain is computed, and a symmetric mapping can be applied to the other half domain. Grid sensitivity tests are also performed. Three grid systems have been system- atically chosen, ranging from 1.8 x 105 nodes at the coarse level and 2.3 x 106 nodes at the fine level. For the coarse grid there are 41 points in the chordwise direction and 31 points in the spanwise direction on the wing surface. The grid distribution 60 z SI I Batten 3 Batten 2 Batten 1 Figure 6-1: Schematic wing with three battens on each half wing. (9c, 5c, 6c B Membrane Wing C (-6c, -5c, Figure 6-2: Computational set-up and boundary conditions. at the coarse level around the wing is shown in Fig. 6-3. The results from the fine grid are used as the reference. The comparison among these three grid systems is summarized in Table 6-1. The coarse grid solution under-predicts the lift-to-drag ratio (CL/CD) by 0.5'. and the difference in form drag is less than 0.7'. A similar conclusion can be drawn for the lift, L. As expected, the skin friction drag, DF is sensitive to the grid density; the coarse mesh overestimates DF by almost 10' An intermediate grid, which has 6.7x105 nodes and more than doubles the surface points on the wing surface comparing to the coarse grid, gives an improved prediction on skin friction drag. From Table 6-1 it can be seen that skin friction drag has much less contribution to the total drag than form drag. This is also valid at high angles of attack where the form drag is the dominant factor in the total drag. Figure 6-3: Schematic wing geometry and structured grid with 1.8x 105 points and 41 x 31 on half wing surface. We further compare the main flow features. Flow separations are observed with all these computational grids at the lower wing surface. We extract the chordwise velocity profile at the quarter chord of the root region. The boundary 1 -i-r velocity Table 6-1: Grid refinement effects on the computed aerodynamic coefficients. GRID Dp(N) L(N) DF(N) CL/CD Coarse 1.8x 10 6.24E-2 5.10E-1 9.70E-3 7.06 Intermediate 6.7x105 6.16E-2 5.08E-1 9.22E-3 7.16 Fine 2.3x 106 6.20E-2 5.04E-1 8.80E-3 7.10 profile is shown in Fig. 6-4. The fine grid has the largest reverse velocity, the coarse grid has the smallest, and the intermediate grid has the thickest separation zone. These results indicate that the grid resolution noticeably affects the sharpness of the velocity profile, as well as the extent of the flow separation. Nevertheless, the coarse grid predicts the main features of the flow which is qualitatively consistent with those on the fine grid. -2 Coarse Fine -1 -0.5- o 0 0.5- , 1.5- .... ..... . 12 -5 0 5 10 Spanwise velocity (m/s) Figure 6-4: Spanwise velocity profiles at the bottom surface at the quarter chord of the root region for rigid wing at a= 6. Comparisons between the intermediate and coarse meshes at other angles of attack are also conducted. The computed aerodynamic coefficients are shown in Fig. 6-5. Both the lift and drag coefficients show good agreements, with a little bit larger difference at higher angles of attack. A few more words should be put for the computational results at high angles of attack. Generally -I'" i1i:- vortex shedding and boundary l-v-r separation occur at large angles of attack, and they introduce unsteadiness to the aerodynamic performance and therefore affect the aerodynamic coefficients. Cummings et al. [10] reported that at large angles of attack the unsteady computations predicted much lower lift coefficients than the steady computations. We perform steady computations when the incidence is less than 15 while conducting both steady and unsteady computations when the angle of attack is larger than that. At high angles of attack, truly steady solutions can not be obtained. The -I. "ly" solutions are obtained with steady computations once certain order is reduced in the residual. The PISO method discussed in C'!I pter 2 is used for the unsteady com- putation. We set the time step to be 5x10-5s. However, for our computations, the difference between the steady and the time-averaged lift and drag coefficients, are found to be less than 1 over all the studied cases. Fig. 6-5 shows the computed aerodynamic coefficients for the rigid wing with both coarse and intermediate meshes. At low angles of attack, the coarse grid predicts slightly larger lift than the interme- diate grid; the trend is reversed at high angles of attack. However, both predict that the stall occurs approximately at 51. 2.2 :3 2 2 *, 2.l.r-r Gr,.J 25 0 0 l 05 0.4t 0. 0 0 10 20 30 40 50 60 0 10 20 30 40 50 60 Angle of Attack Angle of Attack Figure 6-5: Aerodynamic coefficients versus the angle of attack: (left) lift coefficient; (right) drag coefficient. For a problem as complicated as the present one, truly grid independent results are difficult to attain. Considering the available computing resource and the out- come of the aforementioned investigation, we use the intermediate grid for the rigid wing simulation while using the coarse grid to illustrate the key issues in the fluid and flexible structure interaction system, to demonstrate the salient features of mem- brane wing computations, and to illustrate the related aerodynamic and structural characteristics. 6.3.2 Rigid Wing Aerodynamics Vortical Flow Structures of the Rigid Wing Tip vortices exist on a finite wing due to the pressure difference between the upper and lower wing surfaces. The tip vortex establishes a circulatory motion that presents over the wing surfaces and exerts great influence on wing aerodynamics. For example, tip vortices increases the drag. The total drag coefficient on a finite wing at subsonic speeds can be written as [2] C2 CD CD,P+CD,F+ L-A (6.3) where CD,P is the form drag coefficient due to the pressure, CD,F is the friction drag coefficient due to viscosity, e is the span efficiency factor which is less than one, AR is the aspect ratio, and Ct CD,i --r i reAR is the induced drag coefficient due to the existence of tip vortices. Eq. (6.3) demon- strates that induced drag varies as the square of the lift coefficient. Furthermore, Eq. (6.3) illustrates that as AR is decreased, induced drag increases. The MAV wing in our study has a low aspect ratio of 1.4; therefore, it is important to investigate the tip vortex effects on the wing aerodynamics. Besides their effects on the di tip vortices have two-folded impacts on the lift: Tip vortex causes a downwash component that decreases the effective angle of attack and therefore decreases the lift [2]. Tip vortex forms a low pressure region on the upper wing surface, and provides additional lift [67]. Fig. 6-6 shows tip vortices around the wing surface at different chordwise po- sitions together with the streamlines at the angle of attack of 39" [53]. The higher pressure from the lower surface drives the flow toward the upper wing surface where the the pressure is lower. The two vortices on both sides of the wing demonstrate clockwise and anti-clockwise rotations; it seems that the flow is leading from the lower surface to the upper surface. Fig. 6-7 shows the low pressure zone due to the vortical flow. The pressure drop further strengthens the swirl by attracting more flu- ids toward the vortex core; meanwhile, the pressure decreases correspondingly in the vortex core. The low pressure region created by the vortex generates additional lift. Toward downstream, the pressure recovers to its ambient value. Then the swirling motion weakens and the diameter of the vortex core increases. The vortex core losses its coherent structure at downstream. Fig. 6-8 visualizes the evolution of the vortical flow structure with the angle of attack. The pressure distribution on the upper surface is also presented in the same figure. At a= 6 tip vortices are clearly visible even though they only cover a small area and are of modest strength. The flow is attached to the upper surface and follows the chordwise direction. A low pressure region near the tip, which is indicated with green color, is observed. Vortices strengthen with the increase of angle of attack. At a 27, as shown in Fig. 6-8, tip vortices develop a strong swirling motion and entrain the surrounding flow to the vortex core. The low pressure area increases as the angle of attack becomes higher. The pressure drop in the vortex core is demonstrated in Fig. 6-9a where the spanwise pressure coefficient on the upper wing surface at x/c=0.4 is plotted. At Figure 6-6: Streamline and vortices for rigid wing at a= 39. The vortical structures are shown on selected planes. o 10060 rp Figure 6-7: Pressure distribution around the rigid wing in the cross sections with streamlines at angle of attack of 39. (c) 270 "-r7\ \ C (d) 330 (e) 450 (f) 510 Figure 6-8: Evolution of flow pattern for rigid wing versus angles of attack. 'A A/ (a) 60 (b) 15" A: a =6 the spanwise pressure is almost uniform on the upper wing surface, and the pressure drop occurs at approximately 91I' of the half span from the root. With the increase of the angle of attack, the pressure drop moves further towards the root. At a 27, it occurs at the 7 .'- from the root. At lower angles of attack, the vortex core position shows a linear relation with the incidence. This relation diminishes at higher angles of attack when the flow is separated on the upper surface. For example, at a=45, the flow is separated at the leading edge, and the low pressure zone covers more than ,!11' of the wing surface, which helps to maintain the increase in lift. At a=51, considerable spanwise velocity component is seen and the flow is separated from most of the upper surface (Fig. 6-8). The lift by the vortices is not enough to maintain the increase of the lift coefficient with the angle of attack, and stall occurs. Even thought tip vortices affect the pressure distribution at the lower surface, the effected region is mainly located near the tip. Most regions in the lower surface are not affected by the tip vortices. From Fig. 6-9b it can be seen that the pressure is almost uniform in the spanwise direction even at high angles of attack. Fig. 6-9 is illustrative in regard to the pressure distributions versus the vortical structures, however, they are not indicative of the total level of the pressure force. One should also note that Fig. 6-9 indicates that the moment experiences substantial variations as the angle of attack changes. The Laminar Boundary Lv,- r Separation on Rigid Wing The laminar boundary lI,--', especially under low Reynolds number conditions, is prone to separate under an adverse pressure gradient. This separation is usually followed by a quick transition from the laminar to turbulent flow. If the adverse pressure gradient is modest, the resulting turbulent flow, by obtaining energy through entrainment, may reattach to the surface to form a laminar separation bubble and then forms an attached turbulent boundary lV,--r. 69 -2 o I | a=51 -e- a=45 =I .-2 -- a=270 -J =J- .' a=150 -25 a=60 S" -3- 0 02 04 06 08 1 350 04 06 08 z/Z z/Z (a) c, at upper surface (b) c, at lower surface Figure 6-9: Spanwise pressure coefficient distributions at x/c=0.4 for rigid wing at different angles of attack. Ever since its first observation by Jones [42], many experimental investigations on the laminar separation bubble have been conducted, as reviewed by Young and Horton [114]. Typically, a laminar separation bubble has the following characteris- tics. Just downstream of the separation point, there is a "dead-air" region, where the recalculating velocity is significantly less than the freestream velocity and the flow is virtually stationary. The separated flow forms a free-shear li-.-r, which is highly unstable. The separated shear 1i--r quickly undergoes transition from laminar to turbulent flow. The turbulent flow may reattach to the surface behind a vortical structure and forms a turbulent boundary l1-, -r. Hence, a separation bubble forms. The dynamics of a laminar separation bubble depends on the Reynolds number, the pressure distribution, the geometry, the surface roughness, and the freestream tur- bulence intensity. A rough rule to determine the reattachment point is given by Carmichael [7] that -zi- the Reynolds number based on the freestream velocity and the distance from the separation point to the reattachment point is approximately 5x 104. On the one hand, if the Reynolds number is less than 5x 104, an airfoil will experience separation without reattachment; on the other hand, a long separation bubble will occur if the Reynolds number is slightly higher than 5x 04. In the fol- lowing, detailed flow structures around a rigid wing are presented. By contrasting the detailed fluid flow around the rigid and membrane wings, one can better evaluate the MAV technologies. The bubble structures at different angles of attack are demon- strated in Fig. 6-10. Representative instantaneous streamlines at the root section on the rigid wing are plotted. At this moment, the spanwise velocity component is ignored to make a clear representation of these separated structures. At a low angle of attack of 6, the flow primarily attaches to the surface and a tiny separation bubble is observed on the lower wing surface near the leading edge. Beginning at 45'. chord from the leading edge, a weak recirculation zone is seen on the upper wing surface, which produces a reattachment length of x/c=0.5. With the increase of incidence, the separation point moves forward toward the leading edge. At the angle of attack of 150, as shown in Fig. 6-10, the separation oc- curs at I'. of the chord from the leading edge, which is followed by a long "dead-air" zone before it reaches a maximal reverse velocity of 0.26U. In the "dead-air" zone, the stationary reverse flow does not have much effect on the pressure distribution, which is primary determined by the wing curvature. The shear 1-~.-r reattachment happens at I' of the chord. Based on the freestream velocity and the distance between the separation and reattachment points, the Reynolds number is approximately 4.5x 104 which is slightly smaller than the Reynolds number sI-:-, -1 .1 by Carmichael. The reattachment point corresponds to a rapid increase in the surface pressure, it is highly unstable, and moves forward and backward. The streamwise velocity contours at different angles of attack are demonstrated in Fig. 6-11. The velocity is normalized by the freestream velocity. At a= 6, the maximum reverse velocity is less than 0.5'. of the freestream velocity. Under such a condition, experiment shows that a considerable portion of the shear 1-,-r remains laminar [21]. (a) 60 --- -------- (b) 150 (c) 270 (d) 390 (e) 510 Figure 6-10: Streamlines at different angles of attack for rigid wing. Shown are single span shorts of the individual time dependent flows. However, with the increase of angle of attack, the magnitude of the reverse velocity increases as well. At a=27', the reverse flow component of the separation bubble reaches a maximal velocity of 0.49U. Under such a situation, Crompton and Barrett [9] showed that the shear 1~ ,-r is particularly energetic since most of the shear 1lv.,r is turbulent. Even though the flow on the upper surface near the root has separated from a large portion of the surface; however, the lift still increases with the angle of attack. Two reasons might lie behind this. First, tip vortices generate suction areas on the upper wing surface which provide additional lift. Second, even though flow separates near the root, it still attaches to the other region of the wing surface. When the angle of attack is increased further, at a=51", massive separation appears, and stall occurs. Meanwhile, as seen from Fig. 6-12, vortex shedding is observed near the root. A tiny separation bubble emerges at the leading edge first, it then obtains energy from the shear l1v. r and increases its size. This bubble then merges with another bubble downstream to form a larger bubble. The larger bubble is not stable, it breaks into two small bubbles near the maximal camber position, and vortex shedding begins. For low aspect ratio wings, tip vortices have considerable contributions to the lift. This scenario is quite similar to that of delta wings [23]. We have observed that the low aspect ratio wing suffers less from the separation. As seen from Fig. 6-5a, the lift keeps a linear relationship with the angle of attack even at very high angles of attack. Experiments by Torres and Mueller [97] on low aspect ratio wings have similar findings. (a) 6 ^\ 1.30 0.91 .652 0.13/ -0.26 (b) 150 \ 1. 3 1.36 S0.6 -0.7 I 2.40 0.50 ,, /f-"---------Y- I-OO (d) 39 1.45 0.50 [] 0.45 1.40 (e) 51 Figure 6-11: Normalized chordwise velocity u/U contours at different angles of attack. Figure 612: Vortex shedding at a51 Figure 6-12: Vortex shedding at a=51. Unsteady Phenomenon at High Angles of Attack Both vortex shedding and boundary 1-v-r separation occur at large angles of attack, and they introduce unsteadiness to the aerodynamic performance. As men- tioned before, we perform steady computations when the incidence is less than 15 while conducting both steady and unsteady computations when the angle of attack is larger than that. Interestingly, the difference between the steady state computa- tions and the time-averaged unsteady computations are small even at large angles of attack in which unsteady phenomenon such as vortex shedding are prominent. The differences in the lift coefficient and the lift-to-drag ratio are found to be less than 1 (Fig. 6-13). Noticeable difference appears only when the angle of attack becomes sufficiently high. 2.5 -0- Time Averaged -- Steady 2-.5 0 10 20 30 40 50 60 Angle of Attack (degree) Figure 6-13: Comparison of the rigid wing CL between the steady and unsteady computations. The pressure distributions are also compared at the root section, where flow separation usually appears first due to the large camber in the present MAV wing design [34]. Fig. 6-14 shows that at a= 6 the time averaged pressure coefficient 0-1 -08 -0 -1 -04- -1 -06 12 0 02 04 06 08 1 0 02 04 06 08 1 x/c x/c Figure 6-14: The cp comparison on the rigid wing at the root between the steady and unsteady computations: (left) a=6; (right) a=15. matches closely with the steady state result. This finding is consistent with Fig. 6 10 and Fig. 6-11 where a very weak recirculation is seen at this angle of attack. However, at a=15, clear difference is shown after x/c=0.6 which approximately corresponds to the location of the maximal reverse velocity in Fig. 6-11. The time averaged value yields a smooth pressure distribution; the variation in the steady result is apparent from the recirculation zone. Velocity distributions show the same trend as the pressure distributions. At a= 6, as seen from Fig. 6-15, there is almost no difference in the chordwise velocity distribution. However, at a=15, no difference near the leading edge where no separation occurs, but clear difference is shown after the separation bubble (Fig. 6-16). 6.3.3 Membrane Wing Dynamics Time Synchronization Before presenting the membrane wing results, some key issues in coupled membrane- fluid interactions are first addressed. Specifically, the time synchronization between the fluid and structural solutions, and the geometric conservation law in regard to the moving grid method are emphasized. Boundary layer profile at x/c=0 14 for a=60 Time Averaged s-- Steady 0021 0 0205 S002 00195 0019 00185 0 02 04 06 08 u/U 1 12 14 Figure 6-15: Comparison of chordwise velocity profiles on the rigid wing between the steady and time-averaged unsteady computations at the angle of attack of 6. Velocity profiles are at two chordwise locations, and both at the root section. Boundary layer profile atx/c=0 19 fora=150 Time Averaged Steady 1 12 14 Figure 6-16: Comparison of the chordwise velocity profiles on the rigid wing be- tween the steady and time-averaged unsteady computations at angle of attack of 15. Velocity profiles at at two chordwise locations, and both at the root section. 00215 Boundary layer profile at x/c=0 80 for a=60 00215 0021 0 0205 > 002 00195 Boundary layer profile at x/c=0 80 for a=150 0019 00185 0 02 04 06 08 u/U The time synchronization in our study means that iterations should be enforced between the fluid and structural solvers at each time step to avoid phase lag errors. We shall simulate the fluid and structure interaction by integrating the fluid and the structural solvers with subiteration. To demonstrate the importance of the time synchronization between the fluid and structural solvers, We compare the results with and without synchronization in Fig. 6-17 that depicts the displacement history of a trailing edge node. The computation without synchronization fails to continue while that with synchronization keeps on going. The displacement history without synchronization matches well with that with synchronization at the very early stage. However, the 1 ,-.ii-.; errors accumulate with time and eventually result in a much larger displacement than that with synchronization. Gordnier and Visbal [23] have also reported the importance of the time synchronization. The velocity and pressure histories of the same node with and without time synchronization are compared in Fig. 6-18 and Fig. 6-19, respectively. Before the computation diverges, the velocity without synchronization is more than two times larger than that with synchronization. The larger velocity of the membrane node means that the membrane node obtaining more energy from the surrounding flow than that with synchronization. The high velocity is believed to be associated with the sudden pressure drop shown in Fig. 6 19, which causes the wrinkle of the membrane and the failure of the structural solver. The Geometric Conservation Law The geometric conservation law [100] is another important factor when the mov- ing grid technique is employ, l Fig. 6-20a shows the time history of CL/CD for the membrane wing at a= 6. Without satisfying the geometric conservation law, irreg- ular spikes are observed in the course of computations. However, the computations are better regulated once the geometric conservation law is enforced as in Fig. 6-20b. Farhat et al. [27] argued that the geometric conservation law was a necessary con- dition to maintain the stability of a scheme utilized in moving boundary problems. With Synchronization Without Synchronization 0.1 Time (Second) 0.15 Figure 6-17: Effects of time synchronization on the displacement point for membrane wing at a=6 . of a trailing edge With Synchronization Without Synchronization -30 0 0.05 0.1 0.15 Time (Second) Figure 6-18: Effects of time synchronization on vertical velocity point for membrane wing at a= 6. 0.2 of a trailing edge 0.12 I 0.1 I 0.08 - 0.06 0.04 0.02 -0.02L 0 1 80 0 , -1 With Synchronization -2 Without Synchronization -3- -4 -5 -6 -7 -8 0 0.05 0.1 0.15 0.2 Time (Second) Figure 6-19: Effects of time synchronization on pressure distribution on a trailing edge point for membrane wing at a=6. We do not come across instability in our study. Two possibilities lie behind that: one is the small time step we use for the fluid solver, the another is the relatively small displacement in the membrane wing. Nevertheless, the computations behave erratically when the geometric conservation law is not enforced. Self-initiated Membrane Vibration Earlier works in membrane wings [88, 92] were based on considerations that the membrane was massless and there was no time dependent movement in the steady freestream. Experiments [111] showed that the membrane experienced high frequency vibration. By adopting a dynamic membrane model, our computations show that the membrane wing also exerts high frequency vibrations in the steady freestream. Fig. 6-21 demonstrates the displacement history of the trailing edge with time. The maximal displacement during that period occurs near the tip. The displacement is normalized by the maximal camber that is about 0.9 cm. To investigate the vibration frequency, we choose a point between the batten 1 and batten 2 on the trailing edge 5 75 i_! i i. s, , 0 55 0 45 -5r''' 41'' 0 0.02 0.04 0.0B 0.08 D.1 0.12 0 005 01 015 02 Time (sec) Time (sec) Figure 6-20: Effect of the geometric conservation law on CL/CD for membrane wing computation at a= 6: (left) not satisfying the geometric conservation law; (right) satisfying the geometric conservation law as an example. Its deflection history is shown in Fig. 6-17, and the estimated dom- inated frequency is around 120 Hz (Fig. 6-22). Analyses at other points also show the frequencies around 120 Hz that is comparable to the experimental frequency of 140 Hz. [111] Considering that the experimental conditions and the detailed wing configuration between our study and those reported by Waszak et al. [111] are not identical, the qualitative agreement between computations and experiments in this regard is deemed satisfactory. Liu [60] reported that under a typical wing gust situa- tion the energy was mainly located in the low frequency range of several Hz or lower. The membrane fluctuation is not expected to cause sensitive response to the vehicle since the membrane fluctuates in a time scale much faster than either the vehicle control scale or the expected wing gust time scale. When the membrane surface vibrates, the velocity on the wing surface is no longer zero. A considerable velocity component is observed at the wing surface. The vertical velocity contours at two different time instants are plotted in Fig. 6-23. Fig. 6-23a corresponds to a stage when the membrane wing moves up under the aerodynamic forces. At the trailing edge, the vertical velocity is about S'. of the 0.14 0.12 0. 1 0 . E 0.0c CL 0.0-4 N S0.02 z t=o t=3 x 10-3 t=7 v 10-3 1=10 10- I 1,. 10- l=1. 10 -0.02[ -0.04 0 0.2 Figure 6-21: Time history of trailing The camber at the root is 0.90 cm. 0.4 0.6 0.8 edge displacement for membrane wing at a= 6. -^3 x 10-3 7 6- _5- 4- 2- 50 100 150 200 250 Frequency (Hz) 300 350 400 Figure 6-22: Spectrum analysis of the trailing edge point vibration for membrane wing at a=6. 83 v/U contour at t=3.0 x 10-3 v/U contour at t=7.0 x 10-3 0 0 ; 0746 -0.000772 2 2 -0.000772 -,iIJw, 1 JGP 0. 00438 6 r 4 -4o.000772 J} 4 S0.0115 S/ 6 r, ,,, -0.000772 V \I I, I / -2 02 210 12-0 1 '0 12- 12 00592 -6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6 Spanwise (cm) Spanwise (cm) Figure 6-23: Normalized vertical velocity contour at two different time instants for membrane wing at a= 6 on the upper surface. freestream velocity. A negative velocity component is also observed near the leading edge. Fig. 6-23b approximately corresponds to a stage of maximal displacement. A negative velocity near the root indicates that the membrane begins to move down. Comparison between the Membrane and Rigid Wings With the same freestream condition, as shown in Table 6-2, the flexible wing exhibits slightly less lift coefficient than the rigid one at a= 6. The difference in CL/CD is also small. At higher angle of attack of 15, the membrane wing generates a lift coefficient about 2'. less than the rigid wing does; however, its CL/CD is 1.5'. more than the rigid one. These differences lie behind the high frequency of the membrane vibration. Under the external force, the membrane wing changes its shape. This shape change has two effects. On the one hand, it decreases the lift by reducing the effective angle of attack of the membrane wing; on the other hand, it increases the lift by increasing the camber. Our numerical findings show that membrane and rigid wings exhibit comparable aerodynamic performance before stall limit, which was also experimentally observed by Waszak et al. [111]. Fig. 6-24 shows the time averaged vertical displacements of the trailing edge at a=6 and a=15. The displacements are normalized by the maximal camber of the Table 6-2: Aerodynamic comparison between membrane and rigid wings. Membrane Wing Rigid Wing 6 CL/CD 7.05 7.06 CL 0.530 0.532 15 CL/CD 3.94 3.88 CL 0.920 0.954 wing. At a =6 the deflection is about 15'. of the maximal camber and increases to more than 211'. at a=15. Due to the trailing edge deflection the effective angle of attack of the membrane wing is less than that of the rigid wing. The spanwise angles of attack between rigid and membrane wings under the same flow condition and with identical initial geometric configurations are shown in Fig. 6-25. As we can see from Fig. 6-25a, the rigid wing has an incidence of 6" at the root and it monotonically increases to 9.5" at the tip; the membrane wing shares the same angles of attack with the rigid wing in the 3.' of the inner wing; however, the effective angle of attack toward the tip is less than that of the rigid wing. At the tip, the angle of attack of the membrane wing lowers by about 0.8". Fig. 6-25b compares the angle of attack at a=15", and it shows that the effective angle of attack of the membrane wing is about 1 less than that of the rigid wing at the tip. The reduced effective angle of attack causes the decrease in the lift. The membrane wing shape change has another effect. The pressure difference between the upper and lower wing surfaces inflates the membrane wing's camber, which is shown in Fig. 6-26. Two airfoil shapes at different spanwise positions are plotted together with the corresponding rigid wing shapes. The camber increase is more visible in the outer wing than that in the inner wing. This is consistent with Fig. 6-24 where the maximal displacement occurs near the tip. As expected, the increase is larger at a=15" than that at a=6. The wing surface movement affects the overall pressure distribution. Fig. 6 27 compares the time averaged pressure contours between the rigid and membrane 015 E S01 - 8 o Qi 5-_ o Rigid part 0 02 04 w 02 _o 015l S01 E N 0 05 :z Rigid part 06 08 1 0 02 04 06 08 Figure 6-24: Averaged displacement of the membrane wing trailing edge: (left) a=6; (right) a=15. Rigid Wing 1ng Membrane Wing Membrane Wng 95 185 9 18 65 55 6 15 0 02 04 06 08 1 0 02 04 06 08 1 Spanwise z/Z Spanwise z/Z Figure 6-25: Comparison of the spanwise angles of attack between the rigid and membrane wings: (left) a=60; (right) a=15. wings. As expected, the differences are mostly located at the trailing edge where large movement is observed. These differences are also reflected in the chordwise pressure distributions shown in Fig. 6-28. At the root section there is almost no difference between the rigid and membrane wings since the fuselage is fixed. How- ever, the difference becomes clear toward the tip; at z/Z=0.37, the difference is visible. At z/Z=0.69, the time averaged membrane wing shift the lowest pressure point downstream comparing with the rigid wing result. This shift helps to delay flow separation. 025r Spanwise zZ=O 40 a=6 a=15 0 02 04 06 08 1 Figure 6-26: Time averaged airfoil shapes at membrane wing at a= 6 and 15. Pressure Coefficient Contour -0.299 -,n -0.5 -0.567 f0.634 -0.567 -0.433 -0.567 -0. -1.2-, -0.433 Spanwise (cm) Figure 6-27: Comparison of c, contours at a=6: wing; (right) steady rigid wing. (left) time averaged membrane 02 04 06 08 different spanwise positions for the 4-/ -1.1 Spanwise z/Z=O 60 \ Time Averaged Membrane Rigid Wing o- 0 -0.8- 0 0.2 0.4 0.6 0.8 1 x/c Time Averaged Membrane Rigid Wing -0.2 -0.2 0.2 0.4 0.6 0.8 x/c 0.1 0 -0.1 Time Averaged Membrane Rigid Wing 0 0.2 0.4 0.6 0.8 1 x/c Figure 6-28: C'! i Il -. c, comparison at different spanwise locations between time averaged membrane wing and steady rigid wing at a=6: (a) z/Z 0; (b) z/Z=0.37; (c) z/Z 0.69. CHAPTER 7 WING SHAPE OPTIMIZATION In C'! ipter 6 we have performed the aerodynamic study of the wing aerodynamics of both the rigid and flexible wings. In this chapter we will develop a systematic approach to enhance the membrane wing performance. 7.1 Introduction It is desirable not only to understand the pertinent physical phenomena of the membrane wing dynamics, but to improve its performance based on the understand- ing. In C'! lpter 6 we have investigated both rigid and membrane wing aerodynamics. The wing shape adopted is from an existing design. One question then arises natu- rally: is that the best wing design, or, can we improve its performance. To answer this question, it is worthwhile to pause and take a look at how the existing wing is de- signed. The design of the wing is mostly based on observations and a two-dimensional analysis tool to test the performance of different airfoils. Based on the tests, one air- foil shape is then chosen to form a three-dimensional wing. The wing shape is then tested with flight. In this trial and error process, the three-dimensional effects and the membrane material properties are not considered. A systematic approach is necessary to improve the wing performance. There are many existing mature techniques for shape optimization. Could we simply use one of them to achieve our goal? The answer is negative. As discussed in Chapter 6, the MAV wing performs under the low Reynolds number condition, and it exhibits a high frequency vibration even under steady freestream conditions. Therefore, the optimization of a membrane wing introduces multiple challenges. First, coupled, time dependent simulations of the interaction between viscous fluid flow and flexible structure are very expensive. Second, an efficient and automatic grid regeneration is essential. To deal with these difficulties, we propose the following strategies: Adopt the rigid wing as a surrogate model for membrane wing optimization. The surrogate approach was emploil,, in [6, 32, 51, 63], in our study the sur- rogate model is the rigid wing. We shall optimize the rigid wing under a given flow condition, and evaluate the performance of a flexible wing with output from the surrogate model. Automate a moving grid technique by treating the optimization process as a iu1' img boundary piI, i, This technique facilitates not only the optimiza- tion process but the computation of fluid and membrane wing interactions. In the following, we first introduce the optimization problem. This is followed by an optimization algorithm in which we show how to integrate the analysis tool, the moving grid technique, and the optimizer Design Optimization Tools (DOT) [12] into a system to do shape optimization. The optimization process begins with a rigid wing surrogate optimization, then the flexible wing performance is assessed with the optimized shape. Optimization results and the aerodynamic characteristics of both optimized rigid and membrane wings are then highlighted. 7.2 Optimization and Computational Framework The objective is to improve the lift-to-drag ratio of a membrane wing. The baseline membrane wing shape selected is an existing design [33, 34] devised from solutions provided by XFOIL [15] and empirical modification based on flight tests [34]. XFOIL provides two-dimensional fluid flow solutions based on coupled thin l-iv r and inviscid flow models. A schematic baseline wing shape is shown in Fig. 6-1 in C'! lpter 6. The baseline shape has three carbon-fiber-made battens that are labelled also in Fig. 6-1. The shading area corresponds to the fuselage. The fuselage is made of rigid carbon fiber prepreg cloth which does not deform during the flight. A flexible latex membrane covers the top surface which changes its shape under aerodynamic forces. The baseline wing has a variable spanwise camber. The nominal camber of the baseline design is 7.5' at the root, and 2' at the tip. The wing has a maximal chord length of 13.7 cm, a mean chord of 10.5 cm, and a span of 15.2 cm. C is used to indicate the chord length at different span stations, and Z for the half span. The angle of attack is defined in reference to the root geometry. At the design point, the MAV flies with a speed of 10 m/s and an angle of attack of 6. The objective is to maximize the lift-to-drag ratio at such a flight condition. At such a condition the root chord Reynolds number is 9x 104. Six design variables are chosen, they are the y-coordinates at six points on the surface. Three are located at the root whose positions are approximately at 1(1' 7'. ^ and 77'. of the chord, as seen in Fig. 7-1. The other three are located at Batten 2 with the same relative chord positions. The wing surface can be implicitly represented by a function y = f(x, z), where x, y, and z represent the wing surface coordinates. We indicate the design variable with a capital letter Yi whose value is the perturbation of the y coordinate at point i with respect to its baseline value. We interpolate the perturbations from the six design variables over the wing surface with the thin plate spline (TPS) interpolation method discussed in C'!i lter 6. Therefore, the interpolated perturbations over the wing surface is 6y(x,z) G(Y), (7.1) where Y is a vector indicating the design variable values, G is the interpolation operator, and by is the interpolated value over the wing surface. The perturbed wing surface is therefore represented by the summation of the baseline and the interpolated value in Eq. (7.1), ynew baseline + 6y. One advantage of this perturbation approach is that we can recover the original shape when the perturbation vector Y is zero. In 91 (x2,y2) 1.5- (xO,yO) (x3,y3) 0.5- (a) 0- -0.5 0 0.2 0.4 0.6 0.8 1 x/C 1.5- Design point 3 1-) 0.5- (b) Baseline Shape --- Perturbd Shape 0- -0.5 0 0.2 0.4 0.6 0.8 1 x/C Figure 7-1: A cross section of membrane wing along a batten: (a) straight lines indicating how convexity constraint is applied; (b) effect of y coordinate perturbation at a point. The design points are located at 1C('. 27'. and 77' of the chord. our optimization process, both the leading and trailing edges are fixed to maintain the identical angle of attack between the baseline shape and the optimized shape. The optimization problem investigated can be formulated as follows. Maximize CL/CD Subject to 1: CL > CLbaseline (7.2) Yi + y o yo 0 Y2 +y2 2: Convexity constraint: > X1 Xo X2 X0 3:YL Constraint 1 requires that the lift coefficient be no less than that of the baseline wing. Constraint 2 maintains the convexity of the airfoil (see Fig. 7-la) so to eliminate ob- viously inappropriate shapes. Only one of the six constraints is shown for illustration purpose. Constraint 3 gives the lower and upper bounds of the design variables whose values are listed in Table 7-1. Table 7-1: Design variables bounds and their optimal values at 6. Initial Lower Upper Case 1: opti- Case 2: optimal Parameter design bound bound mal values with values with 1000 (cm) (cm) (cm) 600 iterations per iterations per analysis analysis X1 0.00 -0.40 0.00 -0.258 -0.259 X2 0.00 -0.40 0.00 -0.400 -0.400 X3 0.00 -0.20 0.20 -0.152 -0.139 X4 0.00 -0.20 0.20 -0.100 -0.095 Xs 0.00 -0.10 0.20 -0.082 -0.068 X6 0.00 -0.10 0.20 0.135 0.121 Analysis 97 108 Cycle 6 6 CL 0.530 0.529 0.529 0.529 CL/CD 7.06 7.55 7.55 Number of analyses with convexity violations 13 22 |

Full Text |

PAGE 4 First,Iwouldliketoexpressmysinceregratitudetomyadvisor,Dr.WeiShyy.Heprovidedagoodstudyenvironmentandoeredplentyofopportunitiesformetoexploreandpursuemyresearchinterest.Iwouldliketothankhimforhisenduringenthusiasmineducatingmeasaresearcherandaperson.IhopethatIemulatehisneexampleofasuccessfulinternationalstudent,anintelligentandindustriousscholar,aresponsibleandknowledgeableadvisor,andanoptimisticandhumorousperson.Iowethankstoseveralcolleaguesfortheirwillingnesstosharetheirexperiencesandknowledgewithme.Ihavebenetedmuchfromtheircollaborationinvariousareas:Dr.PeterG.Ifju,microairvehicle;Dr.AndrewKurdila,structuraldynamics;andDr.RaphaelT.Haftka,optimization.Dr.SiddharthThakurgenerouslysharedhistimeandexperience,whichmademyworkmucheasier.Ienjoyedthegoodrelationshipswithmyprofessors.Dr.RenweiMeigavememanyhelpfulsuggestionsandimprovedmyknowledgeofadvanceduidmechanics.Dr.MarkSheplakspentmuchtimegatheringclassmaterialsforme.BoththeUniversityofFloridaandtheMechanicalandAerospaceEngineeringDepartmentprovidedapleasantenvironmenttostudyandlive.Iamproudtohavebeenhere;Ihopethattheywillbeproudtohavemehere.MyfamilymembersinChinahavegivenmetremendoussupportformystudyabroad.Theirtrustandlovehavesupportedmeovertheyearsandwillsupportmeinthefuture.Mywife,LihuiBai,hasbeenwithmeandsupportedmeallthetime.Herpatience,understanding,andencouragementaretheinvaluablewealthofmylife.ButnoacknowledgementcouldpossiblystateallthatIowetoher. iv PAGE 5 page ACKNOWLEDGMENTS ............................. iv ABSTRACT .................................... vii CHAPTER 1INTRODUCTION .............................. 1 1.1MAVandMembraneWingAerodynamics ............. 1 1.2MotivationsandChallenges ..................... 4 1.3Outline ................................. 5 2THEORETICALANDCOMPUTATIONALMODELSFORFLUIDFLOWS .................................... 8 2.1GoverningEquationsandFluidSolver ............... 9 2.1.1GoverningEquationsforFluidFlow ............. 9 2.1.2FluidFlowSolvers ....................... 12 2.1.3GeometricConservationLaw ................. 12 2.2TurbulenceModels .......................... 13 2.2.1ReynoldsEquationsandk-"TurbulenceModel ....... 13 2.2.2BoundaryConditionsforTurbulenceModels ........ 15 2.3Laminar-turbulentTransition .................... 18 3MOVINGGRIDTECHNIQUE ....................... 19 3.1PerturbationMethod ......................... 20 3.2Master/SlaveConcepts ........................ 21 3.3NumericalTests ............................ 26 4ACTIVEFLOWCONTROLWITHANADAPTIVELY-SHAPEDWING 27 4.1Introduction .............................. 27 4.2ExperimentalandComputationalSet-up .............. 28 4.3NumericalStudy ........................... 29 4.3.1GridSensitivityAnalysis ................... 29 4.3.2AssessmentoftheTurbulenceModel ............. 30 4.3.3EectsofGapsbetweenFlaps ................ 30 4.3.4EectsofFlappingAmplitude ................ 32 4.4Conclusions .............................. 34 v PAGE 6 .................. 36 5.1Introduction .............................. 36 5.2StructuralSolver ........................... 38 5.2.1Mooney-RivlinModel ..................... 38 5.2.2PrincipleofVirtualWork ................... 39 5.2.3TriangularMembraneElement ................ 40 5.2.4Green-LagrangeStrainTensor ................ 42 5.2.5InternalForces ......................... 43 5.2.6ExternalForces ........................ 45 5.2.7DynamicEquationsfortheMembrane ............ 46 5.2.8SolutionoftheDynamicEquations ............. 47 5.3NumericalTestofMembraneModel ................. 48 6MEMBRANEWINGAERODYNAMICSFORMAVS .......... 53 6.1Introduction .............................. 53 6.2CouplingbetweentheFluidandStructuralSolvers ........ 55 6.3MAVWingAerodynamicsandStructuralResponse ........ 58 6.3.1ComputationalBackground .................. 59 6.3.2RigidWingAerodynamics .................. 64 6.3.3MembraneWingDynamics .................. 76 7WINGSHAPEOPTIMIZATION ...................... 88 7.1Introduction .............................. 88 7.2OptimizationandComputationalFramework ............ 89 7.3ResultsandDiscussion ........................ 93 7.3.1SensitivityAnalysis ...................... 93 7.3.2PerformanceoftheOptimizedRigidWing ......... 95 7.3.3FlexibleWingOptimization ................. 97 7.4SummaryandConclusions ...................... 102 8CONCLUSIONSANDFUTUREWORK .................. 104 APPENDIX AAPPENDIXSTRAIN-DISPLACEMENTMATRIX ............ 107 REFERENCES ................................... 109 BIOGRAPHICALSKETCH ............................ 118 vi PAGE 7 Microairvehicles(MAVs),withwingspanof15cmorlessandightspeedaround10m/s,havemanyapplicationsinbothcivilianandmilitaryareas.TheReynoldsnumberbasedonthegivenparametersisaround104,whichoftenyieldsinsucientlift-to-dragratio.Furthermore,oneexpectstheunsteadyeecttobenoticeableforsuchightvehicles.Theexiblewinghasbeendemonstratedtoexhibitfavorablecharacteristicssuchaspassiveadaptationtotheightenvironmentanddelayedstall. ThepresentstudyfocusesondevelopingcomputationalandmodelingcapabilitiestobetterunderstandtheMAVaerodynamics.Bothexiblewings,utilizingmem-branematerials,andadaptively-shapedwings,utilizingpiezo-actuatedaps,havebeenstudied.Intheadaptively-shapedwingstudy,weusepiezo-actuatedapstoactivelycontroltheow.Weassesstheimpactsoftheapgeometry,appingam-plitude,andturbulencemodelingontheowstructurewithaparallelexperimentaleort.Themembranewingusesapassivecontrolmechanismtodelaythestallangleandtoprovideasmootherightplatform.Ourstudyfocusesonthemutualinterac-tionsbetweenthemembranewinganditssurroundingviscousow.Wecomparethe vii PAGE 8 Besidestheaerodynamicstudy,wealsoperformshapeoptimizationtoimprovethemembranewingperformance.Sincedirectoptimizationofamembranewingistootimeconsumingtobepractical,weoptimizeasurrogaterigidwingmodelbasedonanintegratedoptimizationalgorithm,whichconsistsofaNavier-Stokessolver,anautomaticgridgenerationtool,andagradient-basedoptimizer.Then,weassessthemembranewingperformancebasedontheoutcomefromthesurrogatemodel.Ournumericalresultsconrmthatthemembranewingexhibitsconsistentimprovementinthelift-to-dragratiowiththesurrogatemodel. viii PAGE 9 1.1 MAVandMembraneWingAerodynamics Microairvehicles(MAVs)withamaximaldimensionof15cmorlessandaightspeedaround10m/sareofinteresttobothmilitaryandcivilianapplications.Equippedwithavideocameraorasensor,thesevehiclescanperformsurveillance,reconnaissance,targeting,andbiochemicalsensingataremoteorotherwisehazardouslocation[ 83 ].Withtherapidprogressmadeinstructuralandmaterialtechnologies,miniaturizationofpowerplants,communication,visualization,andcontroldevices,severalgroupshavedevelopedsuccessfulMAVs[ 24 34 ].Asitssizeisreduced,MAVgainsfavorablescalingcharacteristics,includinglowinertiaandreducedstallingspeed[ 83 ].However,intermsofaerodynamics,control,range,andmaneuverability,manyoutstandingissuesarestillpresent.First,anMAViesinalowReynoldsnumberregime(104105)duetoitslowightspeedandsmalldimension.Suchaightenvironmentisoftenaccompaniedbylaminarboundarylayerseparation,transition,andlowlift-to-dragratio.Second,anMAVwingtypicallyhasalowaspectratio,whichcausesstrongvorticalowstructuresandincreasestheinduceddrag.Third,anMAVissusceptibletorollinginstabilities,whichbecomesevenmoreseriousbecauseoftheexistenceoftipvortices.Fourth,anMAVissensitivetowinguctuations,whichcanbecomparabletoMAV'sightspeedandmakesboththeinstantaneousightReynoldsnumberandangleofattackvarysubstantially[ 83 ]. InthedevelopmentofMAVs,twoapproachesexistintheMAVwingdesign:appingwingsandxedwings.Intheappingwingarea,pioneeringworkwasdonebyWeis-Fogh[ 112 ],andLighthill[ 57 ].RecentexperimentalandnumericalstudiesinthisaspectweregivenbyChasmanandChakravarthy[ 8 ],DeLaurier[ 11 ],Smith 1 PAGE 10 [ 90 ],Dickensonetal.[ 13 ],Ellington[ 17 ],JonesandPlatzer[ 43 44 45 ],Katz[ 50 ],Kamakotietal.[ 48 ],LiuandKawachi[ 59 ],VestandKatz[ 106 107 ],Walker[ 109 ],andWang[ 110 ].Shyyetal.[ 83 ]gaveanextensivereviewofappingandxedwingcharacteristics.Templin[ 96 ]presentedtheappingwingspectrumofightanimals. Ourstudyfocusesonaxedwing.Specically,wewillstudyaexiblemembranewingforMAVapplications.Figure 1{1 showsa15cmMAVwithaexiblewingdesignedbyIfjuandcoworkers[ 34 ].TheMAVhasaightspeedbetween20and40km/h.Poweredbyanelectricmotor,itcanyforabout15minutes.Withabuild-invideocameraandatransmitter,thevehiclehasatotalmassaround40g.TheMAVhasarigidfuselageandleadingedgesparwhicharemadeoflaminatematerial.Oneachsideofthewingtherearethreeunidirectionalcarbonberbattens.Thesparandbattensarecoveredwithalatexmembranematerial.Overall,thewinghasaspanof15.2cm,arootchordof13.2cm,ameanchordof10.5cm,andawingareaof160cm2,whichresultinalowaspectratioof1.4.Thewinghasavariablecamberdecreasingfrom7.5%attherootto2%atthetip.Inthisconguration,themaximumcamberislocatedatabout27%chordattheroot. Figure1{1: MAVwithmembranewing. PAGE 11 Flexiblewingshavedierentaerodynamicperformancefromrigidwings.Oneadvantageofexiblemembranewingsisthattheyfacilitatepassiveshapeadaptationandresultinadelayedstall[ 83 111 ].Fig. 1{2 ,adoptedfromWaszaketal.[ 111 ],comparestheliftcurvesversusanglesofattackforrigidandmembranewingswithcongurationssimilartothoseshowninFigure 1{1 .Undermodestanglesofattack,therigidandmembranewingsdemonstrateasimilarliftcurveslopeof2.9,withtherigidwinghavingaslightlyhigherliftcoecient.However,themembranewingshaveamuchhigherstallanglesofattackthantherigidwing,whichisakeyelementinenhancingthestabilityandagilityofMAVs.Therigidstallsat24owhiletheexiblewingshavestallanglesbetween30oand49o,whicharesimilartothatofmuchloweraspectratiorigidwings(AR=0.5to1.0).But,theverylowaspectwingsexhibitlowerliftcurveslopesof1.3to1.7insteadof2.9.Hence,exiblewingsappeartocombinethedesirableperformanceofrigidwingswithhigherandloweraspectratiosbyexhibitingstallbehaviorsimilartothatofrigidwingswithaspectratioof0.5to1.0andtheliftgeneratingcapabilitysimilartorigidwingswithaspectratioof2.0[ 111 ].Anotheradvantageofexiblemembranewingsisthattheycanadapttothewindgustandprovidesmootherightplatforms,whichareveriedbywindtunnelexperiments[ 84 ]. Figure1{2: Liftcoecientversusangleofattack.(FromWaszaketal.[ 111 ]) PAGE 12 1.2 MotivationsandChallenges Thoughexperimentshaveprovidedinformationonthemembranewingaero-dynamics,tofullyunderstandthemembraneaerodynamicsandtoappreciatethemembranewing'spassivecontrolmechanism,weneedtoperformadetailednumer-icalstudy.Anumericalstudyofthemembranewingbringstwochallenges:thelowReynoldsnumberandthelowaspectratiocondition,andtheinteractionbetweenthemembranewinganditssurroundingviscousows. First,thelowReynoldsnumberconditionpresentstremendouschallengesonMAVsstudy.Inourstudy,thechordReynoldsnumberoftheMAVisaround9104.IntherangeofReynoldsnumberof104106,complexowphenomenaoftenoccur.Thelaminarboundarylayerseparation,transition,andreattachmentusuallycoexistontheupperwingsurface.Anaccuratepredicationofthetransitionfromlaminartoturbulentstateisimportanttoevaluatethewingperformance.However,atransi-tionmodelforgeneralapplicationisnotyetavailable.Furthermore,lowaspectratiowingisusuallyaccompaniedbycomplicatedvorticalowstructures.Thevorticalowstructuresintroduceunsteadyphenomena.ValuableinformationregardinglowReynoldsnumberandlowaspectratiowingswasreviewedbyLissaman[ 61 ]andTani[ 95 ],aswellasaddressedinseveralproceedings[ 65 66 ].Gad-el-Hak[ 20 ]discussedapproachestoprevent/alleviateowseparationonMAVswithmicroelectronicsme-chanicalsystems(MEMS).LianandShyy[ 53 ]performedadetailednumericalstudyoftheowseparationontheMAVwinganditsimpactsonaerodynamicperformance. Second,amembranewingexhibitsself-initiatedvibrationsevenundersteadystatefreestreamcondition.Waszaketal.[ 111 ]conductedwindtunnelexperimentsandfoundthatthemembranewingsexhibitedvibrationatO(102)Hzinasteadyfreestream.SimilarobservationwasreportednumericallybyLianandShyy[ 53 ].Suchvibrationsandtheassociatedshapedeformationchangethepressuredistributionon PAGE 13 themembranewing,which,inturn,aectsthemembranedynamics.Thisrecipro-calactionleadstoauidandstructureinteractionproblem.Aeroelasticanalysismethodcouplingacomputationaluiddynamics(CFD)solverwithacomputationalstructuraldynamics(CSD)solverprovidesaneectivetoolforthemembranewingstudy.AlthoughbothCFDandCSDhaveachievedgreatsuccessindividually,acou-pledcomputationisstillachallengingtask[ 80 ].Todate,thestudyofamembranewinganduidowinteractionsislimited.JacksonandChristie[ 38 ]adoptedathree-dimensionalpotentialow-basedsolverandamembranewingmodeltoanalyzetheaeroelasticbehaviorofmarinesails.SmithandShyy[ 92 ]andShyyetal.[ 88 ]pre-sentedacomputationalapproachtomodeltheinteractionbetweenatwo-dimensionalexiblemembranewingandsurroundingviscousows.Acombinednumericalandexperimentalanalysisofaowoveratwo-dimensionalexiblesailwasconductedbyLorilluandHureau[ 62 ].Inthecoupleduidandmembranewinginteractionstudy,weneedauidsolvertoaccuratelycapturethetransientphenomenonoftheuidowsurroundingthemembranewing;weneedastructuralsolvertopredictthein-herentnonlinearbehaviorsofthemembranematerial.Tocouplethesetwosolvers,wealsoneedtohandledierenttimecharacteristicsassociatedwitheachsolver.Be-sidesthese,weneedamovinggridtechniquetofacilitatetheuidandstructureinteractionproblem.Aninterfacetechniqueisalsorequiredtoexchangeinformationbetweenthesetwosolvers. 1.3 Outline Thisworkstudyboththemembranewingandadaptively-shapedwing.First,wewilldevelopacomputationalcapabilitytostudytheuidandstructureinteraction.Second,basedonthiscapability,wewillstudythepiezo-actuatedappingwingthemembranewing.Third,withthegainedknowledge,wewilloptimizethemembranewingperformancewithmodernoptimizationtechnique.Theoutlineofthecurrentworkisasthefollowing: PAGE 14 InChapter 2 wepresenttheuidowsolverandphysicalmodelsoftheuiddynamics.TheNavier-Stokesequationsforincompressibleowoncurvilinearcoordi-natesareemployed,alongwithatwo-equationk-"turbulencemodel.Inthischapter,theuidsolverisconsideredunderthemovingboundaryframework.Abriefin-troductionisgiventothegeometryconservationlaw[ 100 ]thatplaysanimportantroleinmovingboundaryproblems.Eventhoughthepresentworkdoesnotemploylaminar-turbulenttransitionmodel,abriefreviewisoered. InChapter 3 wedevotetomovinggridtechnology,whichisusedbyboththecoupleduid-structurestudyandtheoptimization.Arobust,ecientandeectivemovinggridtechniqueisproposed,whichusestransniteinterpolation-basedpertur-bationmethodinsingleblockgrid[ 56 77 ],andemploysthemaster/slaveconcept[ 29 ]foramultiblockgrid. InChapter 4 weinvestigatetheactiveowcontrolwithadaptively-shapedapsforalow-Reynoldsnumberwing.Resultsoftheinteractionbetweenpiezo-actuatedapsandsurroundingowarepresented. InChapter 5 wefocusonthestructuralsolver,morespecically,anonlineardynamicmembranemodel.Detailedderivationofthemembranemodel,withanemphasisonthenonlinearityaswellasthetimeintegrationwillbepresented.Themembranemodelwillbeassessedusingwelldenedtestproblems. InChapter 6 weproposeacoupleduidandstructureinteractionalgorithm.Thecouplingisrealizedthroughaninterfacetechniqueandtimesynchronization.Thenumericalimplementationofgeometryconservationlawisgivenandtested.ThelowReynoldsandlowaspectratiowingaerodynamics,includingthelaminarboundarylayerseparation,tipvortices,andtheunsteadyphenomenonareinvestigated.Com-putationalresultswillbepresentedtoelucidatethemembranewingaerodynamics. PAGE 15 InChapter 7 weoptimizeofaexiblemembranewingwithagradient-basedmethod.Sinceadirectoptimizationofamembranewingrequiressubstantialcom-putationalresourceandtime,arigidwingisusedasasurrogate.Themembranewingperformanceisthenevaluatedbasedonthesurrogatemodel. Inbrief,thepresentworkismotivatedby,rst,thepracticalneedtounderstandtheaerodynamicsofmembraneandadaptively-shapedwings,second,ourinterestsinmicroairvehicles.Originalworkincludeamovinggridalgorithmformulti-blockgrid,adynamicmembranemodel,numericalstudyofinteractionbetweenuidandaexiblewing,andshapeoptimizationofaexiblewing. PAGE 16 Ourstudyinvestigatestheaerodynamicsofamembranewingaswellasadap-tivelyshapedwingsunderlowReynoldsnumberconditions.Duringight,themem-branewingundergoesshapechangeundertheexternalforces;meanwhile,itsshapevariationaectsthepatternandstructureofitssurroundinguidow.Themem-branewingvibrationwasobservedbothexperimentally[ 111 ]andnumerically[ 53 55 ].Intheadaptivelyshapedwingstudy,thepiezo-actuatedapsvibratedwithdierentfrequenciesandamplitudes.Theapvibrationscauseunsteadyphenomenalikevor-texshedding[ 30 56 ].Itisdesirabletohaveacomputationalcapabilitytoaccuratelycapturethetransientbehaviorofuidowandstructuraldynamics. Eventhoughtsteadyowcomputationshaveachievedsignicantsuccessintheirdevelopmentandapplications,littleattentionhasbeengiventothesystematicanaly-sisofunsteadycomputations[ 79 80 ].Severalfactorsareresponsibleforthissituation.First,steadyCFDisstillthemainstreamaerodynamicanalysis,andthetechniquesusedtoenhancetheperformanceofsteadyCFDcomputationsmaynotbesuitabletounsteadycomputations.Forexample,insteadycomputations,avarietyoftech-niqueshavebeenproposedtoeliminatethelowfrequencytransientstoacceleratetheconvergencetosteadystate;however,lowfrequencyphenomenaareintrinsictomanyuidandstructureproblems,andsuchtechniqueswillbeinappropriatetotheunsteadyanalysis[ 80 ].Second,theapplicabilityoftheturbulencemodeltounsteadyowsisnotwellunderstood[ 79 ]. However,numerousapproacheshavebeenproposedforunsteadyowcomputa-tions.Pulliam[ 74 ]incorporatedsubiterationtoimprovetimeaccuracyofconven-tionalimplicitschemes.Jameson[ 39 ]andRumseyetal.[ 79 ]developedadualtime 8 PAGE 17 steppingtechniqueinthecontextofmultigridmethodologytoenhancetheeciencyandaccuracyoftraditionalfactoredschemes.Issa[ 35 ],andOliveiraandIssa[ 71 ]pro-posedthePressure-ImplicitwithSplittingofOperators(PISO)methodtoimprovetheconvergence.Forrelevantinformationoftheseandsomeofthemostfrequentlyadoptedapproaches,seeShyyandMittal[ 85 ]. TheSIMPLEmethod[ 72 82 ]andthePISOmethod[ 35 71 98 ]aretwopopularmethodstosolvethethree-dimensionalNavier-Stokesequationsforincompressibleows.Bothmethodsbelongtothepressure-basedapproach[ 82 ]whichdevisesanarticiallyderivedpressure(correction)equationbymanipulatingthemasscontinuityandmomentumequations. Three-dimensionalNavier-Stokesequationsforincompressibleowincurvilinearcoordinatesareadoptedinourstudy.Twomethodsareemployedtosolvethegov-erningequations,oneistheSIMPLECmethod[ 102 ]thatisbelongtotheSIMPLEfamily,andtheotheristhePISOmethod[ 35 ].Inthischapterwebrieyintroducetheseuidsolversunderthemovingboundaryframework.Fordetaileddescriptionsofthesemethods,werefertoRefs.[ 35 72 82 98 102 ].Tosolvethemovingbound-aryproblems,wealsoemphasizetheimportanceoftheboundaryconditionsandthegeometricconservationlaw[ 100 ].Anintroductionofturbulencemodels,withanemphasisonthek-"turbulencemodelispresented.Thoughanaccuratetransitionmodelforgeneralapplicationisbeyondthescopeofourstudy,weoerabriefreviewofthetransitionmodelsattheendofthischapterforcompletenesspurpose. 2.1 GoverningEquationsandFluidSolver 2.1.1 GoverningEquationsforFluidFlow Thethree-dimensionalNavier-Stokesequationsforincompressibleuidowsincurvilinearcoordinates[ 82 ]instrongconservativeformcanbewrittenasthefollow-ing, @+@V @+@W @=0;(2.1) PAGE 18 @[ J(q11u+q12u+q13u)]+@ @[ J(q21u+q22u+q23u)]+@ @[ J(q31u+q32u+q33u)][@ @(f11p)+@ @(f21p)+@ @(f31p)]+G1(;;):J; @[ J(q11v+q12v+q13v)]+@ @[ J(q21v+q22v+q23v)]+@ @[ J(q31v+q32v+q33v)][@ @(f12p)+@ @(f22p)+@ @(f32p)]+G2(;;):J; @[ J(q11w+q12w+q13w)]+@ @[ J(q21w+q22w+q23w)]+@ @[ J(q31w+q32w+q33w)][@ @(f13p)+@ @(f23p)+@ @(f33p)]+G3(;;):J; whereu,v,andwarethevelocitycomponentsinthex-,y-,andz-directions,re-spectively,istheviscosity,G1,G2,andG3arethebodyforcecomponentsperunitvolumeinthex-,y-,andz-directions.Thecoecientsappearingintheequations PAGE 19 areexpressedasfollows: Theexpressionsforfij,i=1,3,j=1,3,canbefoundin[ 82 ].TheJacobiantransfor-mationalmatrixisdenedasfollows: ThedeterminantoftheJacobianmatrixisexpressedas Inthegoverningequations,Uistheuxthroughthecontrolsurfacewhosenormaldirectionis;VandWaretheuxesthroughtheandcontrolsurfaces.Whenthegoverningequationsareconsideredunderamovinggridframework,thegridvelocitiesshouldbeconsideredintheuxcomputations.Incurvilinearcoordinates,theyareexpressedasthefollowing: where,forexample,_xisthegridvelocitycomponentinthex-directionthatcanbeapproximatedbytherstorderEulermethod, _x=xx0 wheretisthetimestepandthesuperscript0referstotheprevioustimelevel.Similarly,thevaluesoftheothertwocomponents_yand_zcanbeestimated. PAGE 20 2.1.2 FluidFlowSolvers TheSIMPLECmethod[ 102 ]isavariantoftheSIMPLEmethodoriginallypro-posedbyPatankar[ 72 ]inCartesiancoordinates,anditsextensiontocurvilinearcoordinateswasdiscussedbyShyy[ 82 ]toaccommodatethecomplexgeometrycom-putation.TheSIMPLECmethodusescoordinatedunder-relaxationsforthemomen-tumandpressurecorrectionstoimprovetheconvergencethatisintrinsicallyslowintheoriginalSIMPLEmethod.TheSIMPLECmethodwasalsousedtosolvemov-ingboundaryproblems[ 53 88 ].ContrarytotheSIMPLEfamilymethod,thePISOmethodisanon-iterativeapproachtohandlethepressure-velocitycouplingwhichsplitstheprocessofsolutionintoaseriesofpredictorandcorrectorsteps.ThePISOmethodisgenerallymoreecientfortransientowcomputations[ 35 98 ]. TosolvetheNavier-Stokesequations,werequireproperboundaryconditionsforthecomputations.Attheinterfacebetweentheuidandthesolidstructure,thefollowingboundaryconditionsareapplied, Eq.( 2.10 )expressesthecompatibilitybetweenthewetuidsurfacexfandthestruc-turesurfacexsattheinterface;Eq.( 2.11 )requiresthattheuidvelocityattheinterfacematchesthesurfacevelocityofthestructure.Theseboundaryconditionsareprovidedbycomputingthedisplacementofthesolidstructure. 2.1.3 GeometricConservationLaw Whenbody-ttedcurvilinearcoordinatesareusedinthecomputation,atrans-formationmapsaphysicalowregion(x,y,z)ontoauniformcomputationalspace(,,)whereconservationlawsarecarriedout.Tofacilitatethecoordinatetrans-formation,theJacobianmatrixJisintroduced.TheJacobianmatrixisdenedas PAGE 21 inEq.( 2.6 ).ThedeterminantofJrepresentsavolumeelementinthetransformedcoordinates.Inamovinggridproblem,thecomputationalgridneedstobeupdatedwithtime;meanwhile,theJacobianJneedstobeupdatedsimultaneously.SpecicproceduresarerequiredtocomputetheeectivevalueofJ;otherwise,errorsariseduetotheinconsistencyinevaluatinggeometricquantitiesandtheaddingofnon-physicalquantities.AspointedbyThomasandLombard[ 100 ],intheupdatingprocess,thefollowinggeometricconservationlaw(GCL)neededtobesatised, dtZVJddd=ZV(rWs)ddd;(2.12) whereVisthevolumeboundedbyaclosesurfaceS,WsdenotesthelocalvelocityoftheboundarysurfaceS.ThomasandLombard[ 100 ]proposedtoevaluateJwhilemaintainingthegeometricconservationlawbysetting=1,V=0fromthecontinuityequation,i.e., Implicationsofthegeometricconservationlawwithuid-structureinteractionprob-lemswerediscussedbyShyyetal.[ 88 ]. 2.2 TurbulenceModels 2.2.1 ReynoldsEquationsandk-"TurbulenceModel Forsimplicity,theReynoldsaveragedmomentumandcontinuityequationsforincompressibleowsarepresentedinCartesiancoordinates whereTijisthestresstensor,repeatedindicesinanytermindicateasummationoverallthreevaluesoftheindex,andUiistheReynoldsaveragedmeanvelocity.The PAGE 22 stresstensorinaturbulentowmaybewrittenas 2(@Ui k= Inmanyows,thesenormalstressescontributelittletothetransportofmeanmo-mentum.Theo-diagonalcomponentsofijareshearstresseswhichplayadominantroleinthemeanmomentumtransferbyturbulentmotion. Eq.( 2.14 )containsnineunknowncomponentsofij(ofwhichsixareindependentofeachother)besidesPandthethreecomponentsofUi.Therefore,weneedmoreequationsforijtoclosethesystemofequations.Aturbulenceclosuremodelservesthatpurpose. Dierentversionsofturbulenceclosuremodelshavebeenproposed.Thesuccessofsuchmodelsislargelydependontheapplicationsofthewallfunctionsthatrelatesurfaceboundaryconditionstouidpointsawayfromtheboundaries.Therebyweavoidtheproblemofdirectmodellingtheviscosityinuence.HereweusethewallfunctiontreatmentbyLaunderandSpalding[ 52 ]asoneexampletoillustratetheclosureproblem.Thegoverningequationsforstandardk-"modelscanexpressedin PAGE 23 theformshownbelow, @xi(t @xi)+2tSijSij";(2.20) @xi(t @xi)+2C"1t" kSijSijC"2"2 Theeddyviscositytisdenedas whereCisadimensionlessconstant.Sijisthemeanrateofstrain. Theseequationshaveveadjustableparameters,C,k,",C"1,andC"2.Here,forbrevity,weonlygivethecoecientsforstandardk-"modelsuggestedbyLaunderandSpalding[ 52 ].Theseparametersarevalidforawideofrangeofturbulentows. Toclosethesystemofequationsforturbulencecomputation,weneedanadditionalBoussinesqrelationshipthatrelatestheReynoldsstresstothemeanrateofstrainandturbulencekineticenergy, uiuj=t(@Ui 3kij=2tSij2 3kij:(2.24) 2.2.2 BoundaryConditionsforTurbulenceModels Theellipticcharacteristicofthekand"modelequationsgivesriseofbound-aryconditions.Usuallyweassumezerogradientattheoutletboundary.Atinletboundary,ifthereisnomeasurementofkand"atdisposal,onewayistomakeanapproximationbasedonthefreestreamvelocityU1,turbulenceintensityTi,andthecharacteristiclengthscaleL.Inourstudyweadoptthefollowingformulastodene PAGE 24 thefreestreamboundaryconditionsforkand",k=3 2(U1Ti)2; TheturbulenceintensityTiisusuallysetbetweentherangeof0.02and0.05.AnotherwayistoadjusttheturbulenceReynoldsnumberRet=U1L=t.ForpracticeweusuallysettheturbulenceReynoldsnumberintherangeof100and500. Thesolidwallboundaryconditionsaremorecomplicated.Nearthewallbound-aries,thelocalReynoldsnumbersarelowandtheviscouseectsareimportant.However,thek-"turbulencemodelthatwementionedaboveisbasedonthehigh-Reynolds-numberassumptioninwhichtheturbulentstressforcesaredominated.Moreover,closetothewall,extremelynegridisusuallyrequiredtocapturethesteepvariationsinowprolestosolvetheNavier-Stokesequationsdirectly.TheCPUtimebasedonsuchanegridisalmostbeyondtheabilityofcurrentfacilities.Inordertosolvethesediculties,twoapproachesarecommonlyemployed,oneisthewallfunctiontreatment[ 52 ],theotheristhelow-Reynolds-numberversionsofthek-"[ 46 73 ]. Thewallfunctiontreatmentispopularinpracticebecauseofitssimplicityandcapabilityforcomplexgeometries.UnderthehighReynoldsnumbercondition,owcanbelargelydividedintotwolayers:theouterlayerandthesurfacelayer.IntheouterlayertheReynoldsstressisdominated;whileinthesurfacelayerboththeReynoldsstressandtheshearstressareimportant.Intheouterlayer,theuideldisinuencedbytheretardingeectofwallthroughthevalueofthewallshearstress,butnotbytheviscosityitself. Inthewallfunctiontreatment,velocityproleshavedierentdescriptionsintheoutandsurfacelayers.Thevelocityintheoutlayerisgovernedbythevelocity-defect PAGE 25 lawexpressedas u=g(y );(2.27) whereistheboundarylayerthickness,Umax-Uisthevelocitydecit,gisanunknownfunctionwiththeorderof1.Thevelocityinthesurfacelayerisgivenbythelawofthewall u=f(y+):(2.28) Herey+isanon-dimensionalnormaldistancefromthewall,anditisrepresentedby :(2.29) Intheaboveexpressionyisthedistancefromthesolidwall,anduisthefrictionvelocityatthewalldenedby Theshearstressatthewall,w,isgivenby @y:(2.31) Thelow-Reynolds-numberk-"modelrequiresarelativelynegridresolutionnearthewallregions.FordetaileddescriptionswerefertoRefs.[ 46 73 87 ]. Tomakeappropriateuseofthetwoversionsofturbulencemodels,westrivetoadjustthevalueofy+oftherstpointnexttothewalltosatisfycertainrequirements.Forthelow-Reynolds-numberapproach,y+shouldbelessthan2foranaccurateresult,whiley+intherangebetween30and500[ 103 ]canproduceaccurateresultsforthewallfunctiontreatment.Anequationfortheatplateboundarylayertheoryisusedinthegridgenerationprocesstogetanestimationofpropercellsizefory+, PAGE 26 Italsoshouldberemember,whilethepositionoftherstcelliscritical,itmusthaveenoughcellsinsidetheboundarylayertoresolvethesteepvariationsofowvariables. ThesuccessoftherecentturbulenceclosureisrestrictedtothesucientlylargeenoughReynoldsnumberturbulenceows.Furtherrenementisneededbeforetheyareusedwithcondencetocalculatenear-wallandlowReynoldsnumberows.Amoredetailedreviewaboutturbulencemodelsfornear-wallandlowReynoldsnumberowsisgivenbyPateletal.[ 73 ]. 2.3 Laminar-turbulentTransition Areliablecomputationoftheboundarylayertransitionremainsoneofthemostchallengingdemandsinaerodynamics.Sincethetransitionalseparationbubblesandtheirlossesplayanimportantroleforthepressureanddragdistributions,anac-curaterepresentationofbothlaminarandturbulentseparatedowsiscriticalwhendragpredictionneedstobeaccuratelypredicted.However,thefundamentaluiddy-namicsproblemoftransitiontoturbulencehasresistedanysimplephenomenologicaldescription.Thereisstillnocomprehensivetheoryoftransition. Therangeofexistingtransitionpredictionmethodsextendsfromsimpleempiri-calrelationshipsthroughthosebasedonparallelandlinearstabilitytheories,liketheeNmethods[ 15 75 ],orlinearornonlinearparabolizedstabilityequations(PSE)[ 31 ],orthelow-dimensionalmodels[ 76 ],todirectnumericalsimulation,likelarge-eddysimulationtechniques. PAGE 27 Forproblemsinvolvinginteractionsbetweenuidowsandmovingstructures,thegeometryofthesolidobjectisnotknownapriori.Inthecourseofcomputationitisnecessarytotrackthegeometry,theeldequations,andtheinterfacialconditiontoensurethatallrequirementsaremet.Tofacilitatethesolutionofsuchmovingbound-aryproblems,themovinggridtechnique(ordynamicgridtechnique)isemployedtoadjustthecomputationalgriddynamicallyalongwiththegeometricupdates.Itisdesirablebutdiculttodevelopanautomatedregriddingproceduretoensurethatthenewgridnotonlymatchesthegeometricchangesbutalsomaintainssatisfactorycharacteristicssuchasnon-overlapping,smooth,andnotexcessivelyskewed.Themovinggridtechniquecanalsobeusedbyshapeoptimization[ 54 ]. Formovinggridproblems,severalalternativeapproacheshavebeenproposed.Schusteretal.[ 81 ]usedanalgebraicshearingprocessintheirstaticaeroelasticanal-ysisofaghteraircraft.Thebasicideaistoassumethatthedisplacementofthesubjectisredistributedalongthegridline,whichconnectstheaeroelasticsurfacetotheouterboundary.Thissimplemethodworksquitewellforsingleblockgridwithmodestdeection;however,whenappliedtomultiblockgridarrangements,thisapproachusuallyneedsextensiveinterventions.Thespringanalogymethodwasrstproposedforunstructuredgrids[ 3 ],andlaterwasextendedtostructuredgridsbyRobinsonetal.[ 78 ].Thismethodregardsallgridsasconnectedbysprings,andeachspringhasthestinessinverselyproportionaltothespacingbetweentargetvertexanditsneighbors.Comparedwithothercurrentlyusedgridregenerationmethods,thespringanalogymethodneedsmorememoryandCPUtime.Directtransnite 19 PAGE 28 interpolationbyEriksson[ 18 ]isapopularmethodbecauseofitseciencyforstruc-turedgrid.HartwichandAgrawal[ 29 ]alsoproposedtheMaster/Slaveconcepttoexpeditethegridregenerationprocessandminimizetheuserintervention. 3.1 PerturbationMethod Forsmallamplitudedisplacementononefaceofthecomputationalblock,asimplebutecientone-dimensionalperturbationmethod[ 56 77 ]iscanecientlyregeneratethegrid.Theone-dimensionalperturbationmethodobeysthefollowingformula, wherexirepresentsthelocationofanarbitraryinteriorgridpoint,xsrepresentsthelocationofagridpointonthemovingboundary,andSrepresentsthenormalizedarclengthalongtheradialmeshlinemeasuredfromtheouterdomain,whichisgivenby Tousethisone-dimensionalperturbationmethod,oneneedstoknowthedisplace-mentsofthetwoendingpointsofameshline.ThepositionsoftheinterimpointsaresolelydeterminedbythedisplacementsoftheseendingpointsandtheiroriginalpositionsasshowninEq.( 3.1 ).Thesedisplacementsareconsideredasperturbationsources.Intuitively,theperturbationwillspreadthroughoutthewholedomainwhileexertingmoreeectsonnearerpoints.BasedonEq.( 3.1 ),amorecomplicated3-stageperturbationmethod,analogoustothetransniteinterpolation(TFI)method[ 18 ],wasproposedforthree-dimensionalproblemstreatedbythesingle-blockapproach.UnliketheoriginalTFImethod,thisperturbationmethodconsiderstheoriginalgriddistribution,andpreservestheoriginalgridquality. Takeatwo-dimensionalcaseasanexample,inFig. 3{1 ,arectangleundergoesshapechange.Toregeneratethegridforthechangedgeometry,theperturbationmethodproceedsasa2-stageprocess.Therststageistomatchcornerpoints,i.e., PAGE 29 tomovethecornerpointAtoitsnewpositionA0,andBtoB0,etc.(SeeFig. 3{1 b).Oncethecornerpointsmatchtheirnewpositions,interimedgesaregeneratedbasedonEq.( 3.1 ).ToadjustthepositionofaninternalpointO,eectsfromthefourcornerpointsneedtobeaccountedforinacoordinatedmanner.Inthisapproach,P0,Q0,S0,andT0,whichconnectO0toitsboundary,willberstadjusted.Thedisplacementsofthefouredgepointsin,e.g.,thex-direction(P,Q,S,T)canbecalculatedwithEq.( 3.1 ),andthemovementofO0isevaluatedasacombinationofthesedisplacements.Specically,onecanadoptthefollowingexpression: O0=abs(I)I+abs(J)J where,I=(1Si)P+SiQ,J=(1Sj)S+SjT,SiandSjarethenormalizedarclengthsalongmeshlinesPQandSTrespectively.Displacementsinotherdirectionscanbecomputedinthesimilarway. Aftertherststage,thefourcornerpointshavematchedtheirnewlocations;however,theinterimedgesdonotnecessarilymatchtheirnalpositions.Thesecondstageistomatchtheinterimedgestotheirnaledgelocations.Thedierencebe-tweentheinterimandnaledgelocationsaretakentobethesourceofperturbation.Inthisstage,sincetheperturbationsinbothdirectionsareindependent,thecontri-butionsfromallindividualdirectionsareadded[SeeFig. 3{1 c].Forthree-dimensionalproblems,a3-stageprocessperturbationmethodwasproposedinRefs.[ 56 77 ]. 3.2 Master/SlaveConcepts Theaforementionedperturbationmethodcanbeimplementedecientlyforsin-gleblockandsmalldisplacementproblems.However,whenamulti-blockgridisinvolved,thismethodneedstobeenhancedbyothertechniques.TheMaster/Slaveconceptactsasausefulapproachtomaintaingridqualityandtopreventpotentialgridcross-over. PAGE 30 Figure3{1: Twostageperturbationmethod:(a)ABCDchangestoA0B0C0D0;(b)rststage:matchingcornerpoints;(c)secondstage:matchingedgepoints. PAGE 31 Althoughunstructuredgridhasbeenusedtosolveuiddynamicsproblemsbydierentresearchers[ 3 26 ],here,ourattentionislimitedtostructuredgridonly.Foramulti-blockstructuredgrid,CFDsoftwareoftenrequirespoint-matchedgridblockinterfaces.Abasicrequirementforgridregenerationinmovingboundaryproblemsistomaintainpoint-matchedinterfaces.Whilethisupdatingprocessiseasilysatisedbythoseblocksurfacesthatcoincidewiththebodysurface,theremainingblocksurfacesneedtobehandledwithduecare.Tousetheperturbationmethod,oneneedstospecifythenewpositionsoftheo-bodysurfacepoints,or,atleast,theverticesofsuchblocksurfaces. Whenanobjectchangesitsshape/position,masterpoints,whicharedenedasthepointslocatedatthebodysurface,rstmove,andthenaectdistributionsoftheslavepoints,theo-bodypoints.Whilethetaskofredistributinginteriorgridscanbechallenging,inpractice,itismorediculttomovetheverticesofeachblockifapoint-to-pointmatchedinterfacewithoutoverlapisrequired.Iftwoabuttingblockshavecongruentinterfaces,suchastheinterfaceofblock2andblock3inFig. 3{2 a,oncethecornerverticesaredetermined,theedgepointsandtheinteriorpointscanbeuniquelysettledbasedonEq.( 3.1 ).However,whentwoneighboringgridblocksdonotshareidenticalinterfaces,suchasblock1andblock2showninFig. 3{2 a,thisprocedurecancausediscontinuityattheinterfacebecausetheregridingproceduresattheinterfacearebasedondierentstencilsfordierentblocks.Slavingverticestotheirmasterpointscanavoidcreatingundesirablegriddiscontinuities. TheMaster/Slavecouplingisbasedonthedistancebetweenano-bodypoint(Slavepoint)andasurfacepoint(Masterpoint).Thedistanceisgivenby: wheresubscriptvrepresentsano-bodyvertex,and,b,abodysurfacepoint.Foreacho-bodyvertex,thenearestbodysurfacepointisdenedasitsmasterpoint. PAGE 32 Figure3{2: Testsonthemovinggridscheme:(a)initialgeometry;(b)sideviewofthegrid;(c)deformedshape;(d)newgriddistribution;(e)griddistributionwithhighstinessparameter. PAGE 33 Tosimplifytheconnectionbetweenthemasterandslavepoints,onecanconvertthethree-dimensionalparameterspace(i,j,k)intoaone-dimensionaldatastructure.Eachgridpointhasanidenticationnumberldenedby wherebo(bn)istheidenticationnumberoftherstpointinblockbnintheone-dimensionalarray,imaxandjmaxarethemaximaldimensionsintheiandjdirections,respectively.Inthiswayaslavepointisassociatedwithitsmasterpointbystoringitsmaster'sidenticationnumberintoitsaddress. Oncetherelationisestablished,thenextstepistodeterminehowtheslavepointmovesaccordingtotheinuencefromitsmasterpoint.Intuitively,anearervertexwillhavemoreeect(displacement)thanthosefaraway.Also,toavoidalowersurfacetocrossoveritsoppositefacewhenthemovementisrelativelylarge,oneneedstoscalethemovement.Asimplebuteectiveformulaisexpressedasfollows: ~xs=xs+(~xmxm);(3.6) wheresubscriptsmandsrepresentthemasterandslavepoints,respectively,thetilde()indicatesthenewposition,andisadecayfunction.Inourcomputations,aGaussiandistributionfunctionisused whereisaparameterwhichaectsthestiness. PAGE 34 3.3 NumericalTests Themovinggridtechniqueconsistsoftwoelements,oneistheperturbationmethod,andtheotheristhemaster/slaveconcept.Theproposedmethodistestedonathree-dimensional,multi-block,andstructuredcomputationalgrid.TheinitialcongurationisshowninFig. 3{2 a,block1sharestheinterfacewithblocks2and3.TheinitialgridisshowninFig. 3{2 b,ithasclustereddistributionsatthetopandbottomsurfaces,anditalsohaspoint-to-pointmatchedinterfaces.Thebottomtwofacesarethemovingboundariesthatexperiencelargedeformation.Thecoordinatesofeachvertexarelabelled.ComparingFig. 3{2 aandFig. 3{2 cwecanseethedisplacementisaslargeashalfofthetotalheight.Thenewcomputationalgridiscomputedwiththemethoddiscussedinthischapter.ThenewgridisshowninFig. 3{2 d.Thatguredoesnotshowanycross-overinthenewgridafterthelargedeformation.Thenewgridnotonlykeepsthepoint-to-pointmatchedinterfaces,butmaintainstheclustereddistributionsatthetopandbottomsurfaces.TheCPUtimeforregeneratingthecomputationalgridislessthan10%oftheCPUtimeusedbytheuidsolver. Theparametercontrolsthestiness,wend=1/64workswellforourprob-lems.Alargercausestheblocktobehavemorelikearigidbody(SeeFig. 3{2 e).MoredetailsaboutthemasterandslaveconceptcanbefoundinRefs.[ 29 56 ]. PAGE 35 4.1 Introduction Themainobjectivesforuiddynamistsandaeronauticalengineersthroughoutthehistoryofmodernighthavealwaysbeentoincreaseliftandreducedrag.Microairvehiclesareairbornevehicleswithcharacteristiclengthsunder6inchestravelingatspeedsaround30ft/s.Theyaredesignedtoperformanarrayofdierenttaskssuchasreconnaissance,targetdesignationorbordersurveillance.Toenhancetheperformanceofthesesmallvehicles,improvementsshouldbemadeintherespectsofpowerplants,thrustgenerators,navigationsystems,cameras,datatransmitting,andmoreimportantly,theiraerodynamicperformance.Asfarastheaerodynamicsisconcered,MAVshavetooperateinoneoftheworstightregimesimaginable,oneofthemisthelow-Reynoldsnumberswhichisusuallynorotiouslyaccompaniedwithlaminarboundaryseparation,transition,andlowlift-to-dragratio.FurthermorethesmalldimensionsoftheMAVsgreatlyenhancesthetradeodicultiesmentionedabove.Thesmallsizeinhibitsthedesigner'sfreedomtoachievethedesiredaerody-namicperformance.Thelimitationsinthepowersourcesavailable,atpresenttime,makethisequationevenmorecomplicated. Tosolvetheaerodynamicchallengesoftheseproblems,twomajoreldshavebeenevaluated.Therstapproachistryingtocopynature'sownway;thatofappingwings[ 43 45 59 83 ].Thesecondapproachisusingadynamicalwing.Therearesomedierentimplementationsinthiseld.Onewidelystudiedmethodistheuseofanwingwithaexiblemembraneuppersurface[ 53 55 83 92 ].Anothermethodusesactivedevicestocontroltheowontheuppersurfaceofthewing.Thelattermethodhasgainedmoreattentionwiththesuccessmicroelectromechanicalsystems(MEMS). 27 PAGE 36 Inactiveowcontrolarea,apiezo-actuatedapmountedonalowReynoldsnumberwinghavebeenstudiedbothexperimentally[ 19 ]andcomputationally[ 30 56 ].Gad-el-Hak[ 20 ]hasdiscussedthepossibilityofemployingMEMStodelayorpreventseparationonMAVs. TounderstandthelowReynoldsnumberwingcharacteristicsandimprovethewingperformance,wewillperformeacomputationalinvestigationforowssurround-ingadynamicallyshapedwing,atachordReynoldsnumberof78,800,alongwithaparallelexperimentaleort.Weadoptapiezo-actuatedapontheuppersurfaceofaxedwingforactivecontrol.Theactuationfrequencytheyfocusonis500Hz.Theircomputationalframeworkconsistsofamulti-block,movinggridtechniquetohandlethegeometricvariationsintime,usingtwo-equationturbulenceclosures,andapressure-basedowsolver.Thedynamicgridredistributionemploysthetransniteinterpolationschemewithaspringnetworkapproach.Comparingtheexperimentalandcomputationalresultsforpressureandvelocityelds,eectsofthedetailedge-ometriccongurationoftheap,theappingamplitude,turbulencemodeling,andgriddistributionsontheowstructurewillbeassessedinthischapter. 4.2 ExperimentalandComputationalSet-up ThewinghasmovingapsontheuppersurfaceasshowninFig. 4{1 .Theapcanbeactuatedatprescribedfrequencies[ 19 ].Thewingstudiedintheexperimenthasaspanof1ftandachordof5:400.Intheexperimentalset-up,thereisagapof0:01600widebetweentwoaps,eachis200wideand0:56400long.TheReynoldsnumber,basedonthefreestreamvelocityU=27.5ft/s,is78,800.Theowisinitiatedaslaminar;afterthetransitionpoint,theturbulencemodelisinvoked.Thetransitionpointisidentiedtobeatthetipoftheap.Inthecomputation,theturbulentReynoldsnumberbasedonthechordlength,Ret=UL=t,attheinletoftheturbulentowregionissetto500.Theamplitudeandshapeoftheapareprescribedbyanalyzing PAGE 37 theexperimentaldataobtainedwithPIV[ 19 ]Forthe500Hzcase,twoappingamplitudesareconsidered,namely0:00100and0:007300. Figure4{1: Movingapwingexperimentalsetupinthewind-tunnel Inthischapter,anumberoftwo-dimensionalandthree-dimensionalcases,in-cludingdierentgriddensities,withandwithoutthegap,willbestudied.For3-Dcases,aneorthasbeenmadetoconsidertheeectofthegapbetweenapsasshowninFig. 4{2 .TheoriginalLaunder-Spaldingk-"model[ 52 ]withthewallfunction,aswellasalow-Reynoldsnumberversionk-"model[ 46 ]arecompared. Figure4{2: Flapwingwithrenedthreedimensionalbendingmodel. 4.3 NumericalStudy 4.3.1 GridSensitivityAnalysis Anumberof2-Dgridswithdierentdistributionsandsizeshavebeentested.Thenareference2-Dgridisselectedbyensuringthatsatisfactorysolutionqualityisobtainedonamosteconomicgridsize.Then,the3-Dgridisestablishedbyextendingthereference2-Dgridalongthespanwisedirection.Theresultinggrid,forcaseswithgap,has33nodesinthespanwisedirection,andconsistsof20blockswithapproximately533,000gridpoints.Forcaseswithoutgap,thegridsystemhas PAGE 38 20nodesinthespanwisedirection,andconsistsof8blockswithabout253,000gridpoints.Withthemulti-blockstrategyitisconvenienttoemploythetransitionmodel,i.e.,whenatransitionlocationisidentied,newcomputationalblocksaregeneratedsurroundingthetransitionpoint,sothatthelaminarowequationsaresolvedintheblocksbeforethetransitionlocation,andtheReynolds-averagedturbulentowequationsaresolvedintheblocksafterthetransitionpoint.Inourcomputations,thetransitionpointisidentiedfromtheexperiments[ 19 ],anditisatthetipoftheap,therefore,ourcomputationalblocksaresetupbasedonthistransitionpoint. 4.3.2 AssessmentoftheTurbulenceModel IntheearlyworkbyHeetal.[ 30 ],thewallfunctionapproachisusedexclusively.ForsuchalowReynoldsnumber,thelow-Reynoldsnumberversionofthek-"modelseemsmoresuitable[ 56 ].Fig. 4{3 comparesthetime-averagedpressuredistributionsat500Hzbetweentwodierentturbulencemodelsforthecasewithoutgapbetweenaps.Bothexperimentalandcomputationaldataarecollectedatthemiddleoftheapandalongthechordwisedirection.Asfarasthecomputationaltimeisconcerned,itisobservedthatthelow-ReynoldsnumbermodeluseslessCPUtimethanthatofthewallfunctiontreatment. 4.3.3 EectsofGapsbetweenFlaps Theinclusionofagapobviouslyincreasesthecomplexityoftheoweldinthespanwisedirection.ItcanbeseenfromFig. 4{4 thatstreamlinesenterthegapandmovetowardthemiddleoftheap.Aweakrecirculataryowwithavortexstretchingfromthewingsurfacetotheapispresentonbothsidesofthegap,rotatingintheclockwisedirectionontherightandcounter-clockwisedirectionontheleftofthegapwhenviewedfromabove. Thegapalsoaectstheoverallpressuredistributiononthewing.Fig. 4{5 comparesthepressuredistributionsobtainedwithandwithoutgapbetweenaps.Thelow-Reynoldsnumberturbulencemodelisusedinbothcases.Timeaverageis PAGE 39 Figure4{3: Comparisonbetweentime-averagedpressuredistributionsforcaseswithoutgap,at500Hz,obtainedwithdierentturbulencemodels.(Timeaveragedfortheentirecycle,amplitude:0:00100). Figure4{4: Streamlinefor3-Dcasewithgapat500Hz,amplitudeis0:00100 PAGE 40 madefortheentirecycleforanamplitudeof0:00100andanactuationfrequencyof500Hz.Sincethepressuretapsarelocatedatthemiddleofaap,theeectofthegapmaynotbeclearlydemonstratedbytheguresshown.Fig. 4{6 comparesthevelocityprolesimmediatelydownstreamfromtheaptrailingedge.Thelow-Reynoldsnumberturbulencemodelisusedagainandthecomparisonismadeat0:016200downstreamfromthetrailingedgeoftheap.Casesinvolvingbothstationaryapsandactuatedaps,at500Hz,areexamined.Consistentwiththatalreadyobserved,thevelocityprolesindicatethattheowenterintothegapfromtheuppersurfaceoftheap,maintainingthesamestreamwisedirectionwhilemovinglaterallytowardtheregionbetweentheapandthewingsurface.Arecirculatingowpatternisformedbetweentheapandthewingsurface. Figure4{5: Comparisonoftime-averagedpressuredistributionsforcaseswithandwithoutgapbetweenaps. 4.3.4 EectsofFlappingAmplitude Dependingontheappingamplitude,thereattachmentpointcanmovebackandforth,andtherecirculatingowcaneitherremainattachedtothetipoftheap,orbeshedaway.Inourcomputationthecurvatureoftheapissetbycollecting PAGE 41 Comparisonofvelocityprolesattwodierentspanwisepositionswithgapimplemented:(a)steadycomputation;(b)f=500Hz. theexperimentaldataandthenapproximatedbyasimpliedmodel.Forthe500Hzcaseinourcomputationtwoamplitudesareconsidered,oneis0:00100andtheotheris0:007300.Forthesmalleramplitudecase,themotionofthetrailingedgeismodeledasthefollowing, wherezisthespanwisedistancefromthestudypointtothegap,fisfrequency,andtistime.Thelargeramplitudecasesimplyscalestheaboveequationbasedonthesameexpression. Fig. 4{7 comparesthetime-averagedpressurecoecientsbetweenthetwoap-pingamplitudes.Forthehigherappingamplitude,thetime-averagedreattachmentpointmovesclosertothetrailingedgeoftheap,becausealargermovementoftheapintroducesmorehighmomentumuidsintotheboundarylayer.Intheexperi-mentweobservevortexsheddingasshowninFig. 4{8 a.Inourcomputationswiththesmallerappingamplitude,thevortexcoreisanchoredattheaptip,exhibitingtwoco-rotatingcoresastheapsmoveupwards,asshowninFig. 4{8 b.Ontheotherhand,withthehigheramplitude,vortexsheddingappearsforcaseswithorwithout PAGE 42 Figure4{7: Pressurecoecientatdierentactuationamplitudewithoutgapbetweenaps.Actuationfrequency:500Hz. thegap.Fig. 4{9 showsselectedinstantaneousplotsofthevortexstructureduringasinglecycleoftheapmovement.Asexpected,therecirculationisstrongerwhentheapmovesupwardandbecomesweakerwhentheapmovesdownward. Instantaneoustrailingedgestreamlinesforthefrequencyof500Hzasapismovingupwardsforamplitudeof0:00100.Left:experimentalresult;right:computationalresult. 4.4 Conclusions Baseonourresultswendthatthelow-Reynoldsnumberformulationofthek-"turbulencemodelproduceslongercirculationzone.Thegapbetweentheapsaectsthelocationofthere-attachmentpoint,andcreatesathree-dimensionalowpattern.Aidedbythelimitedexperimentalinformation,thecomputationalmodelcanoer PAGE 43 Trailingedgestreamlinesforthefrequencyof500Hzinthethirdcyclewithanamplitudeof0:007300 Astothephysicalmechanisms,whileitisexpectedthatthevortexsheddingappearswhentheamplitudeishigh,itisnotclearwhyandwhenthisphenomenonceasestoexistwhentheamplitudeisreduced.Similarly,whileitisobviousthatthegapbetweenapsaectstheowreattachment,itishardtopredictthetrendquantitatively. PAGE 44 5.1 Introduction Thedeformationsofmembranestructuresareoftenlargeandthusinherentlynonlinear.Qualitativedescriptionsofthethinmembranebehaviorsareoftenmorecomplexthanthatofthree-dimensionalbodies.Untilrecently,membraneanalyseswereprimarilyreliedontrialanderror.GreenandAdkins[ 25 ]arethersttoderiveageneralnon-lineartheoryformembranes.Themembranetheoryhastwofundamentalassumptionsbecauseoftheirthinness[ 41 ].First,themembranetensionstinessismuchlargerthanitsbendingstiness,therefore,thestresscoupleeectscanoftenbeneglected.Second,theratioofthemembranethicknesstoitsradiusofthecurvatureissmall,whicheectivelydecouplesthestrain-displacementrelationfromthecurvaturetensor.Practicallyspeaking,thereexisttwoapproachestomodeltheresponseofamembranestructure.Therstapproachistoidealizeaplateorshellstructureasamembraneawayfromitsboundaries,inwhichthestresscouplesareneglectedintheinteriorregion.Byformallyequatingtheplatestinesstozero,thismembranemodelcanbedescribedbytheFoppl-vonKarmantheory[ 40 41 ].Inthesecondapproach,thestructurecannotsustainstresscouplesanywhere.ThesecondapproachisadoptedinthecurrentworkbecauseitisappropriateforMAVcomputations. Two-dimensionalmembranewingstudywithxedleadingandtrailingedgeswereoriginallyconductedbyNelsen[ 68 ],Thwaites[ 101 ],Voelz[ 108 ].Theseworksconsideredinextensiblemembranessurroundingbysteady,two-dimensional,andir-rotationalows.Asaconsequenceoftheinextensibleassumption,whencambersandincidenceanglesaresmall,problemscanbelinearizedandexpressedcompactlyina 36 PAGE 45 integralformwhichisreferredasthe"sailequation"byNewman[ 69 ].Thefollowingequationcompletelydenesthelinearizedtheoryofinextensiblemembranewingsinsteadyandinviscidowelds. 1CT wherey(x)denesthemembraneproleasafunctionofthexcoordinate,istheincidenceangle,andCTisthetensioncoecient.BasedontheworkofVoelz[ 108 ],thesailequationhasbeensuccessfullysolvedanalyticallyandnumerically.Newman[ 69 ]gaveacomprehensivereviewoftheearlierworksrelatedtomembranewingaerodynamics. Thedevelopmentofathree-dimensionalmembranemodelintroducesseveralcomplicatingfactorsnotencounteredintwo-dimensionalanalyses.First,unlikethe\sailequation",thetensioninathree-dimensionalmembraneisnolongerasingleconstantvaluebutisatbestdescribedbyastateofbiaxialtensionalongthelinesofprincipalstresses[ 38 ].Second,geometricandmaterialpropertiesvaryalongthespanwisedirection,andtheyhavetobedescribedaccurately.Third,ifoneoftheprincipaltensionsvanishes,themembranewillbecompressedandwrinkled,whichfurthercomplicatestheanalysis.Consideringthecomplexityofthethree-dimensionalproblem,simplicationsaremadeintheearlyworks. OdenandSato[ 70 ]presentedaniteelementmethodfortheanalysisoflargedisplacementsandnitestrainsinelasticmembranesofgeneralshape.Intheirworkthestaticequilibriumofaninatedmembranewasconsidered.Astaticanalysiscanbeconsideredasaspecialcaseofadynamicanalysisinwhichtheinertiaeectisincluded.Worksonthedynamicanalysisofmembranesbefore1991werereviewedbyJenkinsandLeonard[ 41 ];mostrecentworkswerereviewedbyJenkins[ 40 ].Verronetal.[ 105 ]performedacombinednumericalandexperimentalstudyofthedynamicinationofarubber-likemembrane.Membranewrinkleproblemsweretackledin PAGE 46 theworkofDingetal.[ 14 ];theyperformedanumericalstudyofpartialwrinkledmembranes.ThestabilityissuesofauidowpastamembranewerereportedbyThaokarandKumaran[ 99 ]. Inthischapterathree-dimensionalmembranemodelisproposedforthemem-branedynamics.ThemembranematerialconsideredobeysthehyperelasticMooney-Rivlinmodel[ 64 ].WeusetheGreen-Lagrangestraintensorforthedescriptionoflargestrains.Thedynamicresponseofthemembraneisdescribedbyasystemofsecond-orderordinarydierentialequations.Aniteelementmethodwithtriangu-larelementsforspatialdiscretizationisutilized.Fortimeintegration,theimplicitWilson-methodisused. 5.2 StructuralSolver 5.2.1 Mooney-RivlinModel Foraninitiallyisotropicmembrane,GreenandAdkins[ 25 ]showedthatthereexistedastrainenergyfunctionWwhichwasexpressedasthefollowingform whereI1,I2,andI3aretherst,second,andthirdinvariantsoftheGreendeformationtensorC.Ifthematerialisincompressible,namely,I3=1,thenthestrainenergyisafunctionofI1andI2only.Itcanbewrittenas: wherec1andc2aretwomaterialparameters.AmaterialthatobeysEq.( 5.3 )isknownastheMooney-Rivlinmaterial[ 64 ],andthemodeliscalledtheMooney-Rivlinmodel.TheMooney-Rivlinmodelisoneofthemostfrequentlyemployedhyperelasticmodelsbecauseofitsmathematicalsimplicityandrelativelygoodaccuracyforreasonablylargestrains(lessthan150percent)[ 64 ]. PAGE 47 Foraninitiallyisotropicmembrane,thegeneralstress-strainrelationshipiswrit-tenas @C;(5.4) whereSisthesecondPiola-Kirchhostresstensor,pisthehydrostaticpressure.Ifthemembraneisincompressible,thenthestress-strainrelationcanbefurthersimpliedasthefollowing whereIis3-by-3identitymatrix.Inthemembranetheorythestresseldisessentiallytreatedastwo-dimensional,andthestressnormaltothedeformedmembranesurfaceisnegligibleincomparisonwiththestressesinthetangentplane[ 70 ].Underthisconsideration,thedeformationandstressmatrices,C(t)andS(t),aregivenby Therefore,thehydrostaticpressureisdeterminedbyconsideringtheconditionthatS33=0, where3isdenedby inwhichh(t)andh0arethemembranethicknessinthedeformedandundeformedcongurations,respectively. 5.2.2 PrincipleofVirtualWork Consideranelasticmembrane,whichinitsinitialstateoccupiesaniteregioninspace.ThepositionoftheregioncanbespeciedbyanorthogonalcoordinatesystemXYZwhichisreferredtoasglobalcoordinates.Inglobalcoordinates,the PAGE 48 generalformoftheprincipleofvirtualworkintheLagrangiandescriptionisgivenby whereV0and@V0are,respectively,thevolumeandtheboundarysurfaceoftheundeformedmembrane,0isthemembranedensity,D(t)isthedisplacementvectoringlobalcoordinates,D(t)isthevirtualofvectorD,thesuperscriptTreferstothetransverseofavectororatensor,D(t)istheaccelerationvector,P(t)isthegeneralizedexternalsurfaceforce,EistheGreen-Lagrangestraintensor,SisthesecondPiola-Kirchhostresstensor,andE:Sisthescalarproductoftwotensors.IntherectangularCartesiancoordinates,thescalarproductisdenedby withsummationontherepeatedindices. 5.2.3 TriangularMembraneElement Theprincipleofvirtualworkcanbewritteninthefollowingdiscretizedform whereisthetotalvirtualwork,andthesubscripterepresentselementaryquan-tities,Neisthenumberofelementscoveringtheentiremembrane.Inthepresentworkthetriangularelementisemployed.AsshowninFig. 5{1 ,atypicaltriangularelement,e,isdenedbynodes,i,j,k,andthreestraight-lineboundaries.Intheglobalcoordinatesthethreenodeshavecoordinates(Xi,Yi,Zi),(Xj,Yj,Zj),and(Xk,Yk,Zk),respectively. Tofacilitatethedevelopmentofacomputationalprocedure,itisconvenienttoconsiderthemembraneelementanditsdisplacementinthelocalcoordinates.AsshowninFig. 5{1 ,inlocalcoordinates,thetriangularelementliesintheplaneof PAGE 49 Schematicglobalandlocalcoordinates. whereXiistheglobalcoordinateofnodei,andTisanorthogonaltransformationmatrixwhichsatisesTT=T1.Similarlythedisplacementshavethefollowingrelation whereDisthevalueofadisplacementinglobalcoordinates,disitsvalueinlocalcoordinates. Iftheelementissucientlysmall,wecanapproximatethedisplacementatanypointwithintheelementwithalinearcombinationofthedisplacementsofnodes,i,j,andk, PAGE 50 whereN3isaninterpolationmatrixdenedas where01,01,and011.deisavectordenedby inwhichdi,dj,anddkarethenodaldisplacementvectorsofnodesi,j,andk,respectively.Eachnodehasthreedegreesoffreedom, 5.2.4 Green-LagrangeStrainTensor Whenadisplacementisintroducedthelocationofapointoriginallygivenbycoordinatesxnowhasnewcoordinates whereuisthedisplacementasshowninFig. 5{1 .FollowingGreenandAdkins[ 25 ],theGreen-Lagrangestraintensorisdenedbytherelation 2@yk Sincethethicknessofthemembraneissmallcomparedtotheothercharacteristicdimensions,thefollowingsimplicationsareallowedtobemade[ 25 70 ] 2@u 2(21);;=1;2;(5.20) PAGE 51 whereisthesameas3inEq.( 5.8 ).Tosimplifythecomputation,wedividetheGreen-Lagrangestraintensorintotwoparts.OnepartisE0,whichtakesintoaccountsmallstrains, @x1 2(@u @y+@v @x)01 2(@u @y+@v @x)@v @y00003777775;(5.21) andtheotherisEL,whichisthenonlinearpartofthestraintensor, @x@u @x+@v @x@v @x+@w @x@w @x@u @x@u @y+@v @x@v @y+@w @x@w @y0@u @x@u @y+@v @x@v @y+@w @x@w @y@u @y@u @y+@v @y@v @y+@w @y@w @y000213777775:(5.22) TondeachcomponentinthestraintensorsE0andEL,itisnecessarytocomputethedisplacementgradientur(x;y).Aftersomeoperationsweobtaintheexpressionofdisplacementgradient, @x@u @y@v @x@v @y@w @x@w @y3777775=26666664uiuj BasedontheexpressionsinEq.( 5.23 )wecancomputethecomponentsE0andEL. 5.2.5 InternalForces Whenlargedeformationsand/ornonlinearmaterialpropertiesareinvolved,theequivalentinternalelasticresistingforcesofthestructurecanbeexpressedasthefollowing, whereisthenonlinearstressvector,"isthestraintensor,andBisthestrain-displacementmatrixwhichisalsoafunctionofthedisplacementforlargedeforma-tion.Recallthatintheprincipleofvirtualworkthetermrelatedtotheinternal PAGE 52 forcesfromthediscretizedformofprincipleofvirtualworkEq.( 5.11 )reads Eq.( 5.6 )tellsusthatthemembranestressisessentiallytwo-dimensional,therefore,weknowfromEq.( 5.10 )thatitnecessarytoconsiderthereducedstraincomponentsinsteadofallthecomponentswhencalculatingtheinternalforcesbasedonEq.( 5.24 ).ThegeneralmembranestrainvectorcanbeconvenientlydecomposedintotwopartsaccordingtothestraindecompositioninEqs.( 5.21 )and( 5.22 ), where"0isthereducedstrainvectorforthelinearpartofthestraintensor,E0,and"LforthenonlinearpartEL.BasedonEqs.( 5.21 )and( 5.22 )thesetwostrainvectorscanbeexpressedasthefollowing, @x@v @y@u @y+@v @x9>>>>>=>>>>>;;"L=8>>>><>>>>:"Lx"LyLxy9>>>>=>>>>;=1 28>>>>><>>>>>:@u @x@u @x+@v @x@v @x+@w @x@w @x@u @y@u @y+@v @y@v @y+@w @y@w @y2(@u @x@u @y+@v @x@v @y+@w @x@w @y)9>>>>>=>>>>>;:(5.27) Oncethevirtualdisplacementsaregivenatanypointwithintheelement,thestrainvariationatanypointcanbedetermined.Therelationshipsbetweenthevirtualdisplacementandstrainvariationare whereB0andBLarethestrain-displacementmatrices,theirexpressionscanbefoundinAppendixA.BasedonEq.( 5.27 )thevirtualworkineachtriangularelementcan PAGE 53 becomputedas whereB=B0+BL,"T="T0+"TL,andFinteistheinternalforceinglobalcoordinates inwhichTisthesametransformationmatrixasshowninEq.( 5.13 ). 5.2.6 ExternalForces Asthemembranechangesitsshapeandpositionunderthepressure,notonlythedirectionofthepressurechangesbutalsotheareaonwhichitactschanges[ 70 105 ].Toaccountforthissituation,weassumethatthesizeofeachniteelementissucientlysmallthatthepressureisuniformovereachelement.Tofurthersimplifythecomputation,weagainuselocalcoordinates.AsshowninFig. 5{2 ,thelocalcoordinatesliesinthedeformedplaneinsteadofintheundeformedplane.InthelocalcoordinatesthepressurePisgivenby, withn=[001]Tistheexternalnormaldirectionoftheelement. ThevirtualworkbytheexternalforceateachelementisZ@V0eDTPdS; PAGE 54 Schematicoflocalcoordinatesforexternalforces. thenetexternalforceateachnodeisobtainedbysimplyrepresentingitbythreeequalforcesateachnode 3001 3001 3T=DTeFexte:(5.33) Hence,inglobalcoordinatesthenodeforceonelementeduetoauniformpressureoverthatelementis 3001 3001 3T:(5.34) 5.2.7 DynamicEquationsfortheMembrane Thersttermoftheleftsideoftheprincipleofvirtualworkequation( 5.11 )canbewrittenas PAGE 55 whereMeisthemassmatrixofanelement, 120s0h02666642IIII2IIII2I377775;(5.36) inwhichIisa3-by-3identitymatrix. Uptonowwehaveobtainedtheexpressionsfortheinternalforce,externalforce,andmassmatrixonanelemente.Itisnecessarytoassembleallthediscreteelementsintoasingleunit.Foranyelementetheprincipleofvirtualworkis Assemblingalltheforcesinoneelementwehavethefollowingequation, Thisequationisforasingleelement,ifwesumallelementsinglobalcoordinates,thenweobtainthesystemofgoverningequationsformembraneunderexternalloads, whereMisapositivedenitemassmatrix,whichremainsconstant,Fintistheinternalforce,andFextistheexternalload. 5.2.8 SolutionoftheDynamicEquations Eitheranimplicitoranexplicitschemecanbeadoptedtofollowthetimehistoryofthesystemofequations.Theexplicitschemeiscomputationallyeconomicforeachtimestepandrequireslessstorage.Itseciencycanbefurtherimprovedbyspecialtechniques,suchaslumpingtechniqueofthemassmatrix.However,itstimestepneedstobesmalltosatisfythestabilityrequirement.Asfarastheimplicitschemeisconcerned,thematrixsystemisupdatedmorethanoncepertimesteptoaccount PAGE 56 forthenonlineareect,whichrequiresmorecomputationspertimestepandmorestorage.Generallyspeaking,implicitalgorithmsareeectiveforstructuraldynamicsproblemsinwhichtheresponseiscontrolledbyasmallnumberoflowfrequencymodes,whileexplicitalgorithmsareecientforwavepropagationproblemsinwhichthecontributionofintermediateandhighfrequencystructuralmodestotheresponseisimportant[ 94 ]. Therearenumerouswaystointegratethestructuraldynamicsequationsintime;amongtheseweutilizetheWilson-method[ 113 ].TheWilson-methodhasacontrollablealgorithmicnumericaldissipation,itisaone-stepmethod,andityieldssecondorderaccuracyintime.Thesecharacteristicsmakeiteasiertocodeandbehavewellinouruid/structureinteractionsystem.SincetheWilson-methodisdissipative,weemployatimestepsmallerthanthatrequiredbytheoverallaccuracycriteriontoavoidexcessivedissipation.Thisapproach,with=1.4,workswellforthepresentproblems. 5.3 NumericalTestofMembraneModel Theinationofasphericalmembraneunderanimposedpressurestepwasan-alyzedexperimentallybyAlexander[ 1 ]andtheoreticallybyBeaty[ 4 ].Recently,anumericalstudywasperformedbyVerronetal.[ 105 ].Thecorrespondingsemi-analyticalsolutionwasreportedbyVerronetal.[ 104 ].Thesestudiesfocusonanincompressiblesphericalhyperelasticmembrane.ThemembranehasaninitialradiusofR0,anundeformedthicknessofh0,andamassdensityof0,respectively.ThemembranematerialisassumedtoobeytheMooney-Rivlinmodel.Therearetwocon-stitutiveparametersinthismodel,namelyc1andc2.Tosimplifyourcomputation,wedeneanon-dimensionalMooney-Rivlinconstantbasedonthesetwoconstitutiveparameters PAGE 57 TheanalyticalgoverningequationforthesphericalmembraneisdetailedbyVerronetal.[ 104 ].Hereweonlygivethedimensionlessform =p2+(1 whereisthecircumferentialprincipleextension,whichisdenedastheratioofthedeformedtotheundeformedradius.Thecorrespondingnon-dimensionaltime,andthenon-dimensionalimposedpressurestepparegivenasfollows: whereP0isthedimensionalpressurestep.Theinitialconditionsofthemembranemayvary;here,wefocusonthecasethatthemembraneisinitiallyinitsequilibriumstate,i.e.,withthefollowinginitialconditions: Onecanobtainthesemi-analyticalsolutiontoEqs.( 5.40 )and( 5.43 )usinganystan-dardsolver.Thedirectnumericalsolutionsofthemembraneequationareobtainedwiththeniteelementmethoddiscussedbefore.Thecomputationalmeshhas816nodesand1568triangularelementsonasemi-sphericalsurface.Inourworkwechooset=5104whichis500timeslargerthantheexplicitschemeusedbyVerronetal.[ 105 ]wheretheyuseanexplicitschemewithatimestepof106forthesameproblem. Allcomputationsareperformedusingthedouble-precision.InFig. 5{3 wecom-paretheniteelementresultswiththesemi-analyticalsolutions.ExceptforonecaseinFig. 5{3 b,allnumericalsolutionsmatchwellwiththesemi-analyticalsolutions.Thelossoflinearityofoscillationsiswellreproduced.Themismatchedcaseshownin PAGE 58 Fig. 5{3 bcorrespondsto=0.1.Forthisvalueofthematerialparameter,aspointedoutbyVerronetal[ 105 ],thereisanunstableequilibriumpointatp0.687.Neartheunstableequilibriumpoint,theniteelementresultscanonlyapproachqualitativelytherealbehavior.Thisisthereasonwhythenumericalresultsarehighlysensitivetotheassignedpressurestep.Asexpected,thelargertheMooney-Rivlinconstantis,thehigherpressurethemembranecansustain. TwoprofoundissuesshouldbeconsideredintheuseoftheWilson-method.First,theWilson-methodisunconditionalstableforlinearproblems.Thetimestepinnonlinearproblemsisoftendictatedbytheerrorsintroducedintheapproximationtreatmentofthenonlinearities;therefore,duecareshouldbetakentochooseanapproximatetimestep.Second,asmalltimestepisdesirablefortheWilson-methodduetoitsexcessivedissipation.AsshowninFig. 5{4 ,alargertimesteptendstocausemorenumericaldamping.InFig. 5{5 weshowtheerrornormoftheimplicitscheme,thegureexhibitsaslopearound1.Therefore,theoverallaccuracyfortimeintegrationisapproximatelyrstorder. Thepresentmembranemodelgivesagoodagreementwiththeanalyticalsolu-tionsformembranedynamicsunderlargedeformation.Itleavesoneissueoffunda-mentalimportancesinceitcannothandlewrinklephenomenonwhichhappenswhenthemembraneiscompressed.Themembranewingissometimescompressed,espe-ciallyattheleadingedge.Therefore,furtherinvestigationisnecessaryinthisarea. PAGE 59 Figure5{3: DynamicinationofaMooneysphericalmembrane:(a)=0;(b)=0:1;(c)=0:25.Solidline:Semi-analyticalsolution;Circle:Numericalsolution. PAGE 60 Theeectoftimesteptothenumericaldampingforp=0.5and=0. Errornormversustimestepfortheimplicitschemeforp=0.5and=0. PAGE 61 6.1 Introduction Membranestructuresarefoundinmanyengineeringpractices.Traditionally,theyareutilizedinparachutesandsailboatstoimproveperformance.Arecentinterestingapplicationisamembranewing-basedMicroAirVehicle(MAV).MAVswithamaximumdimensionof15cmorlessandightspeedaround10m/sprovideexpendableplatformsforsurveillanceanddatacollectioninsituationswherelargervehiclesarenotsuitableorelsedangeroustoreach.However,MAVsaresensitivetotheenvironmentaldisturbances,anditishardtoavoidconsideringtheirsmalldimensions,lightweights,andlowightspeeds.Ourexperienceindicatesthatanaeroelasticmembranewingcanbetteradapttotheatmosphericdisturbances,providesmootherviewingplatform,andmakethevehicleeasiertoy[ 83 ]. Duringight,thesurroundingowexertsaerodynamicforcesonthemembranewingandcausesstructuraldeformationsofthemembranesurface;meanwhile,thewingdeformationsaecttheowpatternandchangetheaerodynamicforcedistri-butions.Themutualinteractionsbetweentheuidowandthestructurecauseauidandstructureinteractionproblem.Fluidandstructureinteractionproblemshaveattractedinterestofmanyresearches.Numerousstudieshavebeenconductedonthesesubjects.Numericalschemesforuidandstructureinteractionproblemsmayvary.However,ingeneral,therearetwokeyelementsinacoupleduidandstructureinteractionsystem.Oneistheuidsolverthatevaluatestheaerodynamicforces,andtheotheristhestructuralsolverthatcomputestheshapedeformation.Besidethese,amovinggridtechniqueisalsoneededtoregeneratethecomputational 53 PAGE 62 grid;aninterfacealgorithmisutilizedtoexchangeinformationbetweentheuidandstructuralsolvers. Theinteractionbetweenauidandamembraneusuallycausessubstantialstruc-turaldeformationthatisinherentlynonlinear.Themembranebehaviorisingeneralcomplexandearlierstudiesofinteractionsbetweenuidowsandmembranestruc-turesarelargelybasedonempiricalobservations[ 34 ],two-dimensionalowbasedcomputations[ 36 86 92 93 ],orsimpliedthree-dimensionalanalyses[ 37 38 ].Inourstudywepresentacoupleduidandstructureinteractionsystemtoinvestigatetheinterplaybetweenamembraneandincompressibleviscousow.Inthecoupledsystem,theuiddynamicsisgovernedbythethree-dimensionalNavier-Stokesequa-tionsforincompressibleows,andthestructuraldynamicsisgovernedbyanonlineardynamicmembranemodel.TheNavier-Stokesequationsaresolvedusingapressure-basedmethodonastructuredandmulti-blockgrid.ThemembranematerialfollowsthehyperelasticMooney-Rivlinmodel[ 64 ].Themembranemodelisderivedfromtheprincipleofvirtualwork.ItisspatiallydiscretizedwithatriangularniteelementmethodandtemporallyintegratedwiththeWilson-method[ 113 ].Toaccommodatethemembraneshapechange,theowsolverisfurtheraugmentedwithamovinggridtechniquediscussedinChapter 3 Theinteractionsystemiscoupledthroughtheexchangeoftheaerodynamicforcesandshapedeformationsbetweentheuidandstructuresolvers.Sincetheuidsolverandstructuralsolversdonotsharetheidenticalgridattheinterface,thethinplatespline(TPS)interpolation[ 16 ]isusedtomaptheexternalloadandtheshapedeformationtothestructuralanduidsolvers,respectively.Toavoidphaselagbetweentheuidandstructuralsolvers,subiterationsareperformedtosynchronizetheuidandstructuralsolvers. Inthefollowingwerstintroducetheinterfacealgorithmbetweentheuidandstructuralsolvers,thisisfollowedbyastep-by-stepdescriptionoftheuidand PAGE 63 structureinteractionalgorithm.Thiscouplingisachievedthroughexchanginginfor-mationandsynchronizationbetweentheuidandstructuralsolvers.Intheaerody-namicanalysispart,westudyboththerigidandmembranewings.Intherigidwingstudy,wefocusonthecharacteristicsofthelowReynoldsnumberandlowaspectratiorigidwing,includingboundarylayerseparationandreattachment,tipvortices,andunsteadyphenomenonathighangleofattack.Inthemembranewingstudy,westudytheself-initiatedmembranevibration,itsimpactsonthemembranewingperformance,andwealsocomparethemembraneandrigidwingperformances. 6.2 CouplingbetweentheFluidandStructuralSolvers Whenauidowpassesaexiblestructure,thestructurechangesitsshapeduetotheexternalforcessituatedinow;meanwhile,theshapechangeaectstheuidowaroundthestructureandtheforcedistributionsonthestructure.Thisleadstoauidandstructureinteractionproblem.Usually,theuidproblemandthestructuralproblemhavedierentmathematicalandnumericalproperties.De-pendingontheproblemcomplexity(linearproblemornonlinearproblem,simpleorcomplexgeometry,small-scaleorlargescaleproblem),availablesoftware,andcom-putingplatform,numeroustechniqueshavebeendevelopedtoprovidecoupledCFDandCSDmethods.Oneapproachistotreattheowandstructuralequationsasonemonolithicsystemofequations[ 5 ]andsolvethesystemofequationssimultaneouslyona\monolithic"grid.However,sincemultipletimescalespresentinauidandstructuralsystem,auniedtreatmentmaynotbeecienttohandlethecomputa-tionalstiness.Liuetal.[ 58 ]solvedtheowandstructuraleldssimultaneouslybyafullyimplicitmethod.However,unlikethemonolithicmethod,theirCFDsolverandCSDsolveraresolvedondierentgridsystemsinwhichTheuidandstructuralgoverningequationsareregardedasonesinglesystemoftime-dependentequationsinapseudo-time. PAGE 64 Tobetterhandledierentcharacteristicsintheowandstructuraldomains,mostoftheresearcherssolvedtheuidandstructuralequationsseparatelyandthencou-pledthemthroughasynchronizationprocess.SmithandShyy[ 93 ],Kamakotietal.[ 49 ],GordnierandFithen[ 22 ],KalroandTezduyar[ 47 ],ZwaanandPrananta[ 115 ],andLianetal.[ 53 55 ]coupledtheuidandstructuralsolverswithasubiterationapproach.Theuidsolverandthestructuralsolverfunctionindependentlyontheirowncomputationalgridswiththeirowntimesteps,subiterationsbetweenthesetwosolversareemployedateachtimesteptoavoidthephaselag.Alternatively,Farhatetal.[ 26 28 ]formulatedtheuidandstructureinteractionproblemasathree-eldproblem:theuid,thestructure,andthedynamicmeshthatisrepresentedbyapseudo-structuralsystem.Intheirwork,thearbitraryLagrangianEulerian(ALE)nitevolumeschemeisemployedfortheuiddynamicsequations,aniteelementmethodisappliedforthestructuraldynamicsequations.Partitionedproceduresandstaggeredalgorithmsarethenadoptedforthecoupleduidandstructureinteractionproblemsinthetimedomain. Inuid-structureinteractioncomputations,amovinggridtechniquecanbefruit-fullyemployedtoautomaticallyregeneratethenewCFDgridaccordingforthede-formedconguration.Inaddition,sincetheCFDsolverandCSDsolverdonotnecessarilyshareidenticalgridattheinterface,interpolationsareneededtoallowtheinformationexchange.Smithetal.[ 91 ]evaluatedsuitablemethodstotransferinformationbetweenCFDandCSDgrids.Inourstudyathinplateinterpolationmethod[ 16 ]isadoptedforthispurpose.Thethinplatesplineinterpolationisaglobalinterpolationmethod,itisinvariantwithrespecttorotationsandtranslations,andhenceissuitablefortheinterpolationonmovingsurfaces.Aone-dimensionalversionoftheinterpolationis PAGE 65 whereH(x)isthedisplacementdistributionfunctionoverthedomain,iistheunde-terminedcoecient,Nisthenumberofmonitoredstructuralnodesontheinterface,andxi'saretheirlocations.ItcanbeseenfromEq.( 6.1 )thattheinterpolatedfunctionH(x)hasacontinuousrstorderderivative.Oncethesecoecientsiaredeterminedbasedonthestructuralgridlocations,thedisplacementvectordenedontheCSDgridDismappedtotheCFDgridviaamatrixG whereXisthecorrespondingdisplacementvectorontheCFDgrid.SimilarlythesurfaceforcescanalsobemappedfromtheCFDgridtotheCSDgrid. Uponupdatingtheaerodynamicforcesinthestructuralequationsandprovidingthenewboundaryconditionstotheuidsolver,thetemporaldevelopmentbetweentheaerodynamicandstructuralequationscanbecoordinated.Toconducttime-accuratecomputationsfortheuid-membranedynamics,iterationsbetweentheuidandstructuralsolversshouldbedoneateachtimeinstant.Astep-by-stepcoupledalgorithmisdescribedbelow: 1. GeneratetheCFDandCSDgrids. 2. Setn=0,evaluatepressureP0ontheundeformedwingbasedontheinitialvalues0,u0,and0. 3. n=n+1,andtn=tn1+t. 4. SolvetheNavier-StokesequationsandtransferthepressurePn+1tothestruc-turebasedontheTPSinterpolation. 5. UsethestructuralsolvertoevaluatethedisplacementvectorDn+1. 6. Performsubiterationstosynchronizetheuidandstructuralsolvers,andseti=0. (a) i=i+1; PAGE 66 (b) MapthedisplacementsfromtheCSDgridtotheCSDgridwiththeTPSinterpolationXin+1=GDin+1. (c) UpdatetheCFDgridwiththemovinggridtechnique. (d) UpdatetheJacobiantosatisfythegeometricconservationlaw. (e) Updatetheboundaryconditionsfortheuidsolver. (f) ComputethepressurePin+1eldbasedonthenewCFDgridandboundaryconditions. (g) MaptheaerodynamicloadfromtheCFDgridtotheCSDgridwiththeTPSinterpolation. (h) CalculatethestructuraldisplacementDi+1n+1. (i) Checkforconvergenceforthesubiterationprocess. If`No',returntostepa),otherwise,continue. 7. Returntostep3toprocessnexttimestep. Itshouldbenotedthatafullyconvergedsolutionintheprocessishardtoachieve.Thenumberofsubiterationisbasedonthecomputationalresource,problemnonlinearity,andthetimestep.Inourstudywendtwoorthreesubiterationsareenoughtoalleviatethephaselagerrors. 6.3 MAVWingAerodynamicsandStructuralResponse AsshowninFig. 1{1 themembranewingconsistsofaleadingedgesparandthreechordwisebattensoneachsideofthehalfwing.Themembranematerialisbondedtothesparandbattens.Whilethemembraneisexibleanddoesnotsustainbendingmoment,thecarbonberisrigidwhichprovidessupportforthemembraneandcansustainbendingmoment.Tofullymodelthemembranewing,ideally,oneneedstomodelboththerigidbattensandtheexiblemembrane.Twoprincipaldicultiesarisebydoingthat.First,themembranedoesnotbearbendingmoment,ithasthreedegreeoffreedomsateachnode;thebatten,ifmodeledasabeam,hasvedegreesoffreedomateachnode.Aspecialtreatmentisrequiredtotreatthe PAGE 67 interfacepointwhichbelongstodierentregimesandhasdierentdegreesoffreedom.Second,themembraneismuchthinnerthanthebatten,themassmatrixassociatedwiththesetwodierentmaterialsmakestheassembledmassmatrixill-conditioned.Giventhesefactors,wetreatthebattenasaspecialmembranematerialwithalargerdensitywhiletheotherpropertiesarethesame;thedensityratiobetweenthebattenandthemembraneissettobethree. 6.3.1 ComputationalBackground Inourstudyonlythemembranewingisconsideredwhileignoringthefuselageandpropellers.AschematicgeometryofthewingisshowninFig. 6{1 .Theshadedareashowncorrespondstothefuselagewhichisxed.Thephysicalproblemunderconsiderationisviscousowoverawinginanunboundeddomain.Basedonthefreestreamvelocityof10m/s,therootchordReynoldsnumberis9104.Computa-tionsarerstperformedtoassessthesensitivityofthecomputedresultstothewinglocationandthetypesofboundaryconditionsimposedtotheouterowboundaries.Basically,theboundariesshouldbesetsucientlyfarawayfromthemembranesothattheydonothavemuchimpactonthecomputedaerodynamiccoecients.BasedonourtestsaproperwingpositionischosenandshowninFig. 6{2 .Theupstreaminletboundaryis7cfromtheleadingedge,herecisthechordlengthattheroot.AtsurfaceABCDazerogradientboundaryconditionisimposedforthevelocitycompo-nents;alltheotherboundariesareassignedasaDirichlettypeboundarycondition.Sincethefreestreamvelocityisparalleltothechordwisedirectionandnopropellerismodelled,onlyhalfofthedomainiscomputed,andasymmetricmappingcanbeappliedtotheotherhalfdomain. Gridsensitivitytestsarealsoperformed.Threegridsystemshavebeensystem-aticallychosen,rangingfrom1:8105nodesatthecoarseleveland2:3106nodesatthenelevel.Forthecoarsegridthereare41pointsinthechordwisedirectionand31pointsinthespanwisedirectiononthewingsurface.Thegriddistribution PAGE 68 Figure6{1: Schematicwingwiththreebattensoneachhalfwing. Figure6{2: Computationalset-upandboundaryconditions. PAGE 69 atthecoarselevelaroundthewingisshowninFig. 6{3 .Theresultsfromthenegridareusedasthereference.ThecomparisonamongthesethreegridsystemsissummarizedinTable 6{1 .Thecoarsegridsolutionunder-predictsthelift-to-dragratio(CL=CD)by0.5%,andthedierenceinformdragislessthan0.7%.Asimilarconclusioncanbedrawnforthelift,L.Asexpected,theskinfrictiondrag,DFissensitivetothegriddensity;thecoarsemeshoverestimatesDFbyalmost10%.Anintermediategrid,whichhas6:7105nodesandmorethandoublesthesurfacepointsonthewingsurfacecomparingtothecoarsegrid,givesanimprovedpredictiononskinfrictiondrag.FromTable 6{1 itcanbeseenthatskinfrictiondraghasmuchlesscontributiontothetotaldragthanformdrag.Thisisalsovalidathighanglesofattackwheretheformdragisthedominantfactorinthetotaldrag. Figure6{3: Schematicwinggeometryandstructuredgridwith1.8105pointsand4131onhalfwingsurface. Wefurthercomparethemainowfeatures.Flowseparationsareobservedwithallthesecomputationalgridsatthelowerwingsurface.Weextractthechordwisevelocityproleatthequarterchordoftherootregion.Theboundarylayervelocity PAGE 70 Table6{1: Gridrenementeectsonthecomputedaerodynamiccoecients. GRID 5.10E-1 9.70E-3 7.06 Intermediate6.7105 5.08E-1 9.22E-3 7.16 Fine2.3106 5.04E-1 8.80E-3 7.10 proleisshowninFig. 6{4 .Thenegridhasthelargestreversevelocity,thecoarsegridhasthesmallest,andtheintermediategridhasthethickestseparationzone.Theseresultsindicatethatthegridresolutionnoticeablyaectsthesharpnessofthevelocityprole,aswellastheextentoftheowseparation.Nevertheless,thecoarsegridpredictsthemainfeaturesoftheowwhichisqualitativelyconsistentwiththoseonthenegrid. Figure6{4: Spanwisevelocityprolesatthebottomsurfaceatthequarterchordoftherootregionforrigidwingat=6o. Comparisonsbetweentheintermediateandcoarsemeshesatotheranglesofattackarealsoconducted.ThecomputedaerodynamiccoecientsareshowninFig. 6{5 .Boththeliftanddragcoecientsshowgoodagreements,withalittlebitlargerdierenceathigheranglesofattack.Afewmorewordsshouldbeputforthecomputationalresultsathighanglesofattack.Generallyspeaking,vortexshedding PAGE 71 andboundarylayerseparationoccuratlargeanglesofattack,andtheyintroduceunsteadinesstotheaerodynamicperformanceandthereforeaecttheaerodynamiccoecients.Cummingsetal.[ 10 ]reportedthatatlargeanglesofattacktheunsteadycomputationspredictedmuchlowerliftcoecientsthanthesteadycomputations.Weperformsteadycomputationswhentheincidenceislessthan15owhileconductingbothsteadyandunsteadycomputationswhentheangleofattackislargerthanthat.Athighanglesofattack,trulysteadysolutionscannotbeobtained.The\steady"solutionsareobtainedwithsteadycomputationsoncecertainorderisreducedintheresidual.ThePISOmethoddiscussedinChapter 2 isusedfortheunsteadycom-putation.Wesetthetimesteptobe5105s.However,forourcomputations,thedierencebetweenthesteadyandthetime-averagedliftanddragcoecients,arefoundtobelessthan1%overallthestudiedcases.Fig. 6{5 showsthecomputedaerodynamiccoecientsfortherigidwingwithbothcoarseandintermediatemeshes.Atlowanglesofattack,thecoarsegridpredictsslightlylargerliftthantheinterme-diategrid;thetrendisreversedathighanglesofattack.However,bothpredictthatthestalloccursapproximatelyat51o. Aerodynamiccoecientsversustheangleofattack:(left)liftcoecient;(right)dragcoecient. PAGE 72 Foraproblemascomplicatedasthepresentone,trulygridindependentresultsarediculttoattain.Consideringtheavailablecomputingresourceandtheout-comeoftheaforementionedinvestigation,weusetheintermediategridfortherigidwingsimulationwhileusingthecoarsegridtoillustratethekeyissuesintheuidandexiblestructureinteractionsystem,todemonstratethesalientfeaturesofmem-branewingcomputations,andtoillustratetherelatedaerodynamicandstructuralcharacteristics. 6.3.2 RigidWingAerodynamics VorticalFlowStructuresoftheRigidWing Tipvorticesexistonanitewingduetothepressuredierencebetweentheupperandlowerwingsurfaces.Thetipvortexestablishesacirculatorymotionthatpresentsoverthewingsurfacesandexertsgreatinuenceonwingaerodynamics.Forexample,tipvorticesincreasesthedrag.Thetotaldragcoecientonanitewingatsubsonicspeedscanbewrittenas[ 2 ] whereCD;Pistheformdragcoecientduetothepressure,CD;Fisthefrictiondragcoecientduetoviscosity,eisthespaneciencyfactorwhichislessthanone,ARistheaspectratio,andCD;i=C2L 6.3 )demon-stratesthatinduceddragvariesasthesquareoftheliftcoecient.Furthermore,Eq.( 6.3 )illustratesthatasARisdecreased,induceddragincreases.TheMAVwinginourstudyhasalowaspectratioof1.4;therefore,itisimportanttoinvestigatethetipvortexeectsonthewingaerodynamics.Besidestheireectsonthedrag,tipvorticeshavetwo-foldedimpactsonthelift: PAGE 73 2 ]. 67 ]. Fig. 6{6 showstipvorticesaroundthewingsurfaceatdierentchordwisepo-sitionstogetherwiththestreamlinesattheangleofattackof39o[ 53 ].Thehigherpressurefromthelowersurfacedrivestheowtowardtheupperwingsurfacewherethethepressureislower.Thetwovorticesonbothsidesofthewingdemonstrateclockwiseandanti-clockwiserotations;itseemsthattheowisleadingfromthelowersurfacetotheuppersurface.Fig. 6{7 showsthelowpressurezoneduetothevorticalow.Thepressuredropfurtherstrengthenstheswirlbyattractingmoreu-idstowardthevortexcore;meanwhile,thepressuredecreasescorrespondinglyinthevortexcore.Thelowpressureregioncreatedbythevortexgeneratesadditionallift.Towarddownstream,thepressurerecoverstoitsambientvalue.Thentheswirlingmotionweakensandthediameterofthevortexcoreincreases.Thevortexcorelossesitscoherentstructureatdownstream. Fig. 6{8 visualizestheevolutionofthevorticalowstructurewiththeangleofattack.Thepressuredistributionontheuppersurfaceisalsopresentedinthesamegure.At=6otipvorticesareclearlyvisibleeventhoughtheyonlycoverasmallareaandareofmodeststrength.Theowisattachedtotheuppersurfaceandfollowsthechordwisedirection.Alowpressureregionnearthetip,whichisindicatedwithgreencolor,isobserved.Vorticesstrengthenwiththeincreaseofangleofattack.At=27o,asshowninFig. 6{8 ,tipvorticesdevelopastrongswirlingmotionandentrainthesurroundingowtothevortexcore.Thelowpressureareaincreasesastheangleofattackbecomeshigher. ThepressuredropinthevortexcoreisdemonstratedinFig. 6{9 awherethespanwisepressurecoecientontheupperwingsurfaceatx=c=0.4isplotted.At PAGE 74 Figure6{6: Streamlineandvorticesforrigidwingat=39o.Thevorticalstructuresareshownonselectedplanes. Figure6{7: Pressuredistributionaroundtherigidwinginthecrosssectionswithstreamlinesatangleofattackof39o. PAGE 75 Evolutionofowpatternforrigidwingversusanglesofattack. PAGE 76 6{8 ).Theliftbythevorticesisnotenoughtomaintaintheincreaseoftheliftcoecientwiththeangleofattack,andstalloccurs.Eventhoughttipvorticesaectthepressuredistributionatthelowersurface,theeectedregionismainlylocatednearthetip.Mostregionsinthelowersurfacearenotaectedbythetipvortices.FromFig. 6{9 bitcanbeseenthatthepressureisalmostuniforminthespanwisedirectionevenathighanglesofattack.Fig. 6{9 isillustrativeinregardtothepressuredistributionsversusthevorticalstructures,however,theyarenotindicativeofthetotallevelofthepressureforce.OneshouldalsonotethatFig. 6{9 indicatesthatthemomentexperiencessubstantialvariationsastheangleofattackchanges. TheLaminarBoundaryLayerSeparationonRigidWing Thelaminarboundarylayer,especiallyunderlowReynoldsnumberconditions,ispronetoseparateunderanadversepressuregradient.Thisseparationisusuallyfollowedbyaquicktransitionfromthelaminartoturbulentow.Iftheadversepressuregradientismodest,theresultingturbulentow,byobtainingenergythroughentrainment,mayreattachtothesurfacetoformalaminarseparationbubbleandthenformsanattachedturbulentboundarylayer. PAGE 77 (a)cpatuppersurface (b)cpatlowersurface Figure6{9: Spanwisepressurecoecientdistributionsatx/c=0.4forrigidwingatdierentanglesofattack. EversinceitsrstobservationbyJones[ 42 ],manyexperimentalinvestigationsonthelaminarseparationbubblehavebeenconducted,asreviewedbyYoungandHorton[ 114 ].Typically,alaminarseparationbubblehasthefollowingcharacteris-tics.Justdownstreamoftheseparationpoint,thereisa\dead-air"region,wheretherecalculatingvelocityissignicantlylessthanthefreestreamvelocityandtheowisvirtuallystationary.Theseparatedowformsafree-shearlayer,whichishighlyunstable.Theseparatedshearlayerquicklyundergoestransitionfromlaminartoturbulentow.Theturbulentowmayreattachtothesurfacebehindavorticalstructureandformsaturbulentboundarylayer.Hence,aseparationbubbleforms.ThedynamicsofalaminarseparationbubbledependsontheReynoldsnumber,thepressuredistribution,thegeometry,thesurfaceroughness,andthefreestreamtur-bulenceintensity.AroughruletodeterminethereattachmentpointisgivenbyCarmichael[ 7 ]thatsaystheReynoldsnumberbasedonthefreestreamvelocityandthedistancefromtheseparationpointtothereattachmentpointisapproximately5104.Ontheonehand,iftheReynoldsnumberislessthan5104,anairfoilwillexperienceseparationwithoutreattachment;ontheotherhand,alongseparation PAGE 78 bubblewilloccuriftheReynoldsnumberisslightlyhigherthan5104.Inthefol-lowing,detailedowstructuresaroundarigidwingarepresented.Bycontrastingthedetaileduidowaroundtherigidandmembranewings,onecanbetterevaluatetheMAVtechnologies.Thebubblestructuresatdierentanglesofattackaredemon-stratedinFig. 6{10 .Representativeinstantaneousstreamlinesattherootsectionontherigidwingareplotted.Atthismoment,thespanwisevelocitycomponentisignoredtomakeaclearrepresentationoftheseseparatedstructures.Atalowangleofattackof6o,theowprimarilyattachestothesurfaceandatinyseparationbubbleisobservedonthelowerwingsurfaceneartheleadingedge.Beginningat45%chordfromtheleadingedge,aweakrecirculationzoneisseenontheupperwingsurface,whichproducesareattachmentlengthofx/c=0.5. Withtheincreaseofincidence,theseparationpointmovesforwardtowardtheleadingedge.Attheangleofattackof15o,asshowninFig. 6{10 ,theseparationoc-cursat39%ofthechordfromtheleadingedge,whichisfollowedbyalong\dead-air"zonebeforeitreachesamaximalreversevelocityof0.26U.Inthe\dead-air"zone,thestationaryreverseowdoesnothavemucheectonthepressuredistribution,whichisprimarydeterminedbythewingcurvature.Theshearlayerreattachmenthappensat89%ofthechord.Basedonthefreestreamvelocityandthedistancebetweentheseparationandreattachmentpoints,theReynoldsnumberisapproximately4.5104whichisslightlysmallerthantheReynoldsnumbersuggestedbyCarmichael.Thereattachmentpointcorrespondstoarapidincreaseinthesurfacepressure,itishighlyunstable,andmovesforwardandbackward. ThestreamwisevelocitycontoursatdierentanglesofattackaredemonstratedinFig. 6{11 .Thevelocityisnormalizedbythefreestreamvelocity.At=6o,themaximumreversevelocityislessthan0.5%ofthefreestreamvelocity.Undersuchacondition,experimentshowsthataconsiderableportionoftheshearlayerremainslaminar[ 21 ]. PAGE 79 Streamlinesatdierentanglesofattackforrigidwing.Shownaresinglespanshortsoftheindividualtimedependentows. PAGE 80 However,withtheincreaseofangleofattack,themagnitudeofthereversevelocityincreasesaswell.At=27o,thereverseowcomponentoftheseparationbubblereachesamaximalvelocityof0.49U.Undersuchasituation,CromptonandBarrett[ 9 ]showedthattheshearlayerisparticularlyenergeticsincemostoftheshearlayeristurbulent.Eventhoughtheowontheuppersurfaceneartheroothasseparatedfromalargeportionofthesurface;however,theliftstillincreaseswiththeangleofattack.Tworeasonsmightliebehindthis.First,tipvorticesgeneratesuctionareasontheupperwingsurfacewhichprovideadditionallift.Second,eventhoughowseparatesneartheroot,itstillattachestotheotherregionofthewingsurface.Whentheangleofattackisincreasedfurther,at=51o,massiveseparationappears,andstalloccurs.Meanwhile,asseenfromFig. 6{12 ,vortexsheddingisobservedneartheroot.Atinyseparationbubbleemergesattheleadingedgerst,itthenobtainsenergyfromtheshearlayerandincreasesitssize.Thisbubblethenmergeswithanotherbubbledownstreamtoformalargerbubble.Thelargerbubbleisnotstable,itbreaksintotwosmallbubblesnearthemaximalcamberposition,andvortexsheddingbegins. Forlowaspectratiowings,tipvorticeshaveconsiderablecontributionstothelift.Thisscenarioisquitesimilartothatofdeltawings[ 23 ].Wehaveobservedthatthelowaspectratiowingsuerslessfromtheseparation.AsseenfromFig. 6{5 a,theliftkeepsalinearrelationshipwiththeangleofattackevenatveryhighanglesofattack.ExperimentsbyTorresandMueller[ 97 ]onlowaspectratiowingshavesimilarndings. PAGE 81 Normalizedchordwisevelocityu/Ucontoursatdierentanglesofattack. PAGE 82 Vortexsheddingat=51o. PAGE 83 UnsteadyPhenomenonatHighAnglesofAttack Bothvortexsheddingandboundarylayerseparationoccuratlargeanglesofattack,andtheyintroduceunsteadinesstotheaerodynamicperformance.Asmen-tionedbefore,weperformsteadycomputationswhentheincidenceislessthan15owhileconductingbothsteadyandunsteadycomputationswhentheangleofattackislargerthanthat.Interestingly,thedierencebetweenthesteadystatecomputa-tionsandthetime-averagedunsteadycomputationsaresmallevenatlargeanglesofattackinwhichunsteadyphenomenonsuchasvortexsheddingareprominent.Thedierencesintheliftcoecientandthelift-to-dragratioarefoundtobelessthan1%(Fig. 6{13 ).Noticeabledierenceappearsonlywhentheangleofattackbecomessucientlyhigh. ComparisonoftherigidwingCLbetweenthesteadyandunsteadycomputations. Thepressuredistributionsarealsocomparedattherootsection,whereowseparationusuallyappearsrstduetothelargecamberinthepresentMAVwingdesign[ 34 ].Fig. 6{14 showsthatat=6othetimeaveragedpressurecoecient PAGE 84 Figure6{14: Thecpcomparisonontherigidwingattherootbetweenthesteadyandunsteadycomputations:(left)=6o;(right)=15o. matchescloselywiththesteadystateresult.ThisndingisconsistentwithFig. 6{10 andFig. 6{11 whereaveryweakrecirculationisseenatthisangleofattack.However,at=15o,cleardierenceisshownafterx/c=0.6whichapproximatelycorrespondstothelocationofthemaximalreversevelocityinFig. 6{11 .Thetimeaveragedvalueyieldsasmoothpressuredistribution;thevariationinthesteadyresultisapparentfromtherecirculationzone.Velocitydistributionsshowthesametrendasthepressuredistributions.At=6o,asseenfromFig. 6{15 ,thereisalmostnodierenceinthechordwisevelocitydistribution.However,at=15o,nodierenceneartheleadingedgewherenoseparationoccurs,butcleardierenceisshownaftertheseparationbubble(Fig. 6{16 ). 6.3.3 MembraneWingDynamics TimeSynchronization Beforepresentingthemembranewingresults,somekeyissuesincoupledmembrane-uidinteractionsarerstaddressed.Specically,thetimesynchronizationbetweentheuidandstructuralsolutions,andthegeometricconservationlawinregardtothemovinggridmethodareemphasized. PAGE 85 Figure6{15: Comparisonofchordwisevelocityprolesontherigidwingbetweenthesteadyandtime-averagedunsteadycomputationsattheangleofattackof6o.Velocityprolesareattwochordwiselocations,andbothattherootsection. Figure6{16: Comparisonofthechordwisevelocityprolesontherigidwingbe-tweenthesteadyandtime-averagedunsteadycomputationsatangleofattackof15o.Velocityprolesatattwochordwiselocations,andbothattherootsection. PAGE 86 Thetimesynchronizationinourstudymeansthatiterationsshouldbeenforcedbetweentheuidandstructuralsolversateachtimesteptoavoidphaselagerrors.Weshallsimulatetheuidandstructureinteractionbyintegratingtheuidandthestructuralsolverswithsubiteration.Todemonstratetheimportanceofthetimesynchronizationbetweentheuidandstructuralsolvers,WecomparetheresultswithandwithoutsynchronizationinFig. 6{17 thatdepictsthedisplacementhistoryofatrailingedgenode.Thecomputationwithoutsynchronizationfailstocontinuewhilethatwithsynchronizationkeepsongoing.Thedisplacementhistorywithoutsynchronizationmatcheswellwiththatwithsynchronizationattheveryearlystage.However,thelaggingerrorsaccumulatewithtimeandeventuallyresultinamuchlargerdisplacementthanthatwithsynchronization.GordnierandVisbal[ 23 ]havealsoreportedtheimportanceofthetimesynchronization.ThevelocityandpressurehistoriesofthesamenodewithandwithouttimesynchronizationarecomparedinFig. 6{18 andFig. 6{19 ,respectively.Beforethecomputationdiverges,thevelocitywithoutsynchronizationismorethantwotimeslargerthanthatwithsynchronization.Thelargervelocityofthemembranenodemeansthatthemembranenodeobtainingmoreenergyfromthesurroundingowthanthatwithsynchronization.ThehighvelocityisbelievedtobeassociatedwiththesuddenpressuredropshowninFig. 6{19 ,whichcausesthewrinkleofthemembraneandthefailureofthestructuralsolver. TheGeometricConservationLaw Thegeometricconservationlaw[ 100 ]isanotherimportantfactorwhenthemov-inggridtechniqueisemployed.Fig. 6{20 ashowsthetimehistoryofCL/CDforthemembranewingat=6o.Withoutsatisfyingthegeometricconservationlaw,irreg-ularspikesareobservedinthecourseofcomputations.However,thecomputationsarebetterregulatedoncethegeometricconservationlawisenforcedasinFig. 6{20 b.Farhatetal.[ 27 ]arguedthatthegeometricconservationlawwasanecessarycon-ditiontomaintainthestabilityofaschemeutilizedinmovingboundaryproblems. PAGE 87 Figure6{17: Eectsoftimesynchronizationonthedisplacementofatrailingedgepointformembranewingat=6o. Figure6{18: Eectsoftimesynchronizationonverticalvelocityofatrailingedgepointformembranewingat=6o. PAGE 88 Figure6{19: Eectsoftimesynchronizationonpressuredistributiononatrailingedgepointformembranewingat=6o. Wedonotcomeacrossinstabilityinourstudy.Twopossibilitiesliebehindthat:oneisthesmalltimestepweusefortheuidsolver,theanotheristherelativelysmalldisplacementinthemembranewing.Nevertheless,thecomputationsbehaveerraticallywhenthegeometricconservationlawisnotenforced. Self-initiatedMembraneVibration Earlierworksinmembranewings[ 88 92 ]werebasedonconsiderationsthatthemembranewasmasslessandtherewasnotimedependentmovementinthesteadyfreestream.Experiments[ 111 ]showedthatthemembraneexperiencedhighfrequencyvibration.Byadoptingadynamicmembranemodel,ourcomputationsshowthatthemembranewingalsoexertshighfrequencyvibrationsinthesteadyfreestream.Fig. 6{21 demonstratesthedisplacementhistoryofthetrailingedgewithtime.Themaximaldisplacementduringthatperiodoccursnearthetip.Thedisplacementisnormalizedbythemaximalcamberthatisabout0.9cm.Toinvestigatethevibrationfrequency,wechooseapointbetweenthebatten1andbatten2onthetrailingedge PAGE 89 Figure6{20: EectofthegeometricconservationlawonCL/CDformembranewingcomputationat=6o:(left)notsatisfyingthegeometricconservationlaw;(right)satisfyingthegeometricconservationlaw asanexample.ItsdeectionhistoryisshowninFig. 6{17 ,andtheestimateddom-inatedfrequencyisaround120Hz(Fig. 6{22 ).Analysesatotherpointsalsoshowthefrequenciesaround120Hzthatiscomparabletotheexperimentalfrequencyof140Hz.[ 111 ]ConsideringthattheexperimentalconditionsandthedetailedwingcongurationbetweenourstudyandthosereportedbyWaszaketal.[ 111 ]arenotidentical,thequalitativeagreementbetweencomputationsandexperimentsinthisregardisdeemedsatisfactory.Liu[ 60 ]reportedthatunderatypicalwinggustsitua-tiontheenergywasmainlylocatedinthelowfrequencyrangeofseveralHzorlower.Themembraneuctuationisnotexpectedtocausesensitiveresponsetothevehiclesincethemembraneuctuatesinatimescalemuchfasterthaneitherthevehiclecontrolscaleortheexpectedwinggusttimescale. Whenthemembranesurfacevibrates,thevelocityonthewingsurfaceisnolongerzero.Aconsiderablevelocitycomponentisobservedatthewingsurface.TheverticalvelocitycontoursattwodierenttimeinstantsareplottedinFig. 6{23 .Fig. 6{23 acorrespondstoastagewhenthemembranewingmovesupundertheaerodynamicforces.Atthetrailingedge,theverticalvelocityisabout2%ofthe PAGE 90 Figure6{21: Timehistoryoftrailingedgedisplacementformembranewingat=6o.Thecamberattherootis0.90cm. Figure6{22: Spectrumanalysisofthetrailingedgepointvibrationformembranewingat=6o. PAGE 91 Figure6{23: Normalizedverticalvelocitycontourattwodierenttimeinstantsformembranewingat=6oontheuppersurface. freestreamvelocity.Anegativevelocitycomponentisalsoobservedneartheleadingedge.Fig. 6{23 bapproximatelycorrespondstoastageofmaximaldisplacement.Anegativevelocityneartherootindicatesthatthemembranebeginstomovedown. ComparisonbetweentheMembraneandRigidWings Withthesamefreestreamcondition,asshowninTable 6{2 ,theexiblewingexhibitsslightlylessliftcoecientthantherigidoneat=6o.ThedierenceinCL/CDisalsosmall.Athigherangleofattackof15o,themembranewinggeneratesaliftcoecientabout2%lessthantherigidwingdoes;however,itsCL/CDis1.5%morethantherigidone.Thesedierencesliebehindthehighfrequencyofthemembranevibration.Undertheexternalforce,themembranewingchangesitsshape.Thisshapechangehastwoeects.Ontheonehand,itdecreasestheliftbyreducingtheeectiveangleofattackofthemembranewing;ontheotherhand,itincreasestheliftbyincreasingthecamber.Ournumericalndingsshowthatmembraneandrigidwingsexhibitcomparableaerodynamicperformancebeforestalllimit,whichwasalsoexperimentallyobservedbyWaszaketal.[ 111 ]. Fig. 6{24 showsthetimeaveragedverticaldisplacementsofthetrailingedgeat=6oand=15o.Thedisplacementsarenormalizedbythemaximalcamberofthe PAGE 92 Table6{2: Aerodynamiccomparisonbetweenmembraneandrigidwings. MembraneWing RigidWing 6o 7.06 0.532 15o 3.88 0.954 wing.At=6othedeectionisabout15%ofthemaximalcamberandincreasestomorethan20%at=15o.Duetothetrailingedgedeectiontheeectiveangleofattackofthemembranewingislessthanthatoftherigidwing.ThespanwiseanglesofattackbetweenrigidandmembranewingsunderthesameowconditionandwithidenticalinitialgeometriccongurationsareshowninFig. 6{25 .AswecanseefromFig. 6{25 a,therigidwinghasanincidenceof6oattherootanditmonotonicallyincreasesto9:5oatthetip;themembranewingsharesthesameanglesofattackwiththerigidwinginthe36%oftheinnerwing;however,theeectiveangleofattacktowardthetipislessthanthatoftherigidwing.Atthetip,theangleofattackofthemembranewinglowersbyabout0.8o.Fig. 6{25 bcomparestheangleofattackat=15o,anditshowsthattheeectiveangleofattackofthemembranewingisabout1olessthanthatoftherigidwingatthetip.Thereducedeectiveangleofattackcausesthedecreaseinthelift. Themembranewingshapechangehasanothereect.Thepressuredierencebetweentheupperandlowerwingsurfacesinatesthemembranewing'scamber,whichisshowninFig. 6{26 .Twoairfoilshapesatdierentspanwisepositionsareplottedtogetherwiththecorrespondingrigidwingshapes.Thecamberincreaseismorevisibleintheouterwingthanthatintheinnerwing.ThisisconsistentwithFig. 6{24 wherethemaximaldisplacementoccursnearthetip.Asexpected,theincreaseislargerat=15othanthatat=6o. Thewingsurfacemovementaectstheoverallpressuredistribution.Fig. 6{27 comparesthetimeaveragedpressurecontoursbetweentherigidandmembrane PAGE 93 Figure6{24: Averageddisplacementofthemembranewingtrailingedge:(left)=6o;(right)=15o. Figure6{25: Comparisonofthespanwiseanglesofattackbetweentherigidandmembranewings:(left)=6o;(right)=15o. wings.Asexpected,thedierencesaremostlylocatedatthetrailingedgewherelargemovementisobserved.ThesedierencesarealsoreectedinthechordwisepressuredistributionsshowninFig. 6{28 .Attherootsectionthereisalmostnodierencebetweentherigidandmembranewingssincethefuselageisxed.How-ever,thedierencebecomescleartowardthetip;atz/Z=0.37,thedierenceisvisible.Atz/Z=0.69,thetimeaveragedmembranewingshiftthelowestpressurepointdownstreamcomparingwiththerigidwingresult.Thisshifthelpstodelayowseparation. PAGE 94 Figure6{26: Timeaveragedairfoilshapesatdierentspanwisepositionsforthemembranewingat=6oand15o. Figure6{27: Comparisonofcpcontoursat=6o:(left)timeaveragedmembranewing;(right)steadyrigidwing. PAGE 95 (a) (b) (c) Figure6{28: Chordwisecpcomparisonatdierentspanwiselocationsbetweentimeaveragedmembranewingandsteadyrigidwingat=6o:(a)z/Z=0;(b)z/Z=0.37;(c)z/Z=0.69. PAGE 96 InChapter 6 wehaveperformedtheaerodynamicstudyofthewingaerodynamicsofboththerigidandexiblewings.Inthischapterwewilldevelopasystematicapproachtoenhancethemembranewingperformance. 7.1 Introduction Itisdesirablenotonlytounderstandthepertinentphysicalphenomenaofthemembranewingdynamics,buttoimproveitsperformancebasedontheunderstand-ing.InChapter 6 wehaveinvestigatedbothrigidandmembranewingaerodynamics.Thewingshapeadoptedisfromanexistingdesign.Onequestionthenarisesnatu-rally:isthatthebestwingdesign,or,canweimproveitsperformance.Toanswerthisquestion,itisworthwhiletopauseandtakealookathowtheexistingwingisde-signed.Thedesignofthewingismostlybasedonobservationsandatwo-dimensionalanalysistooltotesttheperformanceofdierentairfoils.Basedonthetests,oneair-foilshapeisthenchosentoformathree-dimensionalwing.Thewingshapeisthentestedwithight.Inthistrialanderrorprocess,thethree-dimensionaleectsandthemembranematerialpropertiesarenotconsidered.Asystematicapproachisnecessarytoimprovethewingperformance. Therearemanyexistingmaturetechniquesforshapeoptimization.Couldwesimplyuseoneofthemtoachieveourgoal?Theanswerisnegative.AsdiscussedinChapter 6 ,theMAVwingperformsunderthelowReynoldsnumbercondition,anditexhibitsahighfrequencyvibrationevenundersteadyfreestreamconditions.Therefore,theoptimizationofamembranewingintroducesmultiplechallenges.First,coupled,timedependentsimulationsoftheinteractionbetweenviscousuidowandexiblestructureareveryexpensive.Second,anecientandautomaticgrid 88 PAGE 97 regenerationisessential.Todealwiththesediculties,weproposethefollowingstrategies: 6 32 51 63 ],inourstudythesur-rogatemodelistherigidwing.Weshalloptimizetherigidwingunderagivenowcondition,andevaluatetheperformanceofaexiblewingwithoutputfromthesurrogatemodel. Inthefollowing,werstintroducetheoptimizationproblem.Thisisfollowedbyanoptimizationalgorithminwhichweshowhowtointegratetheanalysistool,themovinggridtechnique,andtheoptimizerDesignOptimizationTools(DOT)[ 12 ]intoasystemtodoshapeoptimization.Theoptimizationprocessbeginswitharigidwingsurrogateoptimization,thentheexiblewingperformanceisassessedwiththeoptimizedshape.Optimizationresultsandtheaerodynamiccharacteristicsofbothoptimizedrigidandmembranewingsarethenhighlighted. 7.2 OptimizationandComputationalFramework Theobjectiveistoimprovethelift-to-dragratioofamembranewing.Thebaselinemembranewingshapeselectedisanexistingdesign[ 33 34 ]devisedfromsolutionsprovidedbyXFOIL[ 15 ]andempiricalmodicationbasedonighttests[ 34 ].XFOILprovidestwo-dimensionaluidowsolutionsbasedoncoupledthinlayerandinviscidowmodels.AschematicbaselinewingshapeisshowninFig. 6{1 inChapter 6 .Thebaselineshapehasthreecarbon-ber-madebattensthatarelabelledalsoinFig. 6{1 .Theshadingareacorrespondstothefuselage.Thefuselageismadeofrigidcarbonberprepregclothwhichdoesnotdeformduringtheight.Aexible PAGE 98 latexmembranecoversthetopsurfacewhichchangesitsshapeunderaerodynamicforces. Thebaselinewinghasavariablespanwisecamber.Thenominalcamberofthebaselinedesignis7.5%attheroot,and2%atthetip.Thewinghasamaximalchordlengthof13.7cm,ameanchordof10.5cm,andaspanof15.2cm.Cisusedtoindicatethechordlengthatdierentspanstations,andZforthehalfspan.Theangleofattackisdenedinreferencetotherootgeometry.Atthedesignpoint,theMAVieswithaspeedof10m/sandanangleofattackof6o.Theobjectiveistomaximizethelift-to-dragratioatsuchaightcondition.AtsuchaconditiontherootchordReynoldsnumberis9104. Sixdesignvariablesarechosen,theyarethey-coordinatesatsixpointsonthesurface.Threearelocatedattherootwhosepositionsareapproximatelyat10%,27%,and77%ofthechord,asseeninFig. 7{1 .TheotherthreearelocatedatBatten2withthesamerelativechordpositions. Thewingsurfacecanbeimplicitlyrepresentedbyafunctiony=f(x;z),wherex,y,andzrepresentthewingsurfacecoordinates.WeindicatethedesignvariablewithacapitalletterYiwhosevalueistheperturbationoftheycoordinateatpointiwithrespecttoitsbaselinevalue.Weinterpolatetheperturbationsfromthesixdesignvariablesoverthewingsurfacewiththethinplatespline(TPS)interpolationmethoddiscussedinChapter 6 .Therefore,theinterpolatedperturbationsoverthewingsurfaceis whereYisavectorindicatingthedesignvariablevalues,Gistheinterpolationoperator,andyistheinterpolatedvalueoverthewingsurface.TheperturbedwingsurfaceisthereforerepresentedbythesummationofthebaselineandtheinterpolatedvalueinEq.( 7.1 ),ynew=ybaseline+y.OneadvantageofthisperturbationapproachisthatwecanrecovertheoriginalshapewhentheperturbationvectorYiszero.In PAGE 99 (a) (b) Figure7{1: Acrosssectionofmembranewingalongabatten:(a)straightlinesindicatinghowconvexityconstraintisapplied;(b)eectofycoordinateperturbationatapoint.Thedesignpointsarelocatedat10%,27%,and77%ofthechord. PAGE 100 ouroptimizationprocess,boththeleadingandtrailingedgesarexedtomaintaintheidenticalangleofattackbetweenthebaselineshapeandtheoptimizedshape. Theoptimizationprobleminvestigatedcanbeformulatedasfollows. MaximizeCL=CDSubjectto1:CLCLbaseline2:Convexityconstraint:Y1+y1y0 Constraint1requiresthattheliftcoecientbenolessthanthatofthebaselinewing.Constraint2maintainstheconvexityoftheairfoil(seeFig. 7{1 a)sotoeliminateob-viouslyinappropriateshapes.Onlyoneofthesixconstraintsisshownforillustrationpurpose.Constraint3givesthelowerandupperboundsofthedesignvariableswhosevaluesarelistedinTable 7{1 Table7{1: Designvariablesboundsandtheiroptimalvaluesat6o. Parameter Initialdesign(cm) Lowerbound(cm) Upperbound(cm) Case1:opti-malvalueswith600iterationsperanalysis Case2:optimalvalueswith1000iterationsperanalysis -0.40 0.00 -0.258 -0.259 -0.40 0.00 -0.400 -0.400 -0.20 0.20 -0.152 -0.139 -0.20 0.20 -0.100 -0.095 -0.10 0.20 -0.082 -0.068 -0.10 0.20 0.135 0.121 Analysis 97 108 Cycle 6 6 0.529 0.529 0.529 7.55 7.55 Numberofanalyseswithconvexityviolations 13 22 PAGE 101 TheoptimizationprocedurebeginswiththebaselinedesignbycallingtheNavier-Stokesequationssolverrst.Theoptimizerthenperturbsoneormoredesignvari-ablesforapotentialdecliningdirection.Afterthat,theperturbationswillbeinter-polatedwiththeTPSmethodshowninEq.( 7.1 ).AnewcomputationalgridisthengeneratedbasedonthenewcongurationwiththemovinggridtechniquediscussedinChapter 3 .TheNavier-Stokesequationssolveriscalledtoevaluatetheobjec-tivefunctionbasedonthenewcomputationalgridagain.Adesigncycleincludestheevaluationofderivativesoftheobjectivefunctionandconstraints(herebynitedierences)andaonedimensionalsearchinadirectionselectedbasedonthederiva-tives.ThesearchterminatesafterafewcycleswhentheoptimizerDOTconvergencecriteriaaresatised.Foranobjectivefunctionwithndesignvariables,atleastnanalysesareneededtoevaluatethegradientoftheobjectivefunctionwithoutcon-straints.Thenumberofanalysisincreasesforconstraintoptimization.WiththeTPSinterpolation,theinterpolatedsurfacehasasecond-ordercontinuity.And,asshowninFig. 7{1 b,theinterpolatedsurfacepassesthroughallthecontrolpoints. Toevaluatetheaerodynamicsperformanceofthemembranewing,weusethecoupledmodeltoperformauidandstructureinteractioncomputation.Thiscou-pledmodelconsistsofthreemodules:aNavier-Stokesequationssolver,amembranestructuralsolver,andaninterfacingtechnique.Theuidsolverisapressure-basedmethodforincompressibleowthatwementionedinChapter 2 ;thestructuralsolverisaniteelementapproachforthedynamicresponseofamembranematerialwhichobeystheMooney-Rivlinhyperelasticmodel[ 64 ]. 7.3 ResultsandDiscussion 7.3.1 SensitivityAnalysis Whenemployingagradient-basedoptimizationmethod,akeyquestionisthesensitivityoftheaerodynamiccoecientsinresponsetotheminutegeometricchange.ToanswerthisquestionwerstrecallthegridsensitivitystudyinChapter 6 .Inthat PAGE 102 studyweconcludethatthefrictiondragissensitivetothegriddistribution.However,theliftandtheformdragonlyhavemodestvariationswiththegridrenement.Sincethefrictiondragisonlyasmallportionofthetotaldrag,thecoarsemeshisusedhereagainforshapeoptimization.Inthefollowing,weinvestigatehowtheuidsolver'sconvergencecriterionaectstheoptimizationresults. Dierentconvergencecriteriaareprobedtoassesstheirimpactsontheopti-mizationoutcome.TheimplicitSIMPLECmethodisusedfortheNavier-Stokesequations.Twoiterationnumbersarechosentointheoptimizationprocess:case1uses600iterationsperanalysis,whilecase2uses1000iterationsperanalysis.Inbothcasestheoverallresidual,denedasthesumoftheabsolutevaluesofthemassandmomentumuximbalanceofallcomputationalcell[ 88 ],isconsistentlyreducedbyanorderoffour,withCase2slightlymore.Bothcasesyieldacceptablesteadystatesolutions,withminordierencesintheresults.Intheoptimizationprocess,itisobservedthataftersixdesigncyclesbothcasesresearchthesamelocaloptimalalthoughtheirconvergencepathsaredierent.Eventhoughtheoptimalvaluesarethesame,thereissmalldierenceinthegeometryatthetip.Weplottheconver-gencehistoriesinFig. 7{2 .TheoptimizationresultsaretabulatedinTable 7{1 .InthefollowingweexclusivelyfocustheoutcomefromCase1asoursurrogatetostudythemembranewingperformance. Inourcomputations,weimposetheconvexityconstraintstoeliminateobviouslyunacceptableshapes.Inthesearchprocess,occasionally,violatingtheconstraintscouldacceleratetheconvergencetowardtheoptimum.However,frequentviolationswilldegradetheeciencyoftheoptimizationprocess.AsshowninTable 7{1 ,thereare13and22analyseswhichviolatethegeometryconstraintsforCase1andCase2,respectively. PAGE 103 Figure7{2: Optimizationdesignhistorywithdierentnumberofiterationsperanal-ysis. 7.3.2 PerformanceoftheOptimizedRigidWing Thebaselineshapehasthelargestcamberattheroot,about7.5%ofthechordlength;itdecreasesto2%atthetip.Theanalysisofthebaselinewingshowsthatthelargecamberneartherootcancauseowrecirculation.Theoptimizationreducesthecamberneartherootwhileincreasingthecambernearthetip,asseeninFig. 7{3 Theoptimizedshapehasamoreuniformcamber,graduallydecreasingfrom4.8%attherootto4%atthetip.Thismodicationresultsintheeliminationoftherecirculationintherootregion.Fig. 7{4 showsthestreamwisevelocityprolesofbothbaselineandoptimizedrigidwingsneartherootonthepressureside.Itdemonstratesthatthereverseowassociatedwiththebaselinewingshapehasnowbeeneliminatedintheoptimizedshape.Fig. 7{5 illustratesthespanwisepressuredistributionsofbothwings.Attheroot,thebaselinewingyieldsacross-overpressuredistributionattheleadingedge,whichdecreasesthetotallift.Theoptimizedshapeeliminatesthecross-over.Inaddition,theadversepressuregradientonthelower PAGE 104 (a) (b) (c) Figure7{3: Comparisonbetweenoptimizedshapeandbaselineshapeatdierentspanpositions.Thecamberdecreasesfrom4.8%attherootto4%atthetip:(a)z/Z=0;(b)z/Z=0.4;(c)z/Z=0.8. PAGE 105 Figure7{4: Comparisonofstreamwisevelocityprolesattherootbetweenbaselineandoptimizedshapes. surfaceissmootherfortheoptimizedshapethanforthebaseline;thishelpsexplainwhytheoptimizedwingdoesnotexperienceowseparationinthatregion. ThespanwiseL/DdistributionisdepictedinFig. 7{6 .Theaerodynamicim-provementislargelyrealizedinthe70%oftheinnerwing.Withareducedcamberintherootregion,thedragcoecientdecreasesnoticeablythere,whichleadstoahigheroveralllift-to-dragratio. Theoptimizationissetattheangleofattackof6o.Toprobetheperformanceoftheoptimizedshapeatotheranglesofattack,theaerodynamiccharacteristicsfordierentanglesofattackareshowninFig. 7{7 .Overall,theimprovementismoresubstantialatthemodestangles.ThespanwiseL/DdistributionsatdierentanglesareshowninFig. 7{8 .Consistently,theimprovementislargelyrealizedfromloweringtheformdrag. 7.3.3 FlexibleWingOptimization Withthesurrogateoptimizationoutput,themembranewingperformanceisevaluatedbasedonthecoupleduid-structureinteractionmodelinChapter 6 .To PAGE 106 (a) (b) (c) Figure7{5: Comparisonofpressurecoecientsonbaselineandoptimizedwingsatdierentspanwiselocationsat=6o:(a)z/Z=0;(b)z/Z=0.4;(c)z/Z=0.8. PAGE 107 Figure7{6: Comparisonofspanwiseaerodynamiccoecientsdistributionbetweenthebaselineandoptimizedwingsat=6o. PAGE 108 Figure7{7: Comparisonsofthebaselineandoptimizedrigidwingsatdierentanglesofattack. PAGE 109 (a) (b) Figure7{8: Comparisonofspanwiselift-to-dragratiobetweenthebaselineandopti-mizedrigidwingsatdierentanglesofattack.(a)=3o;(b)=9o. PAGE 110 investigatethespanwiseforcedistributionatdierenttimeinstants,theoverallaero-dynamicsofthemembranewings,forbothbaselineandoptimizedshapes,issum-marizedinFig. 7{9 .Consistentwiththerigidwingcase,showninFig. 7{7 ,theoptimizedshapeimprovestheL/Dconsistentlywithmoresubstantialimprovementatmodestanglesofattack.Furthermore,betweentherigidandmembranewings,themembranewingvarieslessintheL/Dversustheangleofattack,indicatingthatitcanoeramoresteadyightbehaviorthantherigidwing.Thisnalndingisconsistentwiththatobservedinighttests[ 34 ]. 7.4 SummaryandConclusions Wehaveperformedanoptimizationofaexiblemembranewingusingarigidwingassurrogate.Ananalysisofthenaldesignhasveriedthattheexiblewingcanbeimprovedbyoptimizationofarigidwingwiththesameun-deformedgeometry. Comparedtothebaseline,theoptimizationleadstolowercamberneartherootandhighercambernearthetip,whilestillleavingthecamberslightlyhigherattherootthanatthetip.Thelift-to-dragratioisimprovedoverarangeofanglesofattack.However,atlargeanglesofattack,theimprovementwiththeoptimizedshapediminished,theL/Dofthemembranewingvarieslesswithchangesintheangleofattack,indicatingthatitcanoeramoresteadyightbehaviorthattherigidwing.Theimprovementislargelylocatedwithin70%oftheinnerwing.Theliftoftheoptimizedwingremainsthesameasthebaselinewing,eventhoughitscamberreducesattheroot.Theimprovementinaerodynamicsoftheoptimizedwingislargelyrealizedvisreducingtheformdrag.Forbothrigidandmembranewings,theoptimizationimprovestheL/Dconsistently. PAGE 111 Figure7{9: Theoptimizedmembranewingperformanceatdierentanglesofattack. PAGE 112 WehaveinvestigatedMAVwingaerodynamicsandshapeoptimizationwithapplicationstoMAVs.Eveninthesteadyfreestream,themembranewingexhibitsself-initiatedvibration.Toaccuratelysimulatetheinteractionbetweentheexiblemembranestructureanditssurroundingviscousow,wehaveperformedcoupleduidandstructurecomputations.Intermsofcomputationalrequirementsforthecoupleduidandstructureinteractions,wehaveenforcedthetimesynchronizationbetweentheuidandstructuralsolverstoreducethephaselagerror.Wehavealsoconrmedtheimportanceofthegeometricconservationinthecontextofmovinggridproblems. InregardtoMAVaerodynamics,tipvorticesplayanimportantroleforthelowaspectratiowing.Ontheonehand,tipvorticeslowertheliftbyreducingtheeectiveangle;ontheotherhand,tipvorticesprovideadditionalliftbyforminglowpressureregionsontheupperwingsurface.FlowontheupperwingsurfaceispronetoseparateunderlowReynoldsnumbercondition.Eventhough,MAVsmaintainaveryhighhighstallangleofattack.Tipvorticesandowseparationcauseunsteadybehaviors.Vortexsheddingoccursathighanglesofattack;however,thevortexsheddingdoesnothavemuchimpactontheaerodynamiccoecients. WehaveconductedCFD-basedoptimizationinamethodicalmanner.wehaveoptimizedamembranewingusingtherigidwingasthesurrogatemodel.Comparingwiththebaseline,theoptimizedwinghasalowercamberneartherootandahighercambernearthetip.However,theoptimizedwingstillleavesthecamberslightlyhigherattherootthanatthetip.Eventhoughtheoptimizedwinghaslowercamber 104 PAGE 113 thanthebaselineattheroot,itmaintainsthesameliftasthebaseline.Theim-provementinthelift-to-dragratioislargelyrealizedviareducingtheformdragandabetterpressuredistribution.Ouroptimizationhasorientedatanangleofattackof6o;testsoverarangeofanglesofattackshowdierentlevelsofimprovement.Atlargeanglesofattack,theimprovementwiththeoptimizedshapediminishes.Atlowanglesofattack,theimprovementismostlyachievedwithin70%oftheinnerwing. Muchprogresshasbeenmadeinthisexcitingresearcharea.Withthecontinuousimprovementinsimulation,materials,fabrication,andmeasurementtechniques,aswellasthefastdevelopmentinmicro-systems,signicantpotentialexiststofurtheradvancethemicroairvehicleasavehicleplatformforperformingnumerousmissions,andasatechnologyplatformforevaluatingsystem.ThefutureworkinthisareawillbeacombinednumericalandexperimentalstudyonthelowReynoldsnumberandlowaspectratioMAVaerodynamics.Theprincipalfeaturesoffuturestudyare ThesearethekeyelementsforanintegratedstudyoftheMAVaerodynamics.Thefulllmentofthisworkhasbothpracticalandtheoreticalsignicance.Practically,thisworkwillimproveourunderstandingontheMAVaerodynamics,enhancetheMAVperformance;theoretically,itwillshedlightonthepassivecontroltheory,andprovidedeepinsightforthetransitionmechanism. Thewell-developedandvalidatedcomputationalalgorithmpresentedherecanbeutilizedforthisresearch.ThecurrentnumericalworksinMAVsexclusivelyconsiderthemembranewingitself.Thefuselageandthepropellerarenotconsidered.Inthefuturework,wewillextendourformerworktoconsiderthefuselageandpropellereects. PAGE 114 ExperimentsbyWasazketal.[ 111 ]testedanoldMAVmodel,whichisquitedierentfromthecurrentonewestudyhere.ThecurrentMAVmodelfeaturesastreamlizedbodyandanbuilt-inelectricmotor.Thenewmodelissuperiortotheoldoneintheaerodynamicperformance,endurance,control,andmaneuverability.Therefore,acombinedexperimentalandnumericalstudyonthecurrentMAVmodelisneeded. TheMAViesattherelativelylowReynoldsnumberconditions.Complicatedowstructuresoccurattheupperwingsurface,includinglaminarboundarylayerseparation,laminartoturbulentowtransition,andturbulentowreattachment.Sincethereisnocomprehensivetheorytoaccountforallthese,itisimportanttoinvestigatethesephenomenonwithbothexperimentsandnumericalsimulations.Fornumericalcomputations,duetothelowReynoldsnumbercondition,computationswithrenedresolutionsareneededtoaccuratelypredictlaminar-turbulenttransition,detailedowstructures,andtheresultingaerodynamiccoecients. Astothedevelopmentofoptimizationtechniques,oneshouldexplorenewmethodtooptimizethewingshape.Insteadofconsideringasinglepointsolution,weneedtotakeotheranglesofattackinourconsiderationandconsideramulti-objectiveoptimizationproblem. PAGE 115 Thecomponentsofthenonlinearpartofthestrain-displacementmatrix,BL,canbewrittenasthefollowing, 107 PAGE 117 [1] Alexander,H.,\TheTensileInstabilityofInitiallySphericalBalloons,"Interna-tionalJournalofEngineeringScience,Vol.9,1971,pp.151-162. [2] Anderson,J.D.,Jr.,IntroductiontoFlight,3rded,McGraw-Hill,1989. [3] Batina,J.,\UnsteadyEulerAirfoilSolutionsUsingUnstructuredDynamicMeshes,"AIAAPaper1989-0115. [4] Beaty,M.F.,\TopicsinFiniteElasticity:HyperelasticityofRubber,Elastomers,andBiologicalTissues-withExamples,"AppliedMechanicReview,Vol.40,1987,pp.1699-1734. [5] Bendiksen,O.O.,\ANewApproachtoComputationalAeroelasticity,"1991,AIAAPaper91-0939-CP. [6] Burgee,S.,Giunta,A.A.,Narducci,R.,Watson,L.T.,Grossman,B.,andHaftka,R.T.,\ACoarseGrainedParallelVariable-ComplexityMultidisciplinaryOptimizationParadigm,"TheInternationalJournalofSupercomputerApplica-tionsandHighPerformanceComputing,Vol.10,1996,pp.269-299. [7] Carmichael,B.H.,\LowReynoldsNumberAirfoilSurvey,"NASACR-165803,1981. [8] Chasman,D.,andChakravarthy,S.,\ComputationalandExperimentalStudiesofAsymmetricPitch/PlungeFlapping-TheSecretofBiologicalFlyers,"AIAAPaper2001-0859. [9] Crompton,M.J.,andBarrett,R.V.,\InvestigationoftheSeparationBubbleFormedBehindtheSharpLeadingEdgeofaFlatPlateatIncidence,"JournalofAerospaceEngineering,Vol.214,2000,pp.157-176. [10] Cummings,R.M.,Morton,S.A.,Siegel,S.G.,andBosscher,S.,\NumericalPredictionandWindTunnelExperimentforaPitchingUnmannedCombatAirVehicle,"AIAAPaper2003-0417. [11] DeLaurier,J.D.,\AnAerodynamicModelforFlappingWingFlight,"Aeronau-ticalJournal,1993,pp.125-130. [12] [13] Dickinson,M.H.,Lehmann,F.,andSane,S.P.,\WingRotationandtheAero-dynamicBasisofInsectFlight,"Science,Vol.284,1999,pp.1954-1960. 109 PAGE 118 [14] Ding,H.,Yang,B.,Lou,M.,andFang,H.,\NewNumericalMethodforTwo-DimensionalPartiallyWrinkledMembranes,"AIAAJournal,Vol.41,2003,pp.125-132. [15] Drela,M.,andGiles,M.B.,\Viscous-InviscidAnalysisofTransonicandLowReynoldsNumberAirfoils,"AIAAJournal,Vol.25,1987,pp.1347-1355. [16] Duchon,J.P.,\SplinesMinimizingRotation-InvariantSemi-NormsinSobolevSpaces,"ConstructiveTheoryofFunctionsofSeveralVariables,Oberwolfach1976,editedbySchempp,W.andZeller,K.,Springer-Verlag,Berlin,1977,pp.85-100. [17] Ellington,C.P.,\TheAerodynamicsofHoveringInsectFlight,I.TheQuasi-SteadyAnalysis,"Philos.Trans.R.Soc.LondonSer.A,1984,Vol.35,pp.1-15. [18] Eriksson,L.E.,\GenerationofBoundary-ConrmingGridsAroundWing-BodyCongurationsUsingTransniteInterpolation,"AIAAJournal,Vol.20,1982,pp.1313-1320. [19] Fuentes,C.,He,X.,Carroll,B.,Lian,Y.,andShyy,W.\LowReynoldsNumberFlowsAroundandAirfoilwithaMovableFlap-Part1:Experiments,"AIAAPaper2000-2239. [20] Gad-el-Hak,M.,\Micro-Air-Vehicles:CanTheybeControlledBetter,",tex-titJournalofAircraft,Vol.38,2001,pp.419-429. [21] Gault,D.E.,\AnInvestigationatLowSpeedoftheFlowOveraSimulatedFlatPlateatSmallAnglesofAttackUsingPitotStaticandHot-wireProbes,"1957,NACATN-3876. [22] Gordnier,R.E.,andFithen,R.,\CouplingofaNonlinearFiniteElementStruc-turalMethodwithaNavier-StokesSolver,"ComputersandStructures,Vol.81,2003,pp.75-89. [23] Gordnier,R.E.,andVisbal,M.R.,\DevelopmentofaThree-DimensionalVis-cousAeroelasticSolverforNonlinearPanelFlutter,"JournalofFluidsandStructures,Vol.16,2002,pp.497-527. [24] Grasmeyer,J.M.,andKeennon,M.T.,\DevelopmentoftheBlackWidowMicroAirVehicle,"AIAAPaper2001-0127. [25] Green,A.E.,andAdkins,J.E.,LargeElasticDeformations,TheClarendonPress,Oxford,1960. [26] Farhat,C.,Geuzaine,P.,andBrown,G.,\ApplicationofaThree-eldNonlinearFluid-structureFormulationtothePredictionoftheAeroelasticParametersofanF-16Fighter,"ComputersandFluids,Vol.32,pp.3-29. PAGE 119 [27] Farhat,C.,Geuzaine,P.,andGrandmont,C.,\TheDiscreteGeometricConser-vationLawandtheNonlinearStabilityofALESchemesfortheSolutionofFlowProblemsonMovingGrids,"JournalofComputationalPhysics,Vol.174,2001,pp.669-694. [28] Farhat,C.,andLesoinne,M.,\TwoEcientStaggeredAlgorithmsfortheSe-rialandParallelSolutionofThree-DimensionalNonlinearTransientAeroelasticProblems,"ComputerMethodsinAppliedMechanicsandEngineering,Vol.182,2000,pp.499-515. [29] Hartwich,P.M.,andAgrawal,S.,\MethodforPerturbingMultiblockPatchedGridsinAeroelasticandDesignOptimizationApplications,"AIAAPaper1997-2038. [30] He,X.,Fuentes,C.,Shyy,W.,Lian,Y.,andCarroll,B.,\ComputationofTransitionalFlowsaroundanAirfoilwithaMovableFlap,"AIAAPaper2000-2240. [31] Herbert,T.,\ParabolizedStabilityEquations,"Annu.Rev.Fluid.Mech.,Vol.29,1997,pp.245-283. [32] Hutchison,M.G.,Unger,E.R.,Mason,W.H.,Grossman,B.,andHaftka,R.T.,\Variable-ComplexityAerodynamicOptimizationofaHigh-SpeedCivilTransportWing,"JournalofAircraft,Vol.31,1994,pp.110-116. [33] Ifju,P.G.,Ettinger,S.,Jenkins,D.,andMartinez,L.,\CompositeMaterialsforMicroAirVehicles,"SAMPEJournal,Vol.37,2001,pp.7-12. [34] Ifju,P.,Jenkins,D.,Ettinger,S.,Lian,Y.,Shyy,W.,andWaszak,R.M.,\Flexible-Wing-BasedMicroAirVehicles,"AIAAPaper2002-0705. [35] Issa,R.I.,\SolutionoftheImplicitlyDiscretisedFluidFlowEquationsbyOperator-Splitting,"JournalofComputationalPhysics,Vol.62,1985,pp.40-65. [36] Jackson,P.S.,\ASimpleModelforElasticTwo-DimensionalElasticSails,"AIAAJournal,Vol.21,1983,pp.153-155. [37] Jackson,P.S.,\TheoryforConicalMembraneWingsofHighAspectRatio,"AIAAJournal,Vol.39,2001,pp.781-786. [38] Jackson,P.S.,andChristie,G.W.,\NumericalAnalysisofThree-DimensionalElasticMembraneWings,"AIAAJournal,Vol.25,1987,pp.676-682. [39] Jameson,A.,\TimeDependentCalculationsUsingMultigrid,withApplicationstoUnsteadyFlowspastAirfoilsandWings,"AIAAPaper1991-1596. [40] Jenkins,C.H.,\NonlinearDynamicResponseofMembranes:StateoftheArt{Update,"AppliedMechanicReview,Vol.49,1996,pp.S41-S48. PAGE 120 [41] Jenkins,C.H.,andLeonard,J.W.,\NonlinearDynamicResponseofMem-branes:StateoftheArt,"AppliedMechanicReview,Vol.44,1991,pp.319-328. [42] Jones,B.M.,\Stalling,"J.R.Aero.Soc.Vol.38,1938,pp.747-770. [43] Jones,K.D.,andPlatzer,M.F.,\AnExperimentalandNumericalInvestigationofFlapping-WingPropulsion,"AIAAPaper1999-0995. [44] Jones,K.D.,andPlatzer,M.F.,\Flapping-WingPropulsionforaMicroAirVehicle,"AIAAPaper2000-0897. [45] Jones,K.D.,andPlatzer,M.F.,\ExperimentalInvestigationoftheAerody-namicCharacteristicsofFlapping-WingMicroAirVehicles,"AIAAPaper2003-0418. [46] Jones,W.P.,andLaunder,B.E.,\ThePredictionofLaminarizationwithaTwo-EquationModelofTurbulence,"Int.J.HeatMassTrans.,Vol.15,1972,pp.301-314. [47] Kalro,V.,andTezduyar,T.E.,\AParallel3DComputationalMethodforFluid-StructureInteractionsinParachuteSystems,"ComputerMethodsinAppliedMechanicsandEngineering,Vol.190,2000,pp.321-332. [48] Kamakoti,R.,Berg,M.,Ljungqvist,D.,andShyy,W.,\AComputationalStudyforBiologicalFlappingWingFlight,"TransactionofAeronauticalandAstronau-ticalSocietyoftheRepublicofChina,KeynotePaper,Vol.32,2000,pp.265-279. [49] Kamakoti,R.,Lian,Y.,Regisford,S.,Kurdila,A.,andShyy,W.,\Compu-tationalAeroelasticityUsingaPressure-basedSolver,"ComputerModelinginEngineering&Sciences,Vol.3,2002,pp.773-790. [50] Katz,J.,Low-SpeedAerodynamics:FromWingTheorytoPanelMethods,Mc-Graw-Hill,SanFrancisco,CA,1979. [51] Knill,D.L.,Giunta,A.A.,Baker,C.A.,Grossman,B.,Mason,W.H.,Haftka,R.T.,andWatson,L.T.,\ResponseSurfaceMethodsCombiningLinearandEulerAerodynamicsforSupersonicTransportDesign,"JournalofAircraft,Vol.36,1999,pp.75-86. [52] Launder,B.E.,andSpalding,D.B.,\TheNumericalComputationofTurbulentFlows,"Comp.Meth.Appl.MechEngg.,3,1974,pp.269-289. [53] Lian,Y.,andShyy,W.,\Three-DimensionalFluid-StructureInteractionsofaMembraneWingforMicroAirVehicleApplications,"AIAAPaper2003-1726. [54] Lian,Y.,Shyy,W.,andHaftka,R.,\ShapeOptimizationofaMembraneWingforMicroAirVehicles,"AIAAPaper2003-0106.AcceptedforpublicationatAIAAJournal. PAGE 121 [55] Lian,Y.,Shyy,W.,Ifju,P.,andVerron,E.,\AComputationalModelforCou-pledMembrane-FluidDynamics,"AIAAPaper2002-2972.Acceptedforpubli-cationatAIAAJournal. [56] Lian,Y.,Steen,J.,Trygg-Wilander,M.,andShyy,W.,\LowReynoldsNumberTurbulentFlowsaroundaDynamicallyShapedAirfoil,"ComputersandFluids,Vol.32,2003,pp.287-303. [57] Lighthill,M.J.,\HydrodynamicsofAquaticAnimalPropulsion,"Ann.Rev.FluidMech.,1969,pp.413-445. [58] Liu,F.,Zhu,Y.,Tsai,H.M.,andWong,A.S.F.,\CalculationofWingFlutterbyaCoupledFluid-StructureMethod,"JournalofAircraft,Vol.38,2001,pp.334-342. [59] Liu,H.,andKawachi,K.,\ANumericalStudyofInsectFlight,"Vol.J.Comput.Phys.,146,1998,pp.124-156. [60] Liu,H.T.,\UnsteadyAerodynamicsofaWortmannWingatLowReynoldsNumber,"JournalofAircraft,1992,Vol.29,pp.532-539. [61] Lissaman,P.B.S.,\Low-Reynolds-numberAirfoils,"Annu.Rev.FluidMech.Vol.15,1983,pp.223-239. [62] Lorillu,O.,andHureau,J.,\NumericalandExperimentalAnalysisofTwo-DimensionalSeparatedFlowsoveraFlexibleSail,"J.Fluid.Mech.,Vol.466,pp.319-341. [63] Mason,B.H.,Haftka,R.T.,Johnson,E.R.,andFarley,G.L.,\VariableCom-plexityDesignofCompositeFuselageFramesbyResponseSurfaceTechniques,"ThinWallStructures,Vol.32,1998,pp.235-261. [64] Mooney,M.,\ATheoryofLargeElasticDeformation,"JournalofAppliedPhysics,Vol.11,1940,pp.582-592. [65] Mueller,T.J.(ed.),ProceedingsoftheConferenceonLowReynoldsNumberAerodynamics,Univ.ofNotreDame,NotreDame,IN,1989. [66] Mueller,T.J.(ed.),ProceedingsoftheConferenceonFixed,FlappingandRotaryWingVehiclesatVeryLowReynoldsNumbers,Univ.ofNotreDame,NotreDame,IN,2000. [67] Mueller,T.J.,andDeLaurierJ.D.,\AerodynamicsofSmallVehicles,"AnnualReviewofFluidMechanics,Vol.35,2003,pp.89-111. [68] Nelsen,J.N.,\TheoryofFlexibleAerodynamicsSurfaces,"JournalofAppliedMechanics,Vol.30,1963,pp.435-442. PAGE 122 [69] Newman,B.G.,\AerodynamicsTheoryforMembraneandSails,"ProgressinAerospaceSciences,Vol.24,1987,pp.1-27. [70] Oden,J.T.,andSato,T.,\FiniteStrainsandDisplacementsofElasticMem-branesbytheFiniteElementMethod,"InternationalJournalofSolidsandStructures,3,1967,pp.471-488. [71] Oliveira,P.J.,andIssa,R.I.,\AnImplicitPISOAlgorithmfortheComputationofBuoyancy-DrivenFlows,"NumericalHeatTransfer,PartB,Vol.40,2001,pp.473-493. [72] Patankar,S.V.,andSpalding,D.B.,\ACalculationProcedureforHeat,MassandMomentumTransferinThree-DimensionalParabolicFlows,"Int.J.HeatMassTransf.,Vol.15,1972,pp.1787-1806. [73] Patel,V.C.,Rodi,W.,andScheuerer,G.,\TurbulenceModelsforNear-WallandLowReynoldsNumberFlows,"AIAAJournal,23,1985,pp.1308-1319. [74] Pulliam,T.,\TimeAccuracyandtheUseofImplicitMethods,"1993,AIAAPaper93-3360-CP. [75] Reed,H.L.,andSaric,W.S.,\LinearStabilityTheoryAppliedtoBoundaryLayers,"Annu.Rev.Fluid.Mech.,Vol.28,1996,pp.389-428. [76] Rempfer,D.,\Low-DimensionalModelingandNumericalSimulationofTransi-tioninSimpleShearFlows,"AnnualReviewofFluidMechanics,Vol.35,2003,pp.229-265. [77] Reuther,J.,Alonso,J.J.,Rimlinger,M.J.,andJameson,A.,\AerodynamicShapeOptimizationofSupersonicAircraftCongurationviaanAdjointFormu-lationonDistributedMemoryParallelComputers,"Comput.Fluids,28,1999,pp.675-700. [78] Robinson,B.A.,Batina,J.T.,andYang,H.T.Y.,\AeroelasticAnalysisofWingsUsingtheEulerEquationswithaDeformingMesh,"JournalofAircraft,Vol.28,1991,pp.781-788. [79] Rumsey,C.L.,Sanetrik,M.D.,Biedron,R.T.,Melson,N.D.,andParlette,E.B.,\EciencyandAccuracyofTime-AccurateTurbulentNavier-StokesCom-putations,"ComputersandFluids,Vol.25,1996,pp.217-236. [80] Schuster,D.M.,Liu,D.D.,andHuttsell,L.J.,\ComputationalAeroelasticity:Success,Progress,Challenge,"AIAAPaper2003-1725. [81] Schuster,D.,Vadyak,J.,andAtta,E.,\StaticAeroelasticAnalysisofFighterAircraftUsingaThree-DimensionalNavier-StokesAlgorithm,"J.Aircraft,Vol.27,1990,pp.820-825. PAGE 123 [82] Shyy,W.,ComputationalModelingforFluidFlowandInterfacialTransport,Elsevier,Amsterdam,TheNetherlands,1994. [83] Shyy,W.,Berg,M.,andLjungqvist,D.,\FlappingandFlexibleWingsforBi-ologicalandMicroVehicles,"ProcessinAerospaceSciences,Vol.35,1999,pp.455-506. [84] Shyy,W.,Jenkins,D.A.,andSmith,R.W.,\StudyofAdaptiveShapeAirfoilsatLowReynoldsNumberinOscillatoryFlow,"AIAAJournal,Vol.35,1997,pp.1545-1548. [85] Shyy,W.,andMittal,R.,\SolutionMethodsfortheIncompressibleNavier-StokesEquations,"inR.Johnson(ed.)HandbookofFluidDynamics,CRCPress,BocaRaton,FL,1998,pp.31-1to31-33. [86] Shyy,W.,andSmith,R.W.,\ComputationofLaminarFlowandFlexibleStruc-tureInteraction,"inM.HafezandK.Oshima(eds.)ComputationalFluidDy-namicsReview,JohnWiley&Sons,1995,pp.777-796. [87] Shyy,W.,Thakur,S.S.,Ouyang,H.,Liu,J.,andBlosch,E.,ComputationalTechniquesforComplexTransportPhenomena,CambridgeUniversityPress,1997. [88] Shyy,W.,Udaykumar,H.S.,Rao,M.M.,andSmith,R.W.,ComputationalFluidDynamicswithMovingBoundaries,1996,Taylor&Francis,Washington,DC. [89] Smith,M.J.,\SimulatingMothWingAerodynamics:TowardstheDevelopmentofFlapping-WingTechnology,"AIAAJournal,Vol.34,1996,pp.1348-1355. [90] Smith,M.J.\ComputationalConsiderationsofanEuler/Navier-StokesAeroe-lasticMethodforaHoveringRotor,"JournalofAircraft,Vol.33,1996,pp.429-434. [91] Smith,M.J.,Hodges,D.H.,andCesnik,C.E.S.,\EvaluationofComputationalAlgorithmsSuitableforFluid-StructureInteractions,"JournalofAircraft,Vol.37,2000,pp.282-294. [92] Smith,R.W.,andShyy,W.,\ComputationalModelofFlexibleMembraneWingsinSteadyLaminarFlow,"AIAAJournal,Vol.33,1995,pp.1769-1777. [93] Smith,R.W.,andShyy,W.,\IncrementalPotentialFlowBasedMembraneWingElement,"AIAAJournal,Vol.35,1997,pp.782-788. [94] Subbaraj,K.,andDokainish,M.A.,\ASurveyofDirectTime-IntegrationMeth-odsinComputationalStructuralDynamics-II.ImplicitMethods,"ComputersandStructures,Vol.32,1989,pp.1387-1401. PAGE 124 [95] Tani,I.,\Low-SpeedFlowsInvolvingBubbleSeparations,"ProgressinAeronau-ticalScience,Vol.5,editedbyD.Kuchemannandl.H.G.Sterne,Pergamon,NewYork,1964,pp.70-103. [96] Templin,R.J.,\TheSpectrumofAnimalFlight:InsectstoPterosaurs,"ProgressinAerospaceSciences,Vol.36,2000,pp.393-436. [97] Torres,G.E.,andMueller,T.J.,\AerodynamicCharacteristicsofLowAspectRatioWingsatLowReynoldsNumber."InFixedandFlappingWingAerody-namicsforMicroAirVehicleApplications,ed.TJMueller,Vol.195,2001,pp.191-213.Reston,VA:AIAA. [98] Thakur,S.,Wright,J.,andShyy,W.,STREAM:AComputationalFluidDynam-icsandHeatTransferNavier-StokesSolver.TheoryandApplications,StreamlineNumerics,Inc.,andComputationalThermo-FluidsLaboratory,DepartmentofMechanicalandAerospaceEngineeringTechnicalReport,Gainesville,Florida,2002. [99] Thaokar,R.M.,andKumaran,V.,\StabilityofFluidFlowPastaMembrane",J.FluidMech.,Vol.472,2002,pp.29-50. [100] Thomas,P.D.,andLombard,C.K.,\GeometricConservationLawandItsApplicationtoFlowComputationsonMovingGrids,"AIAAJournal,Vol.17,1979,pp.1030-1037. [101] Thwaites,B.,\AerodynamicsTheoryofSails-PartI,"ProceedingsoftheRoyalSociety,A261,1961,pp.402-442. [102] VanDoormaal,J.P.,andRaithby,G.D.,\EnhancementsoftheSIMPLEMethodforPredictingIncompressibleFluidFlows,"NumericalHeatTransfer,Vol.67,1985,pp.147-163. [103] Versteeg,H.K.,andMalalasekera,W.,AnIntroductiontoComputationalFluidDynamics,LongmanScientic&Technical,1995. [104] Verron,E.,Khayat,R.E.,Derdouri,A.,andPeseux,B.,\DynamicInationofHyperelasticSphericalMembranes,"JournalofRheology,Vol.43,1999,pp.1083-1097. [105] Verron,E.,Marckmann,G.,andPeseux,B.,\DynamicInationofNon-linearElasticandViscoelasticRubber-likeMembranes,"Int.J.Numer.Mech.Engng.,Vol.50,2001,pp.1233-1251. [106] Vest,M.S.,andKata,J.,\UnsteadyAerodynamicModelofFlappingWings,"AIAAJournal,Vol.34,1996,pp.1435-1440. [107] Vest,M.S.,andKatz,J.,\AerodynamicStudyofaFlappingWingMicro-UAV,"AIAA1999-0994. PAGE 125 [108] Voelz,K.,\ProlundLuftriebeinesSegels,"ZAMM,Vol.30,1950,pp.301-317. [109] Walker,J.A.\FunctionalMorphologyandVirtualModels:PhysicalCon-straintsontheDesignofOscillatingWing,Fins,Legs,andFeetatIntermediateReynoldsNumbers,"Integ.andComp.Biol.,Vol.42,2002,pp.232-242. [110] Wang,Z.J.,\VortexSheddingandFrequencySelectioninFlappingFlight,"JournalofFluidMechanics,Vol.410,2000,pp.323-341. [111] Waszak,R.M.,Jenkins,N.L.,andIfju,P.,\StabilityandControlPropertiesofanAeroelasticFixedWingMicroAerialVehicle,"AIAAPaper2001-4005. [112] Weis-Fogh,T.,\QuickEstimatesofFlightFitnessinHoveringAnimals,Includ-ingNovelMechanismsforLiftProduction,"JournalofExperimentalBiology,Vol.59,1973,pp.169-230. [113] Wilson,E.L.,Farhoomand,I.,andBathe,K,J.,\NonlinearDynamicAnalysisofComplexStructures,"EarthquakeEngineeringandStructuralDynamics,Vol.1,1973,pp.241-252. [114] Young,A.D.,andHorton,H.P.,\SomeResultsofInvestigationofSeparationBubbles,"AGARDCP,4,1966,pp.779-811. [115] Zwaan,R.J.,andPrananta,B.B.,\Fluid/structureInteractioninNumericalAeroelasticSimulation,"InternationalJournalofNon-LinearMechanics,Vol.37,2002,pp.987-1002. PAGE 126 YongshengLianwasborninthesmallvillageofZiBo,ShandongProvince,inthePeople'sRepublicofChina.Heearnedhisbachelor'sdegreefromShandongUniversity,majoringinmathematics.BeforehejoinedtheUniversityofFlorida,heearnedtwomaster'sdegrees(bothinmathematics):onefromtheChineseAcademyofSciences;andtheotherfromtheHongKongUniversityofScienceandTechnology.In1999hejoinedtheUniversityofFloridaasagraduatestudent.HewastherecipientoftheAlumniFellowshipattheUniversityofFlorida.Yongshengenjoyssportsevents.HeisafanofGatorsports. 118 |