UFDC Home  myUFDC Home  Help 



Full Text  
MEMBRANE AND ADAPTIVELYSHAPED 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 kE Turbulence Model ...... 13 2.2.2 Boundary Conditions for Turbulence Models . ... 15 2.3 Laminarturbulent 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 ADAPTIVELYSHAPED WING 27 4.1 Introduction ............. . . ...... 27 4.2 Experimental and Computational Setup . . ...... 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 MooneyRivlin Model ................ .. .. 38 5.2.2 Principle of Virtual Work .................. .. 39 5.2.3 Triangular Membrane Element . . ..... 40 5.2.4 GreenLagrange 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 STRAINDISPLACEMENT 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 ADAPTIVELYSHAPED 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 lifttodrag 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 adaptivelyshaped wings, utilizing piezoactuated flaps, have been studied. In the adaptivelyshaped wing study, we use piezoactuated 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 lifttodrag 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 threedimensional NavierStokes 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 iStokes solver, an automatic grid generation tool, and a gradientbased 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 lifttodrag 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 lifttodrag 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 WeisFogh [112], and Lighthill [57]. Recent experimental and numerical studies in this aspect were given by ('!ii ii,, 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 11 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 buildin 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 11: 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. 12, 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 11. 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   6Batten (Monofilm)   2Batten (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 12: 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]. GadelHak [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 selfinitiated 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('! iI1 ,, [38] adopted a three dimensional potential flowbased 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 twodimensional flexible membrane wing and surrounding viscous flows. A combined numerical and experimental ,a i1, ~i; of a flow over a twodimensional 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 adaptivelyshaped wing. First, we will develop a computational capability to study the fluid and structure interaction. Second, based on this capability, we will study the piezoactuated 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 NavierStokes equations for incompressible flow on curvilinear coordi nates are employ, l along with a twoequation ke 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 laminarturbulent transition model, a brief review is offered. In ('!i lpter 3 we devote to moving grid technology, which is used by both the coupled fluidstructure study and the optimization. A robust, efficient and effective moving grid technique is proposed, which uses transfinite interpolationbased 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 adaptivelyshaped flaps for a lowReynolds number wing. Results of the interaction between piezoactuated 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 gradientbased 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 adaptivelyshaped wings, second, our interests in micro air vehicles. Original work include a moving grid algorithm for multiblock 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 piezoactuated 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 PressureImplicit 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 threedimensional NavierStokes equations for incompressible flows. Both methods belong to the pressurebased approach [82] which devises an artificially derived pressure (correction) equation by manipulating the mass continuity and momentum equations. Threedimensional N ,i., rStokes 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 ke 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 threedimensional N ,1i. rStokes 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 + l12 + 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 zdirections, re spectively, p is the viscosity, G1, G2, and G3 are the body force components per unit volume in the x, y, and zdirections. 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 xdirection 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 underrelaxations 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 noniterative approach to handle the pressurevelocity 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 NavierStokes 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 bodyfitted 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 nonphysical 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 fluidstructure interaction prob lems were discussed by Shyy et al. [88]. 2.2 Turbulence Models 2.2.1 Reynolds Equations and kE 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 offdiagonal 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 kE 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 kE 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 kE turbulence model that we mentioned above is based on the high Reynoldsnumber 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 NavierStokes 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 lowReynoldsnumber versions of the kE [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 1v. rs: the outer 1v r and the surface 1i.r. In the outer l?vr the Reynolds stress is dominated; while in the surface l?vr both the Reynolds stress and the shear stress are important. In the outer l1.vr, 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 velocitydefect law expressed as max U g( ), (2.27) U, 6 where 6 is the boundary lir thickness, UmaxU is the velocity deficit, g is an unknown function with the order of 1. The velocity in the surface lvr is given by the law of the wall UT = f(y ). (2.28) Here y+ is a nondimensional 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 lowReynoldsnumber kE 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 lowReynoldsnumber 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 nearwall and low Reynolds number flows. A more detailed review about turbulence models for nearwall and low Reynolds number flows is given by Patel et al. [73]. 2.3 Laminarturbulent 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 lowdimensional models [76], to direct numerical simulation, like largeeddy 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 nonovi. 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 onedimensional perturbation method [56, 77] is can efficiently regenerate the grid. The onedimensional 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 onedimensional 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 3stage perturbation method, analogous to the transfinite interpolation (TFI) method [18], was proposed for threedimensional problems treated by the singleblock approach. Unlike the original TFI method, this perturbation method considers the original grid distribution, and preserves the original grid quality. Take a twodimensional case as an example, in Fig. 31, a rectangle undergoes shape change. To regenerate the grid for the changed geometry, the perturbation method proceeds as a 2stage 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 xdirection (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), 106) 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. 3lc]. For threedimensional problems, a 3stage 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 multiblock 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 crossover. 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 multiblock structured grid, CFD software often requires pointmatched grid block interfaces. A basic requirement for grid regeneration in moving boundary problems is to maintain pointmatched 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 offbody 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 offbody points. While the task of redistributing interior grids can be (d ill. iii] in practice, it is more difficult to move the vertices of each block if a pointtopoint 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. 32a, 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. 32a, 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 offbody 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 offbody vertex, and, b, a body surface point. For each offbody 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 32: 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 threedimensional parameter space (i, j, k) into a onedimensional 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 threedimensional, multiblock, and structured computational grid. The initial configuration is shown in Fig. 32a, block 1 shares the interface with blocks 2 and 3. The initial grid is shown in Fig. 32b, it has clustered distributions at the top and bottom surfaces, and it also has pointtopoint matched interfaces. The bottom two faces are the moving boundaries that experience large deformation. The coordinates of each vertex are labelled. Comparing Fig. 32a and Fig. 32c 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. 32d. That figure does not show any crossover in the new grid after the large deformation. The new grid not only keeps the pointtopoint 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. 32e). More details about the master and slave concept can be found in Refs. [29, 56]. CHAPTER 4 ACTIVE FLOW CONTROL WITH AN ADAPTIVELYSHAPED WING 4.1 Introduction The main objectives for fluid dynamists and aeronautical engineers throughout the history of modern flight have ahlvb, 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 lowReynolds numbers which is usually norotiously accompanied with laminar boundary separation, transition, and low lifttodrag 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 piezoactuated flap mounted on a low Reynolds number wing have been studied both experimentally [19] and computationally [30, 56]. Gad elHak [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 piezoactuated 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 multiblock, moving grid technique to handle the geometric variations in time, using twoequation turbulence closures, and a pressurebased 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 Setup The wing has moving flaps on the upper surface as shown in Fig. 41. 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 setup, 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 41: Moving flap wing experimental setup in the windtunnel In this chapter, a number of twodimensional and threedimensional cases, in cluding different grid densities, with and without the gap, will be studied. For 3D cases, an effort has been made to consider the effect of the gap between flaps as shown in Fig. 42. The original LaunderSpalding kE model [52] with the wall function, as well as a lowReynolds number version kE model [46] are compared. Linear Chords  Figure 42: Flap wing with refined three dimensional bending model. 4.3 Numerical Study 4.3.1 Grid Sensitivity Analysis A number of 2D grids with different distributions and sizes have been tested. Then a reference 2D grid is selected by ensuring that satisfactory solution quality is obtained on a most economic grid size. Then, the 3D grid is established by extending the reference 2D 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 multiblock 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 Reynoldsaveraged 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 lowReynolds number version of the ke model seems more suitable [56]. Fig. 43 compares the timeaveraged 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 lowReynolds 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. 44 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 counterclockwise direction on the left of the gap when viewed from above. The gap also affects the overall pressure distribution on the wing. Fig. 45 compares the pressure distributions obtained with and without gap between flaps. The lowReynolds 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 tineaveraged 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 3D 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. 46 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 45: Comparison of timeaveraged 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 46: 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 104 + 1.6 x 104 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. 47 compares the timeaveraged pressure coefficients between the two flap ping amplitudes. For the higher flapping amplitude, the timeaveraged 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. 48a. In our computations with the smaller flapping amplitude, the vortex core is anchored at the flap tip, exhibiting two corotating cores as the flaps move upwards, as shown in Fig. 48b. On the other hand, with the higher amplitude, vortex shedding appears for cases with or without 0.5 0 Disp=l.0x 103 in SDisp=7.3 x 103 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 47: Pressure coefficient at different actuation amplitude without gap between flaps. Actuation frequency: 500Hz. the gap. Fig. 49 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 48: 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 lowReynolds number formulation of the ke turbulence model produces longer circulation zone. The gap between the flaps affects the location of the reattachment point, and creates a threedimensional 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 8r/6. (e) Phase angle a = 11T/6. Figure 49: 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 threedimensional bodies. Until recently, membrane analyses were primarily relied on trial and error. Green and Adkins [25] are the first to derive a general nonlinear 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 straindisplacement 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 Fopplvon 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. Twodimensional 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, twodimensional, 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 threedimensional membrane model introduces several complicating factors not encountered in twodimensional analyses. First, unlike the  il equation", the tension in a threedimensional 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 threedimensional 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 rubberlike 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 threedimensional membrane model is proposed for the mem brane dynamics. The membrane material considered obeys the hyperelastic Mooney Rivlin model [64]. We use the GreenLagrange strain tensor for the description of large strains. The dynamic response of the membrane is described by a system of secondorder ordinary differential equations. A finite element method with triangu lar elements for spatial discretization is utilized. For time integration, the implicit WilsonO method is used. 5.2 Structural Solver 5.2.1 MooneyRivlin 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 MooneyRivlin material [64], and the model is called the MooneyRivlin model. The MooneyRivlin 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 stressstrain relationship is writ ten as 8W S = pC + 2 (5.4) where S is the second PiolaKirchhoff stress tensor, p is the hydrostatic pressure. If the membrane is incompressible, then the stressstrain relation can be further simplified as the following S pC1 + 2 [(cl + c21)I c2C] (5.5) where I is 3by3 identity matrix. In the membrane theory the stress field is essentially treated as twodimensional, 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 S330, 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 GreenLagrange strain tensor, S is the second PiolaKirchhoff 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. 51, a typical triangular element, e, is defined by nodes, i, j, k, and three straightline 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. 51, in local coordinates, the triangular element lies in the plane of k' y=x+u Y k x 1 z Y X Figure 51: 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= T1. 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 GreenLagrange 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. 51. Following Green and Adkins [25], the GreenLagrange 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 GreenLagrange 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 twodimensional, 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 9i2u 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) 611, 6FLxy where Bo and BL are the straindisplacement 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 9by9 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. 52, 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 52: 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 3by3 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 vv to integrate the structural dynamics equations in time; among these we utilize the WilsonO method [113]. The WilsonO method has a controllable algorithmic numerical dissipation, it is a onestep 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 WilsonO 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 MooneyRivlin model. There are two con stitutive parameters in this model, namely cl and c2. To simplify our computation, we define a nondimensional MooneyRivlin 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 nondimensional time r, and the nondimensional 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 semianalytical 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 semispherical surface. In our work we choose At5 x 104 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 106 for the same problem. All computations are performed using the doubleprecision. In Fig. 53 we com pare the finite element results with the semianalytical solutions. Except for one case in Fig. 53b, all numerical solutions match well with the semianalytical solutions. The loss of linearity of oscillations is well reproduced. The mismatched case shown in Fig. 53b 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 MooneyRivlin constant a is, the higher pressure the membrane can sustain. Two profound issues should be considered in the use of the WilsonO method. First, the WilsonO 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 WilsonO method due to its excessive dissipation. As shown in Fig. 54, a larger time step tends to cause more numerical damping. In Fig. 55 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 53: Dynamic inflation of a Mooney spherical membrane: (a) a = 0; (b) a 0.1; (c) a = 0.25. Solid line: Semianalytical solution; Circle: Numerical solution. A t=2.0e4 1.375 $1 ;~ < 1.25 1.125 10 15 20 25 0 5 10 15 20 25 Figure 54: The effect of time step to the numerical damping for p*=0.5 and a=0. 4.7 0) t4.9 / 0 LU. 5.1 / 3.8 3.6 3.4 Time Step (login) 3.2 Figure 55: 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.0e4 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 wingbased 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], twodimensional flow based computations [36, 86, 92, 93], or simplified threedimensional 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 threedimensional NavierStokes equa tions for incompressible flows, and the structural dynamics is governed by a nonlinear dynamic membrane model. The NavierStokes equations are solved using a pressure based method on a structured and multiblock grid. The membrane material follows the hyperelastic MooneyRivlin 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 WilsonO 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 stepbystep 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 vr separation and reattachment, tip vortices, and unsteady phenomenon at high angle of attack. In the membrane wing study, we study the selfinitiated 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, smallscale 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 timedependent equations in a pseudotime. 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 threefield problem: the fluid, the structure, and the dynamic mesh that is represented by a pseudostructural system. In their work, the arbitrary Lagrangian Eulerian (ALE) finite volume scheme is emploiy, 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 fluidstructure 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 onedimensional 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 fluidmembrane dynamics, iterations between the fluid and structural solvers should be done at each time instant. A stepbystep 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 NavierStokes 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. 11 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 illconditioned. 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. 61. 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. 62. 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 61: Schematic wing with three battens on each half wing. (9c, 5c, 6c B Membrane Wing C (6c, 5c, Figure 62: Computational setup and boundary conditions. at the coarse level around the wing is shown in Fig. 63. The results from the fine grid are used as the reference. The comparison among these three grid systems is summarized in Table 61. The coarse grid solution underpredicts the lifttodrag 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 61 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 63: 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 ir velocity Table 61: Grid refinement effects on the computed aerodynamic coefficients. GRID Dp(N) L(N) DF(N) CL/CD Coarse 1.8x 10 6.24E2 5.10E1 9.70E3 7.06 Intermediate 6.7x105 6.16E2 5.08E1 9.22E3 7.16 Fine 2.3x 106 6.20E2 5.04E1 8.80E3 7.10 profile is shown in Fig. 64. 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 64: 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. 65. 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 lvr 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 5x105s. However, for our computations, the difference between the steady and the timeaveraged lift and drag coefficients, are found to be less than 1 over all the studied cases. Fig. 65 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.rr 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 65: 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+ LA (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 twofolded 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. 66 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 anticlockwise rotations; it seems that the flow is leading from the lower surface to the upper surface. Fig. 67 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. 68 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. 68, 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. 69a where the spanwise pressure coefficient on the upper wing surface at x/c=0.4 is plotted. At Figure 66: Streamline and vortices for rigid wing at a= 39. The vortical structures are shown on selected planes. o 10060 rp Figure 67: 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 68: 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. 68). 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. 69b it can be seen that the pressure is almost uniform in the spanwise direction even at high angles of attack. Fig. 69 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. 69 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 69: 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 "deadair" region, where the recalculating velocity is significantly less than the freestream velocity and the flow is virtually stationary. The separated flow forms a freeshear li.r, which is highly unstable. The separated shear 1ir 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. 610. 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. 610, the separation oc curs at I'. of the chord from the leading edge, which is followed by a long "deadair" zone before it reaches a maximal reverse velocity of 0.26U. In the "deadair" 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. 611. 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 610: 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. 612, 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. 65a, 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 IOO (d) 39 1.45 0.50 [] 0.45 1.40 (e) 51 Figure 611: Normalized chordwise velocity u/U contours at different angles of attack. Figure 612: Vortex shedding at a51 Figure 612: Vortex shedding at a=51. Unsteady Phenomenon at High Angles of Attack Both vortex shedding and boundary 1vr 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 timeaveraged 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 lifttodrag ratio are found to be less than 1 (Fig. 613). 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 613: 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. 614 shows that at a= 6 the time averaged pressure coefficient 01 08 0 1 04 1 06 12 0 02 04 06 08 1 0 02 04 06 08 1 x/c x/c Figure 614: 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. 611 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. 611. 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. 615, 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. 616). 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 615: Comparison of chordwise velocity profiles on the rigid wing between the steady and timeaveraged 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 616: Comparison of the chordwise velocity profiles on the rigid wing be tween the steady and timeaveraged 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. 617 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. 618 and Fig. 619, 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. 620a 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. 620b. 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 617: 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 618: 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 619: 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. Selfinitiated 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. 621 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 620: 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. 617, and the estimated dom inated frequency is around 120 Hz (Fig. 622). 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. 623. Fig. 623a 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.04 N S0.02 z t=o t=3 x 103 t=7 v 103 1=10 10 I 1,. 10 l=1. 10 0.02[ 0.04 0 0.2 Figure 621: 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 103 7 6 _5 4 2 50 100 150 200 250 Frequency (Hz) 300 350 400 Figure 622: Spectrum analysis of the trailing edge point vibration for membrane wing at a=6. 83 v/U contour at t=3.0 x 103 v/U contour at t=7.0 x 103 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 120 1 '0 12 12 00592 6 4 2 0 2 4 6 6 4 2 0 2 4 6 Spanwise (cm) Spanwise (cm) Figure 623: 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. 623b 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 62, 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. 624 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 62: 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. 625. As we can see from Fig. 625a, 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. 625b 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. 626. 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. 624 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 624: 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 625: 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. 628. 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 626: 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 627: 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 628: 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 twodimensional analysis tool to test the performance of different airfoils. Based on the tests, one air foil shape is then chosen to form a threedimensional wing. The wing shape is then tested with flight. In this trial and error process, the threedimensional 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 lifttodrag 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 twodimensional fluid flow solutions based on coupled thin liv r and inviscid flow models. A schematic baseline wing shape is shown in Fig. 61 in C'! lpter 6. The baseline shape has three carbonfibermade battens that are labelled also in Fig. 61. 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 lifttodrag 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 ycoordinates 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. 71. 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 71: 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. 7la) 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 71. Table 71: 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 