Membrane and adaptively-shaped wings for micro air vehicles

Material Information

Membrane and adaptively-shaped wings for micro air vehicles
Lian, Yongsheng ( Author, Primary )
Publication Date:
Copyright Date:


Subjects / Keywords:
Aerodynamics ( jstor )
Aircraft wings ( jstor )
Airfoil camber ( jstor )
Coordinate systems ( jstor )
Low reynolds number ( jstor )
Rigid wings ( jstor )
Shape optimization ( jstor )
Turbulence models ( jstor )
Velocity ( jstor )
Wing flaps ( jstor )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Lian, Yongsheng. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:


This item has the following downloads:

Full Text







Copyright 2003


Yongsheng Lian

To my parents and to my wife Lihui, for their love and support


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


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.


ACKNOWLEDGMENTS ................... ...... iv

ABSTRACT .. .. .. .. ... .. .. .. .. .. .. .. ... .. .. vii


1 INTRODUCTION ........................... 1

1.1 MAV and Membrane Wing Aerodynamics ........ .. .. 1
1.2 Motivations and (C! iII. 'ges ......... ........ ... 4
1.3 Outline ............ ...... ........ ..... 5

F L O W S . . . . . . . . . 8

2.1 Governing Equations and Fluid Solver ............... 9
2.1.1 Governing Equations for Fluid Flow ............. 9
2.1.2 Fluid Flow Solvers ........... ........... 12
2.1.3 Geometric Conservation Law ....... ......... 12
2.2 Turbulence Models . . ..... ....... 13
2.2.1 Reynolds Equations and k-E Turbulence Model ...... 13
2.2.2 Boundary Conditions for Turbulence Models . ... 15
2.3 Laminar-turbulent Transition ............. .. .. .. 18

3 MOVING GRID TECHNIQUE .................. ... .19

3.1 Perturbation Method .................. ..... .. 20
3.2 Master/Slave Concepts .................. ... .. 21
3.3 Numerical Tests .................. ......... .. 26


4.1 Introduction ............. . . ...... 27
4.2 Experimental and Computational Set-up . . ...... 28
4.3 Numerical Study .................. ...... .. .. 29
4.3.1 Grid Sensitivity Ain -i ........... .. ... 29
4.3.2 Assessment of the Turbulence Model . . 30
4.3.3 Effects of Gaps between Flaps ............ .30
4.3.4 Effects of Flapping Amplitude ............ .32
4.4 Conclusions ............... ........... .. 34

5 MEMBRANE STRUCTURAL SOLVER ........ .......... 36

5.1 Introduction ................... ....... 36
5.2 Structural Solver ................... .... 38
5.2.1 Mooney-Rivlin Model ................ .. .. 38
5.2.2 Principle of Virtual Work .................. .. 39
5.2.3 Triangular Membrane Element . . ..... 40
5.2.4 Green-Lagrange Strain Tensor ............... .. 42
5.2.5 Internal Forces .................. .... 43
5.2.6 External Forces .................. .... .. 45
5.2.7 Dynamic Equations for the Membrane . . .... 46
5.2.8 Solution of the Dynamic Equations . . ..... 47
5.3 Numerical Test of Membrane Model ................ .. 48


6.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




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



Yongsheng Lian

August 2003

C('! I: Wei Shyy
Major Department: Mechanical and Aerospace Engineering

Micro air vehicles (ilAVs), with wing span of 15 cm or less and flight speed

around 10 m/s, have many applications in both civilian and military areas. The

Reynolds number based on the given parameters is around 104, which often yields

insufficient lift-to-drag ratio. Furthermore, one expects the unsteady effect to be

noticeable for such flight vehicles. The flexible wing has been demonstrated to exhibit

favorable characteristics such as passive adaptation to the flight environment and

d. 1 .1 stall.

The present study focuses on developing computational and modeling capabilities

to better understand the MAV aerodynamics. Both flexible wings, utilizing mem-

brane materials, and adaptively-shaped wings, utilizing piezo-actuated flaps, have

been studied. In the adaptively-shaped wing study, we use piezo-actuated flaps to

actively control the flow. We assess the impacts of the flap geometry, flapping am-

plitude, and turbulence modeling on the flow structure with a parallel experimental

effort. The membrane wing uses a passive control mechanism to delay the stall angle

and to provide a smoother flight platform. Our study focuses on the mutual interac-

tions between the membrane wing and its surrounding viscous flow. We compare the

lift-to-drag ratio and the flow structure between the flexible wing and the correspond-

ing rigid wing. We also investigate the aerodynamic characteristics associated with

the low Reynolds number and low aspect ratio wing. To assist our study, we propose

an automatic and efficient moving grid technique to facilitate the fluid and structure

interaction computations; we also present a dynamic membrane model to study the

intrinsic large deformation of the flexible membrane wing. Solutions obtained from

the three-dimensional Navier-Stokes equations are presented to highlight the salient

features of the wing aerodynamics.

Besides the aerodynamic study, we also perform shape optimization to improve

the membrane wing performance. Since direct optimization of a membrane wing is

too time consuming to be practical, we optimize a surrogate rigid wing model based

on an integrated optimization algorithm, which consists of a N ,',i i-Stokes solver, an

automatic grid generation tool, and a gradient-based optimizer. Then, we assess the

membrane wing performance based on the outcome from the surrogate model. Our

numerical results confirm that the membrane wing exhibits consistent improvement

in the lift-to-drag ratio with the surrogate model.


1.1 MAV and Membrane Wing Aerodynamics

Micro air vehicles ( \!AVs) with a maximal dimension of 15 cm or less and a

flight speed around 10 m/s are of interest to both military and civilian applications.

Equipped with a video camera or a sensor, these vehicles can perform surveillance,

reconnaissance, targeting, and biochemical sensing at a remote or otherwise hazardous

location [83]. With the rapid progress made in structural and material technologies,

miniaturization of power plants, communication, visualization, and control devices,

several groups have developed successful MAVs [24, 34]. As its size is reduced, MAV

gains favorable scaling characteristics, including low inertia and reduced stalling speed

[83]. However, in terms of aerodynamics, control, range, and maneuverability, many

outstanding issues are still present. First, an MAV flies in a low Reynolds number

regime (104r105) due to its low flight speed and small dimension. Such a flight

environment is often accompanied by laminar boundary 1.-r separation, transition,

and low lift-to-drag ratio. Second, an MAV wing typically has a low aspect ratio,

which causes strong vortical flow structures and increases the induced drag. Third, an

MAV is susceptible to rolling instabilities, which becomes even more serious because

of the existence of tip vortices. Fourth, an MAV is sensitive to wing fluctuations,

which can be comparable to MAV's flight speed and makes both the instantaneous

flight Reynolds number and angle of attack vary substantially [83].

In the development of MAVs, two approaches exist in the MAV wing design:

flapping wings and fixed wings. In the flapping wing area, pioneering work was done

by Weis-Fogh [112], and Lighthill [57]. Recent experimental and numerical studies

in this aspect were given by ('!ii i-i,, i and C'!i ,I.:i varthy [8], DeLaurier [11], Smith

[90], Dickenson et al. [13], Ellington [17], Jones and Platzer [43, 44, 45], Katz [50],
Kamakoti et al. [48], Liu and Kawachi [59], Vest and Katz [106, 107], Walker [109],

and Wang [110]. Shyy et al. [83] gave an extensive review of flapping and fixed wing
characteristics. Templin [96] presented the flapping wing spectrum of flight animals.

Our study focuses on a fixed wing. Specifically, we will study a flexible membrane
wing for MAV applications. Figure 1-1 shows a 15 cm MAV with a flexible wing
designed by Ifju and coworkers [34]. The MAV has a flight speed between 20 and 40
km/h. Powered by an electric motor, it can fly for about 15 minutes. With a build-in
video camera and a transmitter, the vehicle has a total mass around 40 g. The MAV

has a rigid fuselage and leading edge spar which are made of laminate material. On

each side of the wing there are three unidirectional carbon fiber battens. The spar
and battens are covered with a latex membrane material. Overall, the wing has a

span of 15.2 cm, a root chord of 13.2 cm, a mean chord of 10.5 cm, and a wing

area of 160 cm2, which result in a low aspect ratio of 1.4. The wing has a variable
camber decreasing from 7.5'. at the root to 2' at the tip. In this configuration, the
maximum camber is located at about '7'. chord at the root.


Figure 1-1: MAV with membrane wing.

Flexible wings have different aerodynamic performance from rigid wings. One

advantage of flexible membrane wings is that they facilitate passive shape adaptation

and result in a d. 1,. .1 stall [83, 111]. Fig. 1-2, adopted from Waszak et al. [111],

compares the lift curves versus angles of attack for rigid and membrane wings with

configurations similar to those shown in Figure 1-1. Under modest angles of attack,

the rigid and membrane wings demonstrate a similar lift curve slope of 2.9, with the

rigid wing having a slightly higher lift coefficient. However, the membrane wings have

a much higher stall angles of attack than the rigid wing, which is a key element in

enhancing the stability and agility of MAVs. The rigid stalls at 24" while the flexible

wings have stall angles between 30" and 49", which are similar to that of much lower

aspect ratio rigid wings (AR=0.5 to 1.0). But, the very low aspect wings exhibit

lower lift curve slopes of 1.3 to 1.7 instead of 2.9. Hence, flexible wings appear to

combine the desirable performance of rigid wings with higher and lower aspect ratios

by exhibiting stall behavior similar to that of rigid wings with aspect ratio of 0.5 to

1.0 and the lift generating capability similar to rigid wings with aspect ratio of 2.0

[111]. Another advantage of flexible membrane wings is that they can adapt to the

wind gust and provide smoother flight platforms, which are verified by wind tunnel

experiments [84].

0 -- Rigid (Graphite)
-- --- 6-Batten (Monofilm)
-- -- 2-Batten (Latex) -
2.2 2'I
C 1.5

0.75 -

0 .0 L
-10 0 10 20 30 40 50
Angle of Attack (deg)

Figure 1-2: Lift coefficient versus angle of attack. (From Waszak et al. [111])

1.2 Motivations and ('!C iII. iges

Though experiments have provided information on the membrane wing aero-

dynamics, to fully understand the membrane aerodynamics and to appreciate the

membrane wing's passive control mechanism, we need to perform a detailed numer-

ical study. A numerical study of the membrane wing brings two challenges: the low

Reynolds number and the low aspect ratio condition, and the interaction between the

membrane wing and its surrounding viscous flows.

First, the low Reynolds number condition presents tremendous challenges on

MAVs study. In our study, the chord Reynolds number of the MAV is around 9x 104.

In the range of Reynolds number of 104W 106, complex flow phenomena often occur.

The laminar boundary 1lv- r separation, transition, and reattachment usually coexist

on the upper wing surface. An accurate predication of the transition from laminar

to turbulent state is important to evaluate the wing performance. However, a transi-

tion model for general application is not yet available. Furthermore, low aspect ratio

wing is usually accompanied by complicated vortical flow structures. The vortical

flow structures introduce unsteady phenomena. Valuable information regarding low

Reynolds number and low aspect ratio wings was reviewed by Lissaman [61] and Tani

[95], as well as addressed in several proceedings [65, 66]. Gad-el-Hak [20] discussed
approaches to prevent/alleviate flow separation on MAVs with microelectronicsme-

chanical systems (\!I \!S). Lian and Shyy [53] performed a detailed numerical study

of the flow separation on the MAV wing and its impacts on aerodynamic performance.

Second, a membrane wing exhibits self-initiated vibrations even under steady

state freestream condition. Waszak et al. [111] conducted wind tunnel experiments

and found that the membrane wings exhibited vibration at 0(102) Hz in a steady

freestream. Similar observation was reported numerically by Lian and Shyy [53]. Such

vibrations and the associated shape deformation change the pressure distribution on

the membrane wing, which, in turn, affects the membrane dynamics. This recipro-

cal action leads to a fluid and structure interaction problem. Aeroelastic analysis

method coupling a computational fluid dynamics (CFD) solver with a computational

structural dynamics (CSD) solver provides an effective tool for the membrane wing

study. Although both CFD and CSD have achieved great success individually, a cou-

pled computation is still a challenging task [80]. To date, the study of a membrane

wing and fluid flow interactions is limited. Jackson and C('! i-I1 ,, [38] adopted a three-

dimensional potential flow-based solver and a membrane wing model to analyze the

aeroelastic behavior of marine sails. Smith and Shyy [92] and Shyy et al. [88] pre-

sented a computational approach to model the interaction between a two-dimensional

flexible membrane wing and surrounding viscous flows. A combined numerical and

experimental ,a i1, ~i-; of a flow over a two-dimensional flexible sail was conducted by

Lorillu and Hureau [62]. In the coupled fluid and membrane wing interaction study,

we need a fluid solver to accurately capture the transient phenomenon of the fluid

flow surrounding the membrane wing; we need a structural solver to predict the in-

herent nonlinear behaviors of the membrane material. To couple these two solvers,

we also need to handle different time characteristics associated with each solver. Be-

sides these, we need a moving grid technique to facilitate the fluid and structure

interaction problem. An interface technique is also required to exchange information

between these two solvers.

1.3 Outline

This work study both the membrane wing and adaptively-shaped wing. First, we

will develop a computational capability to study the fluid and structure interaction.

Second, based on this capability, we will study the piezo-actuated flapping wing the

membrane wing. Third, with the gained knowledge, we will optimize the membrane

wing performance with modern optimization technique. The outline of the current

work is as the following:

In ('!C ipter 2 we present the fluid flow solver and physical models of the fluid

dynamics. The Navier-Stokes equations for incompressible flow on curvilinear coordi-

nates are employ, l along with a two-equation k-e turbulence model. In this chapter,

the fluid solver is considered under the moving boundary framework. A brief in-

troduction is given to the geometry conservation law [100] that p1 .,i- an important

role in moving boundary problems. Even though the present work does not employ

laminar-turbulent transition model, a brief review is offered.

In ('!i lpter 3 we devote to moving grid technology, which is used by both the

coupled fluid-structure study and the optimization. A robust, efficient and effective

moving grid technique is proposed, which uses transfinite interpolation-based pertur-

bation method in single block grid [56, 77], and employs the master/slave concept

[29] for a multiblock grid.

In ('! Ilpter 4 we investigate the active flow control with adaptively-shaped flaps

for a low-Reynolds number wing. Results of the interaction between piezo-actuated

flaps and surrounding flow are presented.

In C'!i lpter 5 we focus on the structural solver, more specifically, a nonlinear

dynamic membrane model. Detailed derivation of the membrane model, with an

emphasis on the nonlinearity as well as the time integration will be presented. The

membrane model will be assessed using well defined test problems.

In C'!i lpter 6 we propose a coupled fluid and structure interaction algorithm. The

coupling is realized through an interface technique and time synchronization. The

numerical implementation of geometry conservation law is given and tested. The low

Reynolds and low aspect ratio wing aerodynamics, including the laminar boundary

l -.-r separation, tip vortices, and the unsteady phenomenon are investigated. Com-

putational results will be presented to elucidate the membrane wing aerodynamics.


In ('!i ipter 7 we optimize of a flexible membrane wing with a gradient-based

method. Since a direct optimization of a membrane wing requires substantial com-

putational resource and time, a rigid wing is used as a surrogate. The membrane

wing performance is then evaluated based on the surrogate model.

In brief, the present work is motivated by, first, the practical need to understand

the aerodynamics of membrane and adaptively-shaped wings, second, our interests

in micro air vehicles. Original work include a moving grid algorithm for multi-block

grid, a dynamic membrane model, numerical study of interaction between fluid and

a flexible wing, and shape optimization of a flexible wing.


Our study investigates the aerodynamics of a membrane wing as well as adap-

tively shaped wings under low Reynolds number conditions. During flight, the mem-

brane wing undergoes shape change under the external forces; meanwhile, its shape

variation affects the pattern and structure of its surrounding fluid flow. The mem-

brane wing vibration was observed both experimentally [111] and numerically [53, 55].

In the adaptively shaped wing study, the piezo-actuated flaps vibrated with different

frequencies and amplitudes. The flap vibrations cause unsteady phenomena like vor-

tex shedding [30, 56]. It is desirable to have a computational capability to accurately

capture the transient behavior of fluid flow and structural dynamics.

Even thought steady flow computations have achieved significant success in their

development and applications, little attention has been given to the systematic analy-

sis of unsteady computations [79, 80]. Several factors are responsible for this situation.

First, steady CFD is still the mainstream aerodynamic analysis, and the techniques

used to enhance the performance of steady CFD computations may not be suitable

to unsteady computations. For example, in steady computations, a variety of tech-

niques have been proposed to eliminate the low frequency transients to accelerate

the convergence to steady state; however, low frequency phenomena are intrinsic to

many fluid and structure problems, and such techniques will be inappropriate to the

unsteady analysis [80]. Second, the applicability of the turbulence model to unsteady

flows is not well understood [79].

However, numerous approaches have been proposed for unsteady flow computa-

tions. Pulliam [74] incorporated subiteration to improve time accuracy of conven-

tional implicit schemes. Jameson [39] and Rumsey et al. [79] developed a dual time

stepping technique in the context of multigrid methodology to enhance the efficiency

and accuracy of traditional factored schemes. Issa [35], and Oliveira and Issa [71] pro-

posed the Pressure-Implicit with Splitting of Operators (PISO) method to improve

the convergence. For relevant information of these and some of the most frequently

adopted approaches, see Shyy and Mittal [85].

The SIMPLE method [72, 82] and the PISO method [35, 71, 98] are two popular

methods to solve the three-dimensional Navier-Stokes equations for incompressible

flows. Both methods belong to the pressure-based approach [82] which devises an

artificially derived pressure (correction) equation by manipulating the mass continuity

and momentum equations.

Three-dimensional N ,i.-, r-Stokes equations for incompressible flow in curvilinear

coordinates are adopted in our study. Two methods are emplo-, '1 to solve the gov-

erning equations, one is the SIMPLEC method [102] that is belong to the SIMPLE

family, and the other is the PISO method [35]. In this chapter we briefly introduce

these fluid solvers under the moving boundary framework. For detailed descriptions

of these methods, we refer to Refs. [35, 72, 82, 98, 102]. To solve the moving bound-

ary problems, we also emphasize the importance of the boundary conditions and the

geometric conservation law [100]. An introduction of turbulence models, with an

emphasis on the k-e turbulence model is presented. Though an accurate transition

model for general application is beyond the scope of our study, we offer a brief review

of the transition models at the end of this chapter for completeness purpose.

2.1 Governing Equations and Fluid Solver

2.1.1 Governing Equations for Fluid Flow

The three-dimensional N ,1-i. r-Stokes equations for incompressible fluid flows in

curvilinear coordinates [82] in strong conservative form can be written as the follow-

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(

[ (j11iu + qi1- + q13u()]

+ P [(q21 + q__.. : + 923U()]

+ [a (qU31+ 32U+ "+ )

- [(fiP) + (f21P) + (f31p)

+ P (q21V + q._ + q23Vc)]

+ P, (31q +(I + q33V()]

- [(f12P)+ (f22P) )+ (f32P)
+ G2a(, ).J, (2.2)

P (qdjl(w + l-12 + q13w()]
+ [ (q21W+q_, + q23UW)]

+ P[(q3iw + 4 _: ,, + q33')

[(fl1P) + (f2P)+ (f32P)]

+ G3(~ )J, (2.4)
+ [( + 41 _,+ Y1: )1

+ [ (f21p + q(f'p) + i(fq)2

+ G3((, ).J, (2.4)

where u, v, and w are the velocity components in the x-, y-, and z-directions, re-
spectively, p is the viscosity, G1, G2, and G3 are the body force components per unit
volume in the x-, y-, and z-directions. The coefficients appearing in the equations

are expressed as follows:

f2 + f2 + f 12
211 122 13, 912

f221 + f222 + f223, 13

f321 + f322 + f323, 923

fll21 + f12f22 + f13f23,

fll31 + f12f32 + f13f33,

f31f21 + f32f22 + f13f23.

The expressions for fj, i 1, 3, j=1, 3, can be found in [82]. The Jacobian transfor-

mational matrix is defined as follows:


The determinant of the Jacobian matrix is expressed as

J = x0yxz + xyzx, + xyjzz xyZzq xy,'z xyzq,.


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


where, for example, x is the grid velocity component in the x-direction that can be

approximated by the first order Euler method,

x Xo


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.


9(x, y,z)
J 9= ,C

2.1.2 Fluid Flow Solvers

The SIMPLEC method [102] is a variant of the SIMPLE method originally pro-

posed by Patankar [72] in Cartesian coordinates, and its extension to curvilinear

coordinates was discussed by Shyy [82] to accommodate the complex geometry com-

putation. The SIMPLEC method uses coordinated under-relaxations for the momen-

tum and pressure corrections to improve the convergence that is intrinsically slow in

the original SIMPLE method. The SIMPLEC method was also used to solve mov-

ing boundary problems [53, 88]. Contrary to the SIMPLE family method, the PISO

method is a non-iterative approach to handle the pressure-velocity coupling which

splits the process of solution into a series of predictor and corrector steps. The PISO

method is generally more efficient for transient flow computations [35, 98].

To solve the Navier-Stokes equations, we require proper boundary conditions for

the computations. At the interface between the fluid and the solid structure, the

following boundary conditions are applied,

Xf = x, on F, (2.10)

Oxf ax8
t at on F. (2.11)
9t 9t '

Eq. (2.10) expresses the compatibility between the wet fluid surface Xf and the struc-

ture surface x, at the interface F; Eq. (2.11) requires that the fluid velocity at the

interface F matches the surface velocity of the structure. These boundary conditions

are provided by computing the displacement of the solid structure.

2.1.3 Geometric Conservation Law

When body-fitted curvilinear coordinates are used in the computation, a trans-

formation maps a physical flow region (x, y, z) onto a uniform computational space

( T1, () where conservation laws are carried out. To facilitate the coordinate trans-

formation, the Jacobian matrix J is introduced. The Jacobian matrix is defined as

in Eq. (2.6). The determinant of J represents a volume element in the transformed

coordinates. In a moving grid problem, the computational grid needs to be updated

with time; meanwhile, the Jacobian J needs to be updated simultaneously. Specific

procedures are required to compute the effective value of J; otherwise, errors arise due

to the inconsistency in evaluating geometric quantities and the adding of non-physical

quantities. As pointed by Thomas and Lombard [100], in the updating process, the

following geometric conservation law (GCL) needed to be satisfied,

d Jdidqd(= (V W,)ddqd(, (2.12)

where V is the volume bounded by a close surface S, W8 denotes the local velocity

of the boundary surface S. Thomas and Lombard [100] proposed to evaluate J while

maintaining the geometric conservation law by setting p=1, V=0 from the continuity

equation, i.e.,

Jt + (JW) + (Jt)k, + (J(c = 0. (2.13)

Implications of the geometric conservation law with fluid-structure interaction prob-

lems were discussed by Shyy et al. [88].

2.2 Turbulence Models

2.2.1 Reynolds Equations and k-E Turbulence Model

For simplicity, the Reynolds averaged momentum and continuity equations for

incompressible flows are presented in Cartesian coordinates

au, a(uJ U) aTj (2.4)
at + x (2.14)xj

0, (2.15)

where Ti is the stress tensor, repeated indices in any term indicate a summation over

all three values of the index, and Ui is the Reynolds averaged mean velocity. The

stress tensor in a turbulent flow may be written as

T, = -P6,j + 2pSj ,i, i,. (2.16)

6ij is the Kronecker delta, which is equal to one if i= j and zero otherwise. P is the

mean hydrodynamic pressure and p is the dynamic viscosity. The mean rate of strain

Sij is defined by
1 9U, au.
Si ~ + a ) (2.17)
2 Oxj Oxi

ui and uj are the fluctuation velocities in the turbulence flow which contribute to the

mean stress by

T j -7, ,., (2.18)

Tij is the Reynolds stress tensor which is symmetric. The diagonal components of Tij

are normal stresses which make up the turbulent kinetic energy k

k ui (2.19)

In many flows, these normal stresses contribute little to the transport of mean mo-

mentum. The off-diagonal components of ij are shear stresses which pl ', a dominant

role in the mean momentum transfer by turbulent motion.

Eq. (2.14) contains nine unknown components of ij (of which six are independent

of each other) besides P and the three components of Ui. Therefore, we need more

equations for ij to close the system of equations. A turbulence closure model serves

that purpose.

Different versions of turbulence closure models have been proposed. The success

of such models is largely depend on the applications of the wall functions that relate

surface boundary conditions to fluid points away from the boundaries. Thereby we

avoid the problem of direct modelling the viscosity influence. Here we use the wall

function treatment by Launder and Spalding [52] as one example to illustrate the

closure problem. The governing equations for standard k-E models can expressed in

the form shown below,

a(pk) a(pkUi) 9 ( it 9k
t+ ) ax l ) + 2tSijSj pE, (2.20)
at axi 9xi Jk xi

9(pE}) 9(pEU) a Pt E E E2
(p + +( 2CIptSijSij C2p. (2.21)
at + xi ax"i aE xi k k
The eddy viscosity pt is defined as

pt =pC-, (2.22)

where C, is a dimensionless constant. Sij is the mean rate of strain.

These equations have five adjustable parameters, CO, Jk, E, C.1, and C.2. Here,

for brevity, we only give the coefficients for standard k-E model -ir_-'-. -1. I by Launder

and Spalding [52]. These parameters are valid for a wide of range of turbulent flows.

C, 0.09, rk = 1.00, a 1.30, C1 = 1.44, CE2 1.92. (2.23)

To close the system of equations for turbulence computation, we need an additional

Boussinesq relationship that relates the Reynolds stress to the mean rate of strain

and turbulence kinetic energy,

-, Uj = ( + O ) -pkJij 21 tSij pkbij. (2.24)
axj axi 3 3

2.2.2 Boundary Conditions for Turbulence Models

The elliptic characteristic of the k and E model equations gives rise of bound-

ary conditions. Usually we assume zero gradient at the outlet boundary. At inlet

boundary, if there is no measurement of k and E at disposal, one way is to make an

approximation based on the freestream velocity U,, turbulence intensity Ti, and the

characteristic length scale L. In our study we adopt the following formulas to define

the freestream boundary conditions for k and E,

k (UoTi)2, (2.25)
c 4 2 (2.26)

The turbulence intensity Ti is usually set between the range of 0.02 and 0.05. Another

way is to adjust the turbulence Reynolds number Ret=pUL/pft. For practice we

usually set the turbulence Reynolds number in the range of 100 and 500.

The solid wall boundary conditions are more complicated. Near the wall bound-

aries, the local Reynolds numbers are low and the viscous effects are important.

However, the k-E turbulence model that we mentioned above is based on the high-

Reynolds-number assumption in which the turbulent stress forces are dominated.

Moreover, close to the wall, extremely fine grid is usually required to capture the

steep variations in flow profiles to solve the Navier-Stokes equations directly. The

CPU time based on such a fine grid is almost beyond the ability of current facilities.

In order to solve these difficulties, two approaches are commonly emploi- l1 one is

the wall function treatment [52], the other is the low-Reynolds-number versions of

the k-E [46, 73].

The wall function treatment is popular in practice because of its simplicity and

capability for complex geometries. Under the high Reynolds number condition, flow

can be largely divided into two 1-v. rs: the outer 1-v -r and the surface -1i.-r. In the

outer l-?v-r the Reynolds stress is dominated; while in the surface l-?v-r both the

Reynolds stress and the shear stress are important. In the outer l1.v-r, the fluid field

is influenced by the retarding effect of wall through the value of the wall shear stress,

but not by the viscosity itself.

In the wall function treatment, velocity profiles have different descriptions in the

out and surface 1-~v.i--. The velocity in the out -1~ .-r is governed by the velocity-defect

law expressed as
max U g( ), (2.27)
U, 6
where 6 is the boundary l-i-r thickness, Umax-U is the velocity deficit, g is an

unknown function with the order of 1. The velocity in the surface l-v-r is given by

the law of the wall

= f(y ). (2.28)

Here y+ is a non-dimensional normal distance from the wall, and it is represented by

y+ = '-,. (2.29)
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

P (2.31)

The low-Reynolds-number k-E model requires a relatively fine grid resolution near

the wall regions. For detailed descriptions we refer to Refs. [46, 73, 87].

To make appropriate use of the two versions of turbulence models, we strive to

adjust the value of y+ of the first point next to the wall to satisfy certain requirements.

For the low-Reynolds-number approach, y+ should be less than 2 for an accurate

result, while y+ in the range between 30 and 500 [103] can produce accurate results

for the wall function treatment. An equation for the flat plate boundary l -.-r theory

is used in the grid generation process to get an estimation of proper cell size for y+,

+ = 0.172 []Re0O9. (2.32)

It also should be remember, while the position of the first cell is critical, it must have

enough cells inside the boundary l~-.-r to resolve the steep variations of flow variables.

The success of the recent turbulence closure is restricted to the sufficiently large

enough Reynolds number turbulence flows. Further refinement is needed before they

are used with confidence to calculate near-wall and low Reynolds number flows. A

more detailed review about turbulence models for near-wall and low Reynolds number

flows is given by Patel et al. [73].

2.3 Laminar-turbulent Transition

A reliable computation of the boundary 1~.-r transition remains one of the most

challenging demands in aerodynamics. Since the transitional separation bubbles and

their losses p1 i, an important role for the pressure and drag distributions, an ac-

curate representation of both laminar and turbulent separated flows is critical when

drag prediction needs to be accurately predicted. However, the fundamental fluid dy-

namics problem of transition to turbulence has resisted any simple phenomenological

description. There is still no comprehensive theory of transition.

The range of existing transition prediction methods extends from simple empiri-

cal relationships through those based on parallel and linear stability theories, like the

eN methods [15, 75], or linear or nonlinear parabolized stability equations (PSE) [31],

or the low-dimensional models [76], to direct numerical simulation, like large-eddy

simulation techniques.


For problems involving interactions between fluid flows and moving structures,

the geometry of the solid object is not known a priori. In the course of computation it

is necessary to track the geometry, the field equations, and the interfacial condition to

ensure that all requirements are met. To facilitate the solution of such moving bound-

ary problems, the moving grid technique (or dynamic grid technique) is employ, .1 to

adjust the computational grid dynamically along with the geometric updates. It is

desirable but difficult to develop an automated regridding procedure to ensure that

the new grid not only matches the geometric changes but also maintains satisfactory

characteristics such as non-ovi. iI1 pplli, smooth, and not excessively skewed. The

moving grid technique can also be used by shape optimization [54].

For moving grid problems, several alternative approaches have been proposed.

Schuster et al. [81] used an algebraic shearing process in their static aeroelastic anal-

ysis of a fighter aircraft. The basic idea is to assume that the displacement of the

subject is redistributed along the grid line, which connects the aeroelastic surface

to the outer boundary. This simple method works quite well for single block grid

with modest deflection; however, when applied to multiblock grid arrangements, this

approach usually needs extensive interventions. The spring analogy method was first

proposed for unstructured grids [3], and later was extended to structured grids by

Robinson et al. [78]. This method regards all grids as connected by springs, and each

spring has the stiffness inversely proportional to the spacing between target vertex

and its neighbors. Compared with other currently used grid regeneration methods,

the spring analogy method needs more memory and CPU time. Direct transfinite

interpolation by Eriksson [18] is a popular method because of its efficiency for struc-

tured grid. Hartwich and Agrawal [29] also proposed the Master/Slave concept to

expedite the grid regeneration process and minimize the user intervention.

3.1 Perturbation Method

For small amplitude displacement on one face of the computational block, a

simple but efficient one-dimensional perturbation method [56, 77] is can efficiently

regenerate the grid. The one-dimensional perturbation method obeys the following


Xnew X od + old (new Xld), (3.1)

where xi represents the location of an arbitrary interior grid point, x, represents the

location of a grid point on the moving boundary, and S represents the normalized arc

length along the radial mesh line measured from the outer domain, which is given by

e lil + (+1 yi2 + (1+1 Z)2
Si = (3.2)
Z7 (x xI)2 + (y1 yi + (z1+1 () 2

To use this one-dimensional perturbation method, one needs to know the displace-

ments of the two ending points of a mesh line. The positions of the interim points

are solely determined by the displacements of these ending points and their original

positions as shown in Eq.(3.1). These displacements are considered as perturbation

sources. Intuitively, the perturbation will spread throughout the whole domain while

exerting more effects on nearer points. Based on Eq.(3.1), a more complicated 3-stage

perturbation method, analogous to the transfinite interpolation (TFI) method [18],

was proposed for three-dimensional problems treated by the single-block approach.

Unlike the original TFI method, this perturbation method considers the original grid

distribution, and preserves the original grid quality.

Take a two-dimensional case as an example, in Fig. 3-1, a rectangle undergoes

shape change. To regenerate the grid for the changed geometry, the perturbation

method proceeds as a 2-stage process. The first stage is to match corner points, i.e.,

to move the corner point A to its new position A', and B to B', etc. (See Fig. 3 -b).

Once the corner points match their new positions, interim edges are generated based

on Eq. (3.1). To adjust the position of an internal point O, effects from the four

corner points need to be accounted for in a coordinated manner. In this approach,

P', Q', S', and T', which connect 0' to its boundary, will be first adjusted. The

displacements of the four edge points in, e.g., the x-direction (AP, AQ, AS, AT) can

be calculated with Eq. (3.1), and the movement of 0' is evaluated as a combination

of these displacements. Specifically, one can adopt the following expression:

A' abs(AI)AI + abs(AJ)AJ
max(abs(AI) + abs(AJ), 10-6)

where, Al = (1 Si)AP + SiAQ, AJ (1 Sj)AS + SjAT, Si and Sj are the

normalized arc lengths along mesh lines PQ and ST respectively. Displacements in

other directions can be computed in the similar way.

After the first stage, the four corner points have matched their new locations;

however, the interim edges do not necessarily match their final positions. The second

stage is to match the interim edges to their final edge locations. The difference be-

tween the interim and final edge locations are taken to be the source of perturbation.

In this stage, since the perturbations in both directions are independent, the contri-

butions from all individual directions are added [See Fig. 3-lc]. For three-dimensional

problems, a 3-stage process perturbation method was proposed in Refs. [56, 77].

3.2 Master/Slave Concepts

The aforementioned perturbation method can be implemented efficiently for sin-

gle block and small displacement problems. However, when a multi-block grid is

involved, this method needs to be enhanced by other techniques. The Master/Slave

concept acts as a useful approach to maintain grid quality and to prevent potential

grid cross-over.

I- SC S'

-P .. Q------- -- (a)
o/ Q

XA B A' B'


S c '

P I(b)

*o Q' Q

T f (
A B A' '



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.


Although unstructured grid has been used to solve fluid dynamics problems by

different researchers [3, 26], here, our attention is limited to structured grid only. For

a multi-block structured grid, CFD software often requires point-matched grid block

interfaces. A basic requirement for grid regeneration in moving boundary problems is

to maintain point-matched interfaces. While this updating process is easily satisfied

by those block surfaces that coincide with the body surface, the remaining block

surfaces need to be handled with due care. To use the perturbation method, one

needs to specify the new positions of the off-body surface points, or, at least, the

vertices of such block surfaces.

When an object changes its shape/position, master points, which are defined as

the points located at the body surface, first move, and then affect distributions of the

slave points, the off-body points. While the task of redistributing interior grids can

be (d ill. i-ii] in practice, it is more difficult to move the vertices of each block if a

point-to-point matched interface without overlap is required. If two abutting blocks

have congruent interfaces, such as the interface of block 2 and block 3 in Fig. 3-2a,

once the corner vertices are determined, the edge points and the interior points can

be uniquely settled based on Eq. (3.1). However, when two neighboring grid blocks

do not share identical interfaces, such as block 1 and block 2 shown in Fig. 3-2a, this

procedure can cause discontinuity at the interface because the regriding procedures

at the interface are based on different stencils for different blocks. Slaving vertices to

their master points can avoid creating undesirable grid discontinuities.

The Master/Slave coupling is based on the distance between an off-body point

(Slave point) and a surface point (\! ,Iter point). The distance is given by:

|r| = [(X, Xb)2 + (Y, b) + (, )2]0.5, (3.4)

where subscript v represents an off-body vertex, and, b, a body surface point. For

each off-body vertex, the nearest body surface point is defined as its master point.

(a) (b)


BloBlock ck 3

Block 1 (2,0.9,1)
(0,0 B lock 2


(c) (d)



Figure 3-2: Tests on the moving grid scheme: (a) initial geometry; (b) side view of
the grid; (c) deformed shape; (d) new grid distribution; (e) grid distribution with
high stiffness parameter.

To simplify the connection between the master and slave points, one can convert the

three-dimensional parameter space (i, j, k) into a one-dimensional data structure.

Each grid point has an identification number 1 defined by

I = boff(bn) + i + imax X (j 1) + imax X max X (k 1), (3.5)

where boff(bn) is the identification number of the first point in block bn in the one-

dimensional array, imax and jmax are the maximal dimensions in the i and j directions,

respectively. In this way a slave point is associated with its master point by storing

its master's identification number into its address.

Once the relation is established, the next step is to determine how the slave

point moves according to the influence from its master point. Intuitively, a nearer

vertex will have more effect (displacement) than those far away. Also, to avoid a

lower surface to cross over its opposite face when the movement is relatively large,

one needs to scale the movement. A simple but effective formula is expressed as


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.3 Numerical Tests

The moving grid technique consists of two elements, one is the perturbation

method, and the other is the master/slave concept. The proposed method is tested

on a three-dimensional, multi-block, and structured computational grid. The initial

configuration is shown in Fig. 3-2a, block 1 shares the interface with blocks 2 and 3.

The initial grid is shown in Fig. 3-2b, it has clustered distributions at the top and

bottom surfaces, and it also has point-to-point matched interfaces. The bottom two

faces are the moving boundaries that experience large deformation. The coordinates

of each vertex are labelled. Comparing Fig. 3-2a and Fig. 3-2c we can see the

displacement is as large as half of the total height. The new computational grid

is computed with the method discussed in this chapter. The new grid is shown in

Fig. 3-2d. That figure does not show any cross-over in the new grid after the large

deformation. The new grid not only keeps the point-to-point matched interfaces, but

maintains the clustered distributions at the top and bottom surfaces. The CPU time

for regenerating the computational grid is less than 10' of the CPU time used by

the fluid solver.

The parameter 3 controls the stiffness, we find /3=1/64 works well for our prob-

lems. A larger 3 causes the block to behave more like a rigid body (See Fig. 3-2e).

More details about the master and slave concept can be found in Refs. [29, 56].


4.1 Introduction

The main objectives for fluid dynamists and aeronautical engineers throughout

the history of modern flight have ahlv-b-, been to increase lift and reduce drag. Micro

air vehicles are airborne vehicles with characteristic lengths under 6 inches traveling

at speeds around 30 ft/s. They are designed to perform an array of different tasks

such as reconnaissance, target designation or border surveillance. To enhance the

performance of these small vehicles, improvements should be made in the respects

of power plants, thrust generators, navigation systems, cameras, data transmitting,

and more importantly, their aerodynamic performance. As far as the aerodynamics

is concerned, MAVs have to operate in one of the worst flight regimes imaginable, one

of them is the low-Reynolds numbers which is usually norotiously accompanied with

laminar boundary separation, transition, and low lift-to-drag ratio. Furthermore the

small dimensions of the MAVs greatly enhances the trade off difficulties mentioned

above. The small size inhibits the designer's freedom to achieve the desired aerody-

namic performance. The limitations in the power sources available, at present time,

make this equation even more complicated.

To solve the aerodynamic challenges of these problems, two in ii r fields have

been evaluated. The first approach is trying to copy nature's own way; that of flapping

wings [43, 45, 59, 83]. The second approach is using a dynamical wing There are

some different implementations in this field. One widely studied method is the use of

an wing with a flexible membrane upper surface [53, 55, 83, 92]. Another method uses

active devices to control the flow on the upper surface of the wing. The latter method

has gained more attention with the success microelectromechanical systems ( \!i \!S).

In active flow control area, a piezo-actuated flap mounted on a low Reynolds number

wing have been studied both experimentally [19] and computationally [30, 56]. Gad-

el-Hak [20] has discussed the possibility of employing MEMS to delay or prevent

separation on MAVs.

To understand the low Reynolds number wing characteristics and improve the

wing performance, we will perform a computational investigation for flows surround-

ing a dynamically shaped wing, at a chord Reynolds number of 78,800, along with

a parallel experimental effort. We adopt a piezo-actuated flap on the upper surface

of a fixed wing for active control. The actuation frequency they focus on is 500 Hz.

Their computational framework consists of a multi-block, moving grid technique to

handle the geometric variations in time, using two-equation turbulence closures, and

a pressure-based flow solver. The dynamic grid redistribution employs the transfinite

interpolation scheme with a spring network approach. Comparing the experimental

and computational results for pressure and velocity fields, effects of the detailed ge-

ometric configuration of the flap, the flapping amplitude, turbulence modeling, and

grid distributions on the flow structure will be assessed in this chapter.

4.2 Experimental and Computational Set-up

The wing has moving flaps on the upper surface as shown in Fig. 4-1. The flap

can be actuated at prescribed frequencies [19]. The wing studied in the experiment has

a span of 1 ft and a chord of 5.4". In the experimental set-up, there is a gap of 0.016"

wide between two flaps, each is 2" wide and 0.564" long. The Reynolds number, based

on the freestream velocity U=27.5 ft/s, is 78,800. The flow is initiated as laminar;

after the transition point, the turbulence model is invoked. The transition point is

identified to be at the tip of the flap. In the computation, the turbulent Reynolds

number based on the chord length, Ret=pUL/pt, at the inlet of the turbulent flow

region is set to 500. The amplitude and shape of the flap are prescribed by analyzing

the experimental data obtained with PIV [19] For the 500 Hz case, two flapping

amplitudes are considered, namely 0.001" and 0.0073".

region of interest piezo actuator (flaps)
for PIV \ /
\o / pressure taps, 0.25 in. apart



Figure 4-1: Moving flap wing experimental setup in the wind-tunnel

In this chapter, a number of two-dimensional and three-dimensional cases, in-

cluding different grid densities, with and without the gap, will be studied. For 3-D

cases, an effort has been made to consider the effect of the gap between flaps as shown

in Fig. 4-2. The original Launder-Spalding k-E model [52] with the wall function, as

well as a low-Reynolds number version k-E model [46] are compared.

Linear Chords -

Figure 4-2: Flap wing with refined three dimensional bending model.

4.3 Numerical Study

4.3.1 Grid Sensitivity Analysis

A number of 2-D grids with different distributions and sizes have been tested.

Then a reference 2-D grid is selected by ensuring that satisfactory solution quality

is obtained on a most economic grid size. Then, the 3-D grid is established by

extending the reference 2-D grid along the spanwise direction. The resulting grid,

for cases with gap, has 33 nodes in the spanwise direction, and consists of 20 blocks

with approximately 533,000 grid points. For cases without gap, the grid system has

20 nodes in the spanwise direction, and consists of 8 blocks with about 253,000 grid

points. With the multi-block strategy it is convenient to employ the transition model,

i.e., when a transition location is identified, new computational blocks are generated

surrounding the transition point, so that the laminar flow equations are solved in

the blocks before the transition location, and the Reynolds-averaged turbulent flow

equations are solved in the blocks after the transition point. In our computations,

the transition point is identified from the experiments [19], and it is at the tip of the

flap, therefore, our computational blocks are set up based on this transition point.

4.3.2 Assessment of the Turbulence Model

In the early work by He et al. [30], the wall function approach is used exclusively.

For such a low Reynolds number, the low-Reynolds number version of the k-e model

seems more suitable [56]. Fig. 4-3 compares the time-averaged pressure distributions

at 500 Hz between two different turbulence models for the case without gap between

flaps. Both experimental and computational data are collected at the middle of the

flap and along the chordwise direction. As far as the computational time is concerned,

it is observed that the low-Reynolds number model uses less CPU time than that of

the wall function treatment.

4.3.3 Effects of Gaps between Flaps

The inclusion of a gap obviously increases the complexity of the flow field in the

spanwise direction. It can be seen from Fig. 4-4 that streamlines enter the gap and

move toward the middle of the flap. A weak recirculatary flow with a vortex stretching

from the wing surface to the flap is present on both sides of the gap, rotating in the

clockwise direction on the right and counter-clockwise direction on the left of the gap

when viewed from above.

The gap also affects the overall pressure distribution on the wing. Fig. 4-5

compares the pressure distributions obtained with and without gap between flaps.

The low-Reynolds number turbulence model is used in both cases. Time average is

- Low Re Model
-- Wall function model



0.3 -



1.5 2 2.5 3 3.5 4

Figure 4 3: Comparison between tine-averaged pressure distributions for cases
without gap, at 500 11z, obtained with different turbulence models. (: : :
for the entire cycle, .i.' ude: I01").

-\ '

Figure 4 4: Streamline for 3-D case with gap at 500 Hz, amplitude is 0.'-i1"

made for the entire cycle for an amplitude of 0.001" and an actuation frequency of

500 Hz. Since the pressure taps are located at the middle of a flap, the effect of

the gap may not be clearly demonstrated by the figures shown. Fig. 4-6 compares

the velocity profiles immediately downstream from the flap trailing edge. The low-

Reynolds number turbulence model is used again and the comparison is made at

0.0162" downstream from the trailing edge of the flap. Cases involving both stationary

flaps and actuated flaps, at 500 Hz, are examined. Consistent with that already

observed, the velocity profiles indicate that the flow enter into the gap from the

upper surface of the flap, maintaining the same streamwise direction while moving

laterally toward the region between the flap and the wing surface. A recirculating

flow pattern is formed between the flap and the wing surface.

No gap between flap
0.65 Experimental
.--- Gap between flap



O 0.45




0.25 -

1 1.5 2 2.5 3 3.5 4

Figure 4-5: Comparison of time-averaged pressure distributions for cases with and
without gap between flaps.

4.3.4 Effects of Flapping Amplitude

Depending on the flapping amplitude, the reattachment point can move back

and forth, and the recirculating flow can either remain attached to the tip of the flap,

or be shed away. In our computation the curvature of the flap is set by collecting

10 5 10' 105
10- 10-
95 95 /
9 9-
S85 / At the gap 85
S- Middle of the flap 0

- 75- 75-
7 /7
655- 65
6 5 6 At the gap
Middle of the flap
55 55
2 0 2 4 6 8 10 12 -5 0 5 10 15
U velocity foot/s U velocity foot/s

Figure 4-6: Comparison of velocity profiles at two different spanwise positions with
gap implemented: (a) steady computation; (b) f=500 Hz.

the experimental data and then approximated by a simplified model. For the 500 Hz

case in our computation two amplitudes are considered, one is 0.001" and the other

is 0.0073". For the smaller amplitude case, the motion of the trailing edge is modeled

as the following,

y(z, t) = (8.6 x 10-4 + 1.6 x 10-4 sin( )) sin(27f t), (4.1)

where z is the spanwise distance from the study point to the gap, f is frequency, and

t is time. The larger amplitude case simply scales the above equation based on the

same expression.

Fig. 4-7 compares the time-averaged pressure coefficients between the two flap-

ping amplitudes. For the higher flapping amplitude, the time-averaged reattachment

point moves closer to the trailing edge of the flap, because a larger movement of the

flap introduces more high momentum fluids into the boundary 1 r. In the experi-

ment we observe vortex shedding as shown in Fig. 4-8a. In our computations with

the smaller flapping amplitude, the vortex core is anchored at the flap tip, exhibiting

two co-rotating cores as the flaps move upwards, as shown in Fig. 4-8b. On the other

hand, with the higher amplitude, vortex shedding appears for cases with or without

0 Disp=l.0x 10-3 in
SDisp=7.3 x 10-3 in


0 0.35-



1 1.5 2 2.5 3 3.5
X position(in)

Figure 4-7: Pressure coefficient at different actuation amplitude without gap between
flaps. Actuation frequency: 500Hz.

the gap. Fig. 4-9 shows selected instantaneous plots of the vortex structure during

a single cycle of the flap movement. As expected, the recirculation is stronger when

the flap moves upward and becomes weaker when the flap moves downward.

Figure 4-8: Instantaneous trailing edge streamlines for the frequency of 500 Hz as
flap is moving upwards for amplitude of 0.001". Left: experimental result; right:
computational result.

4.4 Conclusions

Base on our results we find that the low-Reynolds number formulation of the k-e

turbulence model produces longer circulation zone. The gap between the flaps affects

the location of the re-attachment point, and creates a three-dimensional flow pattern.

Aided by the limited experimental information, the computational model can offer

Trailing edge of iap

(a) Phase angle a = -/6.

(c) Phase angle a 7i /6

(b) Phase angle a = r3.

(d) Phase angle ao 8-r/6.

(e) Phase angle a = 11T/6.

Figure 4-9: Trailing edge streamlines
with an ..:.. i'tudc of 0.C" '

(f) Phase angle a = 2-.

. the f- : ..-- of 500 Hz in the third cycle

insight into the fluid physics. Vortex dynamics changes substantially, in accordance

with the flapping amplitude. Vortex shedding occurs when the amplitude is higher,

while anchored vortex is observed with the lower amplitude.

As to the physical mechanisms, while it is expected that the vortex shedding

appears when the amplitude is high, it is not clear why and when this phenomenon

ceases to exist when the amplitude is reduced. Similarly, while it is obvious that

the gap between flaps affects the flow reattachment, it is hard to predict the trend



5.1 Introduction

The deformations of membrane structures are often large and thus inherently

nonlinear. Qualitative descriptions of the thin membrane behaviors are often more

complex than that of three-dimensional bodies. Until recently, membrane analyses

were primarily relied on trial and error. Green and Adkins [25] are the first to

derive a general non-linear theory for membranes. The membrane theory has two

fundamental assumptions because of their thinness [41]. First, the membrane tension

stiffness is much larger than its bending stiffness, therefore, the stress couple effects

can often be neglected. Second, the ratio of the membrane thickness to its radius

of the curvature is small, which effectively decouples the strain-displacement relation

from the curvature tensor. Practically speaking, there exist two approaches to model

the response of a membrane structure. The first approach is to idealize a plate or

shell structure as a membrane away from its boundaries, in which the stress couples

are neglected in the interior region. By formally equating the plate stiffness to zero,

this membrane model can be described by the Foppl-von Kdrmdn theory [40, 41].

In the second approach, the structure cannot sustain stress couples anywhere. The

second approach is adopted in the current work because it is appropriate for MAV


Two-dimensional membrane wing study with fixed leading and trailing edges

were originally conducted by N. 1- [68], Thwaites [101], Voelz [108]. These works

considered inextensible membranes surrounding by steady, two-dimensional, and ir-

rotational flows. As a consequence of the inextensible assumption, when cambers and

incidence angles are small, problems can be linearized and expressed compactly in a

integral form which is referred as the "sail equation" by N. i.--in ,, [69]. The following

equation completely defines the linearized theory of inextensible membrane wings in

steady and inviscid flow fields.

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


The development of a three-dimensional membrane model introduces several

complicating factors not encountered in two-dimensional analyses. First, unlike the

- il equation", the tension in a three-dimensional membrane is no longer a single

constant value but is at best described by a state of biaxial tension along the lines

of principal stresses [38]. Second, geometric and material properties vary along the

spanwise direction, and they have to be described accurately. Third, if one of the

principal tensions vanishes, the membrane will be compressed and wrinkled, which

further complicates the analysis. Considering the complexity of the three-dimensional

problem, simplifications are made in the early works.

Oden and Sato [70] presented a finite element method for the analysis of large

displacements and finite strains in elastic membranes of general shape. In their work

the static equilibrium of an inflated membrane was considered. A static analysis can

be considered as a special case of a dynamic analysis in which the inertia effect is

included. Works on the dynamic analysis of membranes before 1991 were reviewed by

Jenkins and Leonard [41]; most recent works were reviewed by Jenkins [40]. Verron

et al. [105] performed a combined numerical and experimental study of the dynamic

inflation of a rubber-like membrane. Membrane wrinkle problems were tackled in

the work of Ding et al. [14]; they performed a numerical study of partial wrinkled

membranes. The stability issues of a fluid flow past a membrane were reported by

Thaokar and Kumaran [99].

In this chapter a three-dimensional membrane model is proposed for the mem-

brane dynamics. The membrane material considered obeys the hyperelastic Mooney-

Rivlin model [64]. We use the Green-Lagrange strain tensor for the description of

large strains. The dynamic response of the membrane is described by a system of

second-order ordinary differential equations. A finite element method with triangu-

lar elements for spatial discretization is utilized. For time integration, the implicit

Wilson-O method is used.

5.2 Structural Solver

5.2.1 Mooney-Rivlin Model

For an initially isotropic membrane, Green and Adkins [25] showed that there

existed a strain energy function W which was expressed as the following form

W =W(II, 12,1), (5.2)

where 1i, 12, and 13 are the first, second, and third invariants of the Green deformation

tensor C. If the material is incompressible, namely, 13 1, then the strain energy is a

function of IA and I2 only. It can be written as:

W -= ci( 3) + c2(2 3), (5.3)

where ci and c2 are two material parameters. A material that obeys Eq. (5.3) is known

as the Mooney-Rivlin material [64], and the model is called the Mooney-Rivlin model.

The Mooney-Rivlin model is one of the most frequently employ, ,1 hyperelastic models

because of its mathematical simplicity and relatively good accuracy for reasonably

large strains (less than 150 percent) [64].

For an initially isotropic membrane, the general stress-strain relationship is writ-

ten as
S = -pC- + 2 (5.4)

where S is the second Piola-Kirchhoff stress tensor, p is the hydrostatic pressure. If the

membrane is incompressible, then the stress-strain relation can be further simplified

as the following

S -pC-1 + 2 [(cl + c21)I c2C] (5.5)

where I is 3-by-3 identity matrix. In the membrane theory the stress field is essentially

treated as two-dimensional, and the stress normal to the deformed membrane surface

is negligible in comparison with the stresses in the tangent plane [70]. Under this

consideration, the deformation and stress matrices, C(t) and S(t), are given by

CiI(t) C12(t) 0 S11(t) S12(t) 0
C(t) c 12() 22 () 0 S(t) 12(t) S22(t) 0 (5.6)

0 0 C33(t) 0 0 0

Therefore, the hydrostatic pressure is determined by considering the condition that

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


f DT (t)poD(t)dV + E: SdV- f o DT(t)P(t)dS = 0, VOU, (5.9)

where Vo and iOVo are, respectively, the volume and the boundary surface of the

undeformed membrane, po is the membrane density, D(t) is the displacement vector

in global coordinates, 6D(t) is the virtual of vector D, the superscript T refers to

the transverse of a vector or a tensor, D(t) is the acceleration vector, P(t) is the

generalized external surface force, E is the Green-Lagrange strain tensor, S is the

second Piola-Kirchhoff stress tensor, and 6E:S is the scalar product of two tensors.

In the rectangular Cartesian coordinates, the scalar product is defined by

E : S =6EijSji, (5.10)

with summation on the repeated indices.

5.2.3 Triangular Membrane Element

The principle of virtual work can be written in the following discretized form

61I (f DTpoUdV + 6E: SdV / DTPdS =0, V6D, (5.11)
i \JVo0 V0 ,,

where 61 is the total virtual work, and the subscript e represents elementary quan-

tities, N, is the number of elements covering the entire membrane. In the present

work the triangular element is employ, As shown in Fig. 5-1, a typical triangular

element, e, is defined by nodes, i, j, k, and three straight-line boundaries. In the

global coordinates the three nodes have coordinates (Xi, Y1, Z.), (Xj, Y, Zj), and

(Xk, Yk, Zk), respectively.
To facilitate the development of a computational procedure, it is convenient to

consider the membrane element and its displacement in the local coordinates. As

shown in Fig. 5-1, in local coordinates, the triangular element lies in the plane of

k' y=x+u





Figure 5-1: Schematic global and local coordinates.

;, Nodes i, j, and k, have local coordinates of (xi, yi, 0), (xj, yj, 0), and (xk, yk,

0), respectively. For each point within the triangular element, its local coordinates x

and global coordinates X have a relation expressed as follows:

X X + Tx, (5.12)

where Xi is the global coordinate of node i, and T is an orthogonal transformation

matrix which satisfies TT= T-1. Similarly the displacements have the following


D =Td, (5.13)

where D is the value of a displacement in global coordinates, d is its value in local


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,


where N3 is an interpolation matrix defined as

1- -I 0 0 0 0 r 0 0
N3(,r)= 0 1 0 0 00 0 (5.5)

0 0 1- 0 0 0 0 rI

where 0 < ~< 1, 0

de di dj dk (5.16)

in which di, dj, and dk are the nodal displacement vectors of nodes i, j, and k,

respectively. Each node has three degrees of freedom,

di= u (5.17)

5.2.4 Green-Lagrange Strain Tensor

When a displacement is introduced the location of a point originally given by

coordinates x now has new coordinates

y x + u, (5.18)

where u is the displacement as shown in Fig. 5-1. Following Green and Adkins [25],

the Green-Lagrange strain tensor is defined by the relation

1 ( Oyk Yk (5.9)
E =i -xi -xj 5i (5.19)

Since the thickness of the membrane is small compared to the other characteristic

dimensions, the following simplifications are allowed to be made [25, 70]

_1 (Qu 8 9uP 9ui 9ui\
E t ( -, + U-- + N O
2 x, 3 dx dOx, dxo
E3 = 0, (5.20)
E33 2 (A 1), C, = 1, 2,

where A is the same as A3 in Eq. (5.8). To simplify the computation, we divide the

Green-Lagrange strain tensor into two parts. One part is Eo, which takes into account

small strains,

u 1 iu av
+ )O
ax 2 O y ax
1 au av av
Eo (= t + )v 0
2 ay Ox ay
0 0 0

and the other is EL, which is the nonlinear part of the strain tensor,

ou u Ovu v Owv w Ouw u Ovu v Ov w Ow
Ox Oxox ox x Oxx Ox Oy Ox 9y Ox Oy
SOuu + Ovv v Oww Ouw u u Ovv + Oww 0w
EL = + + + +0
Ox ay Ox y Ox y OyOy ay y ay y y
0 0 A2



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


F'nt j BT((d))aT(e)dV, (5.24)
where a is the nonlinear stress vector, e is the strain tensor, and B is the strain-

displacement matrix which is also a function of the displacement for large deforma-

tion. Recall that in the principle of virtual work the term related to the internal

forces from the discretized form of principle of virtual work Eq. (5.11) reads

S 6E : SdV. (5.25)

Eq. (5.6) tells us that the membrane stress is essentially two-dimensional, therefore,

we know from Eq. (5.10) that it necessary to consider the reduced strain components

instead of all the components when calculating the internal forces based on Eq. (5.24).

The general membrane strain vector can be conveniently decomposed into two parts

according to the strain decomposition in Eqs. (5.21) and (5.22),

E = o + L, (5.26)

where so is the reduced strain vector for the linear part of the strain tensor, Eo, and

EL for the nonlinear part EL. Based on Eqs.(5.21) and (5.22) these two strain vectors
can be expressed as the following,

S u 'u Ou Ovv Ov aw w

Once the virtual displacements are given at any point within the element, the strain

variation at any point can be determined. The relationships between the virtual
oxOx Ox Ox Ox Ox OX OX

displacement and strain variation arev v

< 60 = Fgy B0de, 6L { j Ly } BLd, (5.28)
9u +9Ov 9i2-u u +9v v +9Ow Ow

6x,,, 5ELay
f< y Ox Ox 9y Ox 9y Ox 9y ,

where and B are the stradisplacements are given at any point wicesthin the expressions can be straifound
variation Appendix A. Based on Eq. (5.27) the virtual work in each triangular element canelationships between the virtual
displacement and strain variation are

6Eo 6FOy Bo6de, 6L 6ey BMde, (5.28)

6-11, 6FLxy

where Bo and BL are the strain-displacement matrices, their expressions can be found
in Appendix A. Based on Eq. (5.27) the virtual work in each triangular element can

be computed as

/ 6E : SdV = 6DT3 f BTa dV = D TFnt, (5.29)

where B = Bo + BL, 6eT = Le + bg, and Fnt is the internal force in global


Fi't = TS BTa dV. (5.30)

T3 is the a 9-by-9 transformation matrix defined as

T 0 0
T = 0 T 0 (5.31)

S0 T

in which T is the same transformation matrix as shown in Eq. (5.13).

5.2.6 External Forces

As the membrane changes its shape and position under the pressure, not only

the direction of the pressure changes but also the area on which it acts changes

[70, 105]. To account for this situation, we assume that the size of each finite element

is sufficiently small that the pressure is uniform over each element. To further simplify

the computation, we again use local coordinates. As shown in Fig. 5-2, the local

coordinates lies in the deformed plane instead of in the undeformed plane. In the

local coordinates the pressure P is given by,

P = Pn, (5.32)

with n [ 0 0 1 ]T is the external normal direction of the element.

The virtual work by the external force at each element is


k X

1 z

/ Y


Figure 5-2: Schematic of local coordinates for external forces.

the net external force at each node is obtained by simply representing it by three

equal forces at each node

SDTPdS = DT3oP 0 0 0 0 3 0 0 1 6D Ft. (5.33)
J n' ,, L3 3 3

Hence, in global coordinates the node force on element e due to a uniform pressure

over that element is

F T3OP o0 03 O 0 1 0 0 1 (5.34)

5.2.7 Dynamic Equations for the Membrane

The first term of the left side of the principle of virtual work equation (5.11) can

be written as

VI 6DTpoDdVo D= DMeD, (5.35)

where Me is the mass matrix of an element,

21 I I
Me = posoho I 21 I (5.36)
I I 21

in which I is a 3-by-3 identity matrix.

Up to now we have obtained the expressions for the internal force, external force,

and mass matrix on an element e. It is necessary to assemble all the discrete elements

into a single unit. For any element e the principle of virtual work is

6H1e = JDTpoDdV + 6E: SdV 6DTPdS. (5.37)
Vo IVo ,,

Assembling all the forces in one element we have the following equation,

Se = 6D> (MD + F'nt Fext (5.38)

This equation is for a single element, if we sum all elements in global coordinates,

then we obtain the system of governing equations for membrane under external loads,

MD(t) + F'nt Fext, (5.39)

where M is a positive definite mass matrix, which remains constant, Fint is the

internal force, and Fext is the external load.

5.2.8 Solution of the Dynamic Equations

Either an implicit or an explicit scheme can be adopted to follow the time history

of the system of equations. The explicit scheme is computationally economic for each

time step and requires less storage. Its efficiency can be further improved by special

techniques, such as lumping technique of the mass matrix. However, its time step

needs to be small to satisfy the stability requirement. As far as the implicit scheme

is concerned, the matrix system is updated more than once per time step to account

for the nonlinear effect, which requires more computations per time step and more

storage. Generally speaking, implicit algorithms are effective for structural dynamics

problems in which the response is controlled by a small number of low frequency

modes, while explicit algorithms are efficient for wave propagation problems in which

the contribution of intermediate and high frequency structural modes to the response

is important [94].

There are numerous v--v to integrate the structural dynamics equations in time;

among these we utilize the Wilson-O method [113]. The Wilson-O method has a

controllable algorithmic numerical dissipation, it is a one-step method, and it yields

second order accuracy in time. These characteristics make it easier to code and

behave well in our fluid/structure interaction system. Since the Wilson-O method is

dissipative, we employ a time step smaller than that required by the overall accuracy

criterion to avoid excessive dissipation. This approach, with 0=1.4, works well for

the present problems.

5.3 Numerical Test of Membrane Model

The inflation of a spherical membrane under an imposed pressure step was an-

alyzed experimentally by Alexander [1] and theoretically by Beaty [4]. Recently, a

numerical study was performed by Verron et al. [105]. The corresponding semi-

analytical solution was reported by Verron et al. [104]. These studies focus on an

incompressible spherical hyperelastic membrane. The membrane has an initial radius

of Ro, an undeformed thickness of ho, and a mass density of po, respectively. The

membrane material is assumed to obey the Mooney-Rivlin model. There are two con-

stitutive parameters in this model, namely cl and c2. To simplify our computation,

we define a non-dimensional Mooney-Rivlin constant based on these two constitutive


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)

where A is the circumferential principle extension, which is defined as the ratio of the

deformed to the undeformed radius. The corresponding non-dimensional time r, and

the non-dimensional imposed pressure step p* are given as follows:

t 2 c, (5.42)
Ro V Po

P* = o* 4 (5.43)
where Po is the dimensional pressure step. The initial conditions of the membrane

may vary; here, we focus on the case that the membrane is initially in its equilibrium

state, i.e., with the following initial conditions:

A(r 0) = 1 and A(r 0) = 0. (5.44)

One can obtain the semi-analytical solution to Eqs. (5.40) and (5.43) using any stan-

dard solver. The direct numerical solutions of the membrane equation are obtained

with the finite element method discussed before. The computational mesh has 816

nodes and 1568 triangular elements on a semi-spherical surface. In our work we

choose At5 x 10-4 which is 500 times larger than the explicit scheme used by Verron

et al. [105] where they use an explicit scheme with a time step of 10-6 for the same


All computations are performed using the double-precision. In Fig. 5-3 we com-

pare the finite element results with the semi-analytical solutions. Except for one case

in Fig. 5-3b, all numerical solutions match well with the semi-analytical solutions.

The loss of linearity of oscillations is well reproduced. The mismatched case shown in

Fig. 5-3b corresponds to a=0.1. For this value of the material parameter, as pointed

out by Verron et al [105], there is an unstable equilibrium point at p*W0.687. Near the

unstable equilibrium point, the finite element results can only approach qualitatively

the real behavior. This is the reason why the numerical results are highly sensitive

to the assigned pressure step. As expected, the larger the Mooney-Rivlin constant a

is, the higher pressure the membrane can sustain.

Two profound issues should be considered in the use of the Wilson-O method.

First, the Wilson-O method is unconditional stable for linear problems. The time step

in nonlinear problems is often dictated by the errors introduced in the approximation

treatment of the nonlinearities; therefore, due care should be taken to choose an

approximate time step. Second, a small time step is desirable for the Wilson-O method

due to its excessive dissipation. As shown in Fig. 5-4, a larger time step tends to

cause more numerical damping. In Fig. 5-5 we show the error norm of the implicit

scheme, the figure exhibits a slope around 1. Therefore, the overall accuracy for time

integration is approximately first order.

The present membrane model gives a good agreement with the analytical solu-

tions for membrane dynamics under large deformation. It leaves one issue of funda-

mental importance since it cannot handle wrinkle phenomenon which happens when

the membrane is compressed. The membrane wing is sometimes compressed, espe-

cially at the leading edge. Therefore, further investigation is necessary in this area.

p =0.6

I0 5
0 5

15 20

p=1 2


10 15 20

Figure 5-3: Dynamic inflation of a Mooney spherical membrane: (a) a = 0; (b) a
0.1; (c) a = 0.25. Solid line: Semi-analytical solution; Circle: Numerical solution.

A t=2.0e-4




< 1.25


10 15 20 25 0 5

10 15 20 25

Figure 5-4: The effect of time step to the numerical damping for p*=0.5 and a=0.


t-4.9 /

-5.1 /


-3.6 -3.4
Time Step (login)


Figure 5-5: Error norm versus time step for the implicit scheme for p*=0.5 and a=0.


z 1.25


0 5

A t=5.0e-4


6.1 Introduction

Membrane structures are found in many engineering practices. Traditionally,

they are utilized in parachutes and sail boats to improve performance. A recent

interesting application is a membrane wing-based Micro Air Vehicle (\!A V). MAVs

with a maximum dimension of 15 cm or less and flight speed around 10 m/s provide

expendable platforms for surveillance and data collection in situations where larger

vehicles are not suitable or else dangerous to reach. However, MAVs are sensitive

to the environmental disturbances, and it is hard to avoid considering their small

dimensions, light weights, and low flight speeds. Our experience indicates that an

aeroelastic membrane wing can better adapt to the atmospheric disturbances, provide

smoother viewing platform, and make the vehicle easier to fly [83].

During flight, the surrounding flow exerts aerodynamic forces on the membrane

wing and causes structural deformations of the membrane surface; meanwhile, the

wing deformations affect the flow pattern and change the aerodynamic force distri-

butions. The mutual interactions between the fluid flow and the structure cause a

fluid and structure interaction problem. Fluid and structure interaction problems

have attracted interest of many researches. Numerous studies have been conducted

on these subjects. Numerical schemes for fluid and structure interaction problems

may vary. However, in general, there are two key elements in a coupled fluid and

structure interaction system. One is the fluid solver that evaluates the aerodynamic

forces, and the other is the structural solver that computes the shape deformation.

Beside these, a moving grid technique is also needed to regenerate the computational

grid; an interface algorithm is utilized to exchange information between the fluid and

structural solvers.

The interaction between a fluid and a membrane usually causes substantial struc-

tural deformation that is inherently nonlinear. The membrane behavior is in general

complex and earlier studies of interactions between fluid flows and membrane struc-

tures are largely based on empirical observations [34], two-dimensional flow based

computations [36, 86, 92, 93], or simplified three-dimensional analyses [37, 38]. In

our study we present a coupled fluid and structure interaction system to investigate

the interplay between a membrane and incompressible viscous flow. In the coupled

system, the fluid dynamics is governed by the three-dimensional Navier-Stokes equa-

tions for incompressible flows, and the structural dynamics is governed by a nonlinear

dynamic membrane model. The Navier-Stokes equations are solved using a pressure-

based method on a structured and multi-block grid. The membrane material follows

the hyperelastic Mooney-Rivlin model [64]. The membrane model is derived from the

principle of virtual work. It is spatially discretized with a triangular finite element

method and temporally integrated with the Wilson-O method [113]. To accommodate

the membrane shape change, the flow solver is further augmented with a moving grid

technique discussed in C'! lpter 3.

The interaction system is coupled through the exchange of the aerodynamic

forces and shape deformations between the fluid and structure solvers. Since the

fluid solver and structural solvers do not share the identical grid at the interface, the

thin plate spline (TPS) interpolation [16] is used to map the external load and the

shape deformation to the structural and fluid solvers, respectively. To avoid phase lag

between the fluid and structural solvers, subiterations are performed to synchronize

the fluid and structural solvers.

In the following we first introduce the interface algorithm between the fluid

and structural solvers, this is followed by a step-by-step description of the fluid and

structure interaction algorithm. This coupling is achieved through exchanging infor-

mation and synchronization between the fluid and structural solvers. In the aerody-

namic analysis part, we study both the rigid and membrane wings. In the rigid wing

study, we focus on the characteristics of the low Reynolds number and low aspect

ratio rigid wing, including boundary 1 -v-r separation and reattachment, tip vortices,

and unsteady phenomenon at high angle of attack. In the membrane wing study,

we study the self-initiated membrane vibration, its impacts on the membrane wing

performance, and we also compare the membrane and rigid wing performances.

6.2 Coupling between the Fluid and Structural Solvers

When a fluid flow passes a flexible structure, the structure changes its shape

due to the external forces situated in flow; meanwhile, the shape change affects the

fluid flow around the structure and the force distributions on the structure. This

leads to a fluid and structure interaction problem. Usually, the fluid problem and

the structural problem have different mathematical and numerical properties. De-

pending on the problem complexity (linear problem or nonlinear problem, simple or

complex geometry, small-scale or large scale problem), available software, and com-

puting platform, numerous techniques have been developed to provide coupled CFD

and CSD methods. One approach is to treat the flow and structural equations as one

monolithic system of equations [5] and solve the system of equations simultaneously

on a "monolithic" grid. However, since multiple time scales present in a fluid and

structural system, a unified treatment may not be efficient to handle the computa-

tional stiffness. Liu et al. [58] solved the flow and structural fields simultaneously by

a fully implicit method. However, unlike the monolithic method, their CFD solver

and CSD solver are solved on different grid systems in which The fluid and structural

governing equations are regarded as one single system of time-dependent equations

in a pseudo-time.

To better handle different characteristics in the flow and structural domains, most

of the researchers solved the fluid and structural equations separately and then cou-

pled them through a synchronization process. Smith and Shyy [93], Kamakoti et al.

[49], Gordnier and Fithen [22], Kalro and Tezduyar [47], Zwaan and Prananta [115],

and Lian et al. [53, 55] coupled the fluid and structural solvers with a subiteration

approach. The fluid solver and the structural solver function independently on their

own computational grids with their own time steps, subiterations between these two

solvers are employ. .1 at each time step to avoid the phase lag. Alternatively, Farhat

et al. [26, 28] formulated the fluid and structure interaction problem as a three-field

problem: the fluid, the structure, and the dynamic mesh that is represented by a

pseudo-structural system. In their work, the arbitrary Lagrangian Eulerian (ALE)

finite volume scheme is emploi-y, for the fluid dynamics equations, a finite element

method is applied for the structural dynamics equations. Partitioned procedures and

-I I,.-:- red algorithms are then adopted for the coupled fluid and structure interaction

problems in the time domain.

In fluid-structure interaction computations, a moving grid technique can be fruit-

fully employ. 1 to automatically regenerate the new CFD grid according for the de-

formed configuration. In addition, since the CFD solver and CSD solver do not

necessarily share identical grid at the interface, interpolations are needed to allow

the information exchange. Smith et al. [91] evaluated suitable methods to transfer

information between CFD and CSD grids. In our study a thin plate interpolation

method [16] is adopted for this purpose. The thin plate spline interpolation is a global

interpolation method, it is invariant with respect to rotations and translations, and

hence is suitable for the interpolation on moving surfaces. A one-dimensional version

of the interpolation is

H(x) = Y aix X2log Ix Xi, (6.1)
i= 1

where H(x) is the displacement distribution function over the domain, ai is the unde-

termined coefficient, N is the number of monitored structural nodes on the interface,

and xi's are their locations. It can be seen from Eq. (6.1) that the interpolated

function H(x) has a continuous first order derivative. Once these coefficients ai are

determined based on the structural grid locations, the displacement vector defined

on the CSD grid D is mapped to the CFD grid via a matrix G

6X = GD. (6.2)

where 6X is the corresponding displacement vector on the CFD grid. Similarly the

surface forces can also be mapped from the CFD grid to the CSD grid.

Upon updating the aerodynamic forces in the structural equations and providing

the new boundary conditions to the fluid solver, the temporal development between

the aerodynamic and structural equations can be coordinated. To conduct time-

accurate computations for the fluid-membrane dynamics, iterations between the fluid

and structural solvers should be done at each time instant. A step-by-step coupled

algorithm is described below:

1. Generate the CFD and CSD grids.

2. Set n 0, evaluate pressure Po on the undeformed wing based on the initial

values po, no, and po.

3. n=n+l, and t,=t,_l + At.

4. Solve the Navier-Stokes equations and transfer the pressure P,+i to the struc-

ture based on the TPS interpolation.

5. Use the structural solver to evaluate the displacement vector D,+i.

6. Perform subiterations to synchronize the fluid and structural solvers, and set

i 0.

(a) i i+1;

(b) Map the displacements from the CSD grid to the CSD grid with the TPS

interpolation 6X+ = GD+ .

(c) Update the CFD grid with the moving grid technique.

(d) Update the Jacobian to satisfy the geometric conservation law.

(e) Update the boundary conditions for the fluid solver.

(f) Compute the pressure PT4+1 field based on the new CFD grid and boundary


(g) Map the aerodynamic load from the CFD grid to the CSD grid with the

TPS interpolation.

(h) Calculate the structural displacement D'i.

(i) C('! 1: for convergence for the subiteration process.

If 'No', return to step a), otherwise, continue.

7. Return to step 3 to process next time step.

It should be noted that a fully converged solution in the process is hard to

achieve. The number of subiteration is based on the computational resource, problem

nonlinearity, and the time step. In our study we find two or three subiterations are

enough to alleviate the phase lag errors.

6.3 MAV Wing Aerodynamics and Structural Response

As shown in Fig. 1-1 the membrane wing consists of a leading edge spar and

three chordwise battens on each side of the half wing. The membrane material is

bonded to the spar and battens. While the membrane is flexible and does not sustain

bending moment, the carbon fiber is rigid which provides support for the membrane

and can sustain bending moment. To fully model the membrane wing, ideally, one

needs to model both the rigid battens and the flexible membrane. Two principal

difficulties arise by doing that. First, the membrane does not bear bending moment,

it has three degree of freedoms at each node; the batten, if modeled as a beam, has

five degrees of freedom at each node. A special treatment is required to treat the

interface point which belongs to different regimes and has different degrees of freedom.

Second, the membrane is much thinner than the batten, the mass matrix associated

with these two different materials makes the assembled mass matrix ill-conditioned.

Given these factors, we treat the batten as a special membrane material with a larger

density while the other properties are the same; the density ratio between the batten

and the membrane is set to be three.

6.3.1 Computational Background

In our study only the membrane wing is considered while ignoring the fuselage

and propellers. A schematic geometry of the wing is shown in Fig. 6-1. The shaded

area shown corresponds to the fuselage which is fixed. The physical problem under

consideration is viscous flow over a wing in an unbounded domain. Based on the

freestream velocity of 10 m/s, the root chord Reynolds number is 9x104. Computa-

tions are first performed to assess the sensitivity of the computed results to the wing

location and the types of boundary conditions imposed to the outer flow boundaries.

Basically, the boundaries should be set sufficiently far away from the membrane so

that they do not have much impact on the computed aerodynamic coefficients. Based

on our tests a proper wing position is chosen and shown in Fig. 6-2. The upstream

inlet boundary is 7c from the leading edge, here c is the chord length at the root. At

surface ABCD a zero gradient boundary condition is imposed for the velocity compo-

nents; all the other boundaries are assigned as a Dirichlet type boundary condition.

Since the freestream velocity is parallel to the chordwise direction and no propeller

is modelled, only half of the domain is computed, and a symmetric mapping can be

applied to the other half domain.

Grid sensitivity tests are also performed. Three grid systems have been system-

atically chosen, ranging from 1.8 x 105 nodes at the coarse level and 2.3 x 106 nodes

at the fine level. For the coarse grid there are 41 points in the chordwise direction

and 31 points in the spanwise direction on the wing surface. The grid distribution



SI I Batten 3

Batten 2

Batten 1

Figure 6-1: Schematic wing with three battens on each half wing.

(9c, 5c, 6c

Membrane Wing

(-6c, -5c,

Figure 6-2: Computational set-up and boundary conditions.

at the coarse level around the wing is shown in Fig. 6-3. The results from the fine

grid are used as the reference. The comparison among these three grid systems is

summarized in Table 6-1. The coarse grid solution under-predicts the lift-to-drag

ratio (CL/CD) by 0.5'. and the difference in form drag is less than 0.7'. A similar

conclusion can be drawn for the lift, L. As expected, the skin friction drag, DF is

sensitive to the grid density; the coarse mesh overestimates DF by almost 10' An

intermediate grid, which has 6.7x105 nodes and more than doubles the surface points

on the wing surface comparing to the coarse grid, gives an improved prediction on

skin friction drag. From Table 6-1 it can be seen that skin friction drag has much

less contribution to the total drag than form drag. This is also valid at high angles

of attack where the form drag is the dominant factor in the total drag.

Figure 6-3: Schematic wing geometry and structured grid with 1.8x 105 points and
41 x 31 on half wing surface.

We further compare the main flow features. Flow separations are observed with

all these computational grids at the lower wing surface. We extract the chordwise

velocity profile at the quarter chord of the root region. The boundary 1 -i-r velocity

Table 6-1: Grid refinement effects on the computed aerodynamic coefficients.

Coarse 1.8x 10 6.24E-2 5.10E-1 9.70E-3 7.06
Intermediate 6.7x105 6.16E-2 5.08E-1 9.22E-3 7.16
Fine 2.3x 106 6.20E-2 5.04E-1 8.80E-3 7.10

profile is shown in Fig. 6-4. The fine grid has the largest reverse velocity, the coarse

grid has the smallest, and the intermediate grid has the thickest separation zone.

These results indicate that the grid resolution noticeably affects the sharpness of the

velocity profile, as well as the extent of the flow separation. Nevertheless, the coarse

grid predicts the main features of the flow which is qualitatively consistent with those

on the fine grid.




o 0

0.5- ,

1.5- .... ..... .

-5 0 5 10
Spanwise velocity (m/s)

Figure 6-4: Spanwise velocity profiles at the bottom surface at the quarter chord of
the root region for rigid wing at a= 6.

Comparisons between the intermediate and coarse meshes at other angles of

attack are also conducted. The computed aerodynamic coefficients are shown in

Fig. 6-5. Both the lift and drag coefficients show good agreements, with a little bit

larger difference at higher angles of attack. A few more words should be put for the

computational results at high angles of attack. Generally -I'" i1i:- vortex shedding

and boundary l-v-r separation occur at large angles of attack, and they introduce

unsteadiness to the aerodynamic performance and therefore affect the aerodynamic

coefficients. Cummings et al. [10] reported that at large angles of attack the unsteady

computations predicted much lower lift coefficients than the steady computations. We

perform steady computations when the incidence is less than 15 while conducting

both steady and unsteady computations when the angle of attack is larger than that.

At high angles of attack, truly steady solutions can not be obtained. The -I. "ly"

solutions are obtained with steady computations once certain order is reduced in the

residual. The PISO method discussed in C'!I pter 2 is used for the unsteady com-

putation. We set the time step to be 5x10-5s. However, for our computations, the

difference between the steady and the time-averaged lift and drag coefficients, are

found to be less than 1 over all the studied cases. Fig. 6-5 shows the computed

aerodynamic coefficients for the rigid wing with both coarse and intermediate meshes.

At low angles of attack, the coarse grid predicts slightly larger lift than the interme-

diate grid; the trend is reversed at high angles of attack. However, both predict that

the stall occurs approximately at 51.

2.2 :3
2 *, 2.l.r-r Gr,.J

0 l

0. 0
0 10 20 30 40 50 60 0 10 20 30 40 50 60
Angle of Attack Angle of Attack

Figure 6-5: Aerodynamic coefficients versus the angle of attack: (left) lift coefficient;
(right) drag coefficient.

For a problem as complicated as the present one, truly grid independent results

are difficult to attain. Considering the available computing resource and the out-

come of the aforementioned investigation, we use the intermediate grid for the rigid

wing simulation while using the coarse grid to illustrate the key issues in the fluid

and flexible structure interaction system, to demonstrate the salient features of mem-

brane wing computations, and to illustrate the related aerodynamic and structural


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]

CD CD,P+CD,F+ L-A (6.3)

where CD,P is the form drag coefficient due to the pressure, CD,F is the friction drag

coefficient due to viscosity, e is the span efficiency factor which is less than one, AR

is the aspect ratio, and
CD,i --r
i reAR

is the induced drag coefficient due to the existence of tip vortices. Eq. (6.3) demon-

strates that induced drag varies as the square of the lift coefficient. Furthermore, Eq.

(6.3) illustrates that as AR is decreased, induced drag increases. The MAV wing in

our study has a low aspect ratio of 1.4; therefore, it is important to investigate the

tip vortex effects on the wing aerodynamics. Besides their effects on the di tip

vortices have two-folded impacts on the lift:

Tip vortex causes a downwash component that decreases the effective angle of

attack and therefore decreases the lift [2].

Tip vortex forms a low pressure region on the upper wing surface, and provides

additional lift [67].

Fig. 6-6 shows tip vortices around the wing surface at different chordwise po-

sitions together with the streamlines at the angle of attack of 39" [53]. The higher

pressure from the lower surface drives the flow toward the upper wing surface where

the the pressure is lower. The two vortices on both sides of the wing demonstrate

clockwise and anti-clockwise rotations; it seems that the flow is leading from the

lower surface to the upper surface. Fig. 6-7 shows the low pressure zone due to the

vortical flow. The pressure drop further strengthens the swirl by attracting more flu-

ids toward the vortex core; meanwhile, the pressure decreases correspondingly in the

vortex core. The low pressure region created by the vortex generates additional lift.

Toward downstream, the pressure recovers to its ambient value. Then the swirling

motion weakens and the diameter of the vortex core increases. The vortex core losses

its coherent structure at downstream.

Fig. 6-8 visualizes the evolution of the vortical flow structure with the angle of

attack. The pressure distribution on the upper surface is also presented in the same

figure. At a= 6 tip vortices are clearly visible even though they only cover a small

area and are of modest strength. The flow is attached to the upper surface and follows

the chordwise direction. A low pressure region near the tip, which is indicated with

green color, is observed. Vortices strengthen with the increase of angle of attack.

At a 27, as shown in Fig. 6-8, tip vortices develop a strong swirling motion and

entrain the surrounding flow to the vortex core. The low pressure area increases as

the angle of attack becomes higher.

The pressure drop in the vortex core is demonstrated in Fig. 6-9a where the

spanwise pressure coefficient on the upper wing surface at x/c=0.4 is plotted. At

Figure 6-6: Streamline and vortices for rigid wing at a= 39. The vortical structures
are shown on selected planes.

o 10060


Figure 6-7: Pressure distribution around the rigid wing in the cross sections with
streamlines at angle of attack of 39.

(c) 270


\ C

(d) 330

(e) 450

(f) 510

Figure 6-8: Evolution of flow pattern for rigid wing versus angles of attack.


(a) 60

(b) 15"


a =6 the spanwise pressure is almost uniform on the upper wing surface, and the

pressure drop occurs at approximately 91I' of the half span from the root. With the

increase of the angle of attack, the pressure drop moves further towards the root. At

a 27, it occurs at the 7 .'- from the root. At lower angles of attack, the vortex core

position shows a linear relation with the incidence. This relation diminishes at higher

angles of attack when the flow is separated on the upper surface. For example, at

a=45, the flow is separated at the leading edge, and the low pressure zone covers

more than ,!11' of the wing surface, which helps to maintain the increase in lift. At

a=51, considerable spanwise velocity component is seen and the flow is separated

from most of the upper surface (Fig. 6-8). The lift by the vortices is not enough to

maintain the increase of the lift coefficient with the angle of attack, and stall occurs.

Even thought tip vortices affect the pressure distribution at the lower surface, the

effected region is mainly located near the tip. Most regions in the lower surface are

not affected by the tip vortices. From Fig. 6-9b it can be seen that the pressure

is almost uniform in the spanwise direction even at high angles of attack. Fig. 6-9

is illustrative in regard to the pressure distributions versus the vortical structures,

however, they are not indicative of the total level of the pressure force. One should

also note that Fig. 6-9 indicates that the moment experiences substantial variations

as the angle of attack changes.

The Laminar Boundary Lv,- r Separation on Rigid Wing

The laminar boundary lI,--', especially under low Reynolds number conditions,

is prone to separate under an adverse pressure gradient. This separation is usually

followed by a quick transition from the laminar to turbulent flow. If the adverse

pressure gradient is modest, the resulting turbulent flow, by obtaining energy through

entrainment, may reattach to the surface to form a laminar separation bubble and

then forms an attached turbulent boundary lV,--r.



o I | a=51
-e- a=45
=I .-2 -- a=270
-J =J- .' a=150
-25 a=60
S" -3-

0 02 04 06 08 1 350 04 06 08
z/Z z/Z
(a) c, at upper surface (b) c, at lower surface

Figure 6-9: Spanwise pressure coefficient distributions at x/c=0.4 for rigid wing at
different angles of attack.

Ever since its first observation by Jones [42], many experimental investigations

on the laminar separation bubble have been conducted, as reviewed by Young and

Horton [114]. Typically, a laminar separation bubble has the following characteris-

tics. Just downstream of the separation point, there is a "dead-air" region, where the

recalculating velocity is significantly less than the freestream velocity and the flow

is virtually stationary. The separated flow forms a free-shear li-.-r, which is highly

unstable. The separated shear 1i--r quickly undergoes transition from laminar to

turbulent flow. The turbulent flow may reattach to the surface behind a vortical

structure and forms a turbulent boundary l1-, -r. Hence, a separation bubble forms.

The dynamics of a laminar separation bubble depends on the Reynolds number, the

pressure distribution, the geometry, the surface roughness, and the freestream tur-

bulence intensity. A rough rule to determine the reattachment point is given by

Carmichael [7] that -zi- the Reynolds number based on the freestream velocity and

the distance from the separation point to the reattachment point is approximately

5x 104. On the one hand, if the Reynolds number is less than 5x 104, an airfoil will

experience separation without reattachment; on the other hand, a long separation

bubble will occur if the Reynolds number is slightly higher than 5x 04. In the fol-

lowing, detailed flow structures around a rigid wing are presented. By contrasting the

detailed fluid flow around the rigid and membrane wings, one can better evaluate the

MAV technologies. The bubble structures at different angles of attack are demon-

strated in Fig. 6-10. Representative instantaneous streamlines at the root section

on the rigid wing are plotted. At this moment, the spanwise velocity component is

ignored to make a clear representation of these separated structures. At a low angle

of attack of 6, the flow primarily attaches to the surface and a tiny separation bubble

is observed on the lower wing surface near the leading edge. Beginning at 45'. chord

from the leading edge, a weak recirculation zone is seen on the upper wing surface,

which produces a reattachment length of x/c=0.5.

With the increase of incidence, the separation point moves forward toward the

leading edge. At the angle of attack of 150, as shown in Fig. 6-10, the separation oc-

curs at I'. of the chord from the leading edge, which is followed by a long "dead-air"

zone before it reaches a maximal reverse velocity of 0.26U. In the "dead-air" zone, the

stationary reverse flow does not have much effect on the pressure distribution, which

is primary determined by the wing curvature. The shear 1-~.-r reattachment happens

at I' of the chord. Based on the freestream velocity and the distance between the

separation and reattachment points, the Reynolds number is approximately 4.5x 104

which is slightly smaller than the Reynolds number sI-:-, -1 .1 by Carmichael. The

reattachment point corresponds to a rapid increase in the surface pressure, it is highly

unstable, and moves forward and backward.

The streamwise velocity contours at different angles of attack are demonstrated

in Fig. 6-11. The velocity is normalized by the freestream velocity. At a= 6, the

maximum reverse velocity is less than 0.5'. of the freestream velocity. Under such a

condition, experiment shows that a considerable portion of the shear 1-,-r remains

laminar [21].

(a) 60

--- --------

(b) 150

(c) 270

(d) 390

(e) 510

Figure 6-10: Streamlines at different angles of attack for rigid wing. Shown are single
span shorts of the individual time dependent flows.

However, with the increase of angle of attack, the magnitude of the reverse

velocity increases as well. At a=27', the reverse flow component of the separation

bubble reaches a maximal velocity of 0.49U. Under such a situation, Crompton and

Barrett [9] showed that the shear 1~ ,-r is particularly energetic since most of the

shear 1lv.,r is turbulent. Even though the flow on the upper surface near the root

has separated from a large portion of the surface; however, the lift still increases with

the angle of attack. Two reasons might lie behind this. First, tip vortices generate

suction areas on the upper wing surface which provide additional lift. Second, even

though flow separates near the root, it still attaches to the other region of the wing

surface. When the angle of attack is increased further, at a=51", massive separation

appears, and stall occurs. Meanwhile, as seen from Fig. 6-12, vortex shedding is

observed near the root. A tiny separation bubble emerges at the leading edge first,

it then obtains energy from the shear l1v. r and increases its size. This bubble then

merges with another bubble downstream to form a larger bubble. The larger bubble

is not stable, it breaks into two small bubbles near the maximal camber position, and

vortex shedding begins.

For low aspect ratio wings, tip vortices have considerable contributions to the

lift. This scenario is quite similar to that of delta wings [23]. We have observed that

the low aspect ratio wing suffers less from the separation. As seen from Fig. 6-5a,

the lift keeps a linear relationship with the angle of attack even at very high angles

of attack. Experiments by Torres and Mueller [97] on low aspect ratio wings have

similar findings.

(a) 6
^\ 1.30



(b) 150

\ 1. 3


I 2.40
,, /f-"---------Y- I-OO

(d) 39

[] 0.45

(e) 51

Figure 6-11: Normalized chordwise velocity u/U contours at different angles of attack.

Figure 612: Vortex shedding at a51

Figure 6-12: Vortex shedding at a=51.

Unsteady Phenomenon at High Angles of Attack

Both vortex shedding and boundary 1-v-r separation occur at large angles of

attack, and they introduce unsteadiness to the aerodynamic performance. As men-

tioned before, we perform steady computations when the incidence is less than 15

while conducting both steady and unsteady computations when the angle of attack

is larger than that. Interestingly, the difference between the steady state computa-

tions and the time-averaged unsteady computations are small even at large angles of

attack in which unsteady phenomenon such as vortex shedding are prominent. The

differences in the lift coefficient and the lift-to-drag ratio are found to be less than

1 (Fig. 6-13). Noticeable difference appears only when the angle of attack becomes

sufficiently high.

-0- Time Averaged
-- Steady


0 10 20 30 40 50 60
Angle of Attack (degree)

Figure 6-13: Comparison of the rigid wing CL between the steady and unsteady

The pressure distributions are also compared at the root section, where flow

separation usually appears first due to the large camber in the present MAV wing

design [34]. Fig. 6-14 shows that at a= 6 the time averaged pressure coefficient

-0 -1
-04- -1
-06 12

0 02 04 06 08 1 0 02 04 06 08 1
x/c x/c

Figure 6-14: The cp comparison on the rigid wing at the root between the steady
and unsteady computations: (left) a=6; (right) a=15.

matches closely with the steady state result. This finding is consistent with Fig. 6

10 and Fig. 6-11 where a very weak recirculation is seen at this angle of attack.

However, at a=15, clear difference is shown after x/c=0.6 which approximately

corresponds to the location of the maximal reverse velocity in Fig. 6-11. The time

averaged value yields a smooth pressure distribution; the variation in the steady result

is apparent from the recirculation zone. Velocity distributions show the same trend

as the pressure distributions. At a= 6, as seen from Fig. 6-15, there is almost no

difference in the chordwise velocity distribution. However, at a=15, no difference

near the leading edge where no separation occurs, but clear difference is shown after

the separation bubble (Fig. 6-16).

6.3.3 Membrane Wing Dynamics

Time Synchronization

Before presenting the membrane wing results, some key issues in coupled membrane-

fluid interactions are first addressed. Specifically, the time synchronization between

the fluid and structural solutions, and the geometric conservation law in regard to

the moving grid method are emphasized.

Boundary layer profile at x/c=0 14 for a=60
Time Averaged
s-- Steady


0 0205





0 02 04 06 08

1 12 14

Figure 6-15: Comparison of chordwise velocity profiles on the rigid wing between

the steady and time-averaged unsteady computations at the angle of attack of 6.

Velocity profiles are at two chordwise locations, and both at the root section.

Boundary layer profile atx/c=0 19 fora=150
Time Averaged

1 12 14

Figure 6-16: Comparison of the chordwise velocity profiles on the rigid wing be-

tween the steady and time-averaged unsteady computations at angle of attack of 15.

Velocity profiles at at two chordwise locations, and both at the root section.


Boundary layer profile at x/c=0 80 for a=60



0 0205

> 002


Boundary layer profile at x/c=0 80 for a=150



0 02 04 06 08

The time synchronization in our study means that iterations should be enforced

between the fluid and structural solvers at each time step to avoid phase lag errors.

We shall simulate the fluid and structure interaction by integrating the fluid and

the structural solvers with subiteration. To demonstrate the importance of the time

synchronization between the fluid and structural solvers, We compare the results

with and without synchronization in Fig. 6-17 that depicts the displacement history

of a trailing edge node. The computation without synchronization fails to continue

while that with synchronization keeps on going. The displacement history without

synchronization matches well with that with synchronization at the very early stage.

However, the 1 ,-.ii-.; errors accumulate with time and eventually result in a much

larger displacement than that with synchronization. Gordnier and Visbal [23] have

also reported the importance of the time synchronization. The velocity and pressure

histories of the same node with and without time synchronization are compared in

Fig. 6-18 and Fig. 6-19, respectively. Before the computation diverges, the velocity

without synchronization is more than two times larger than that with synchronization.

The larger velocity of the membrane node means that the membrane node obtaining

more energy from the surrounding flow than that with synchronization. The high

velocity is believed to be associated with the sudden pressure drop shown in Fig. 6

19, which causes the wrinkle of the membrane and the failure of the structural solver.

The Geometric Conservation Law

The geometric conservation law [100] is another important factor when the mov-

ing grid technique is employ, l Fig. 6-20a shows the time history of CL/CD for the

membrane wing at a= 6. Without satisfying the geometric conservation law, irreg-

ular spikes are observed in the course of computations. However, the computations

are better regulated once the geometric conservation law is enforced as in Fig. 6-20b.

Farhat et al. [27] argued that the geometric conservation law was a necessary con-

dition to maintain the stability of a scheme utilized in moving boundary problems.

With Synchronization
Without Synchronization

Time (Second)


Figure 6-17: Effects of time synchronization on the displacement
point for membrane wing at a=6 .

of a trailing edge

With Synchronization
Without Synchronization


0.05 0.1 0.15
Time (Second)

Figure 6-18: Effects of time synchronization on vertical velocity
point for membrane wing at a= 6.


of a trailing edge

0.12 I

0.1 I

0.08 -







0 ,

With Synchronization
-2 Without Synchronization






0 0.05 0.1 0.15 0.2
Time (Second)

Figure 6-19: Effects of time synchronization on pressure distribution on a trailing
edge point for membrane wing at a=6.

We do not come across instability in our study. Two possibilities lie behind that:

one is the small time step we use for the fluid solver, the another is the relatively

small displacement in the membrane wing. Nevertheless, the computations behave

erratically when the geometric conservation law is not enforced.

Self-initiated Membrane Vibration

Earlier works in membrane wings [88, 92] were based on considerations that the

membrane was massless and there was no time dependent movement in the steady

freestream. Experiments [111] showed that the membrane experienced high frequency

vibration. By adopting a dynamic membrane model, our computations show that

the membrane wing also exerts high frequency vibrations in the steady freestream.

Fig. 6-21 demonstrates the displacement history of the trailing edge with time. The

maximal displacement during that period occurs near the tip. The displacement is

normalized by the maximal camber that is about 0.9 cm. To investigate the vibration

frequency, we choose a point between the batten 1 and batten 2 on the trailing edge

5 75
i_! i i. s, ,

0 55
-5r''' 41''
0 0.02 0.04 0.0B 0.08 D.1 0.12 0 005 01 015 02
Time (sec) Time (sec)

Figure 6-20: Effect of the geometric conservation law on CL/CD for membrane wing
computation at a= 6: (left) not satisfying the geometric conservation law; (right)
satisfying the geometric conservation law

as an example. Its deflection history is shown in Fig. 6-17, and the estimated dom-

inated frequency is around 120 Hz (Fig. 6-22). Analyses at other points also show

the frequencies around 120 Hz that is comparable to the experimental frequency of

140 Hz. [111] Considering that the experimental conditions and the detailed wing

configuration between our study and those reported by Waszak et al. [111] are not

identical, the qualitative agreement between computations and experiments in this

regard is deemed satisfactory. Liu [60] reported that under a typical wing gust situa-

tion the energy was mainly located in the low frequency range of several Hz or lower.

The membrane fluctuation is not expected to cause sensitive response to the vehicle

since the membrane fluctuates in a time scale much faster than either the vehicle

control scale or the expected wing gust time scale.

When the membrane surface vibrates, the velocity on the wing surface is no

longer zero. A considerable velocity component is observed at the wing surface.

The vertical velocity contours at two different time instants are plotted in Fig. 6-23.

Fig. 6-23a corresponds to a stage when the membrane wing moves up under the

aerodynamic forces. At the trailing edge, the vertical velocity is about S'. of the



0. 1

0 .
E 0.0c

CL 0.0-4


t=3 x 10-3
t=7 v 10-3
1=10 10-
I 1,. 10-
l=1. 10


0 0.2

Figure 6-21: Time history of trailing
The camber at the root is 0.90 cm.

0.4 0.6 0.8

edge displacement for membrane wing at a= 6.

x 10-3






50 100 150 200 250
Frequency (Hz)


350 400

Figure 6-22: Spectrum analysis of the trailing edge point vibration for membrane
wing at a=6.


v/U contour at t=3.0 x 10-3 v/U contour at t=7.0 x 10-3
0 0
; 0746 -0.000772
2 2 -0.000772
-,iIJw, 1
JGP 0. 00438
6 r 4 -4o.000772
J} 4
S/ 6 r, ,,, -0.000772
V \I I, I /
-2 02

210 12-0 1 '0

12- 12 00592

-6 -4 -2 0 2 4 6 -6 -4 -2 0 2 4 6
Spanwise (cm) Spanwise (cm)

Figure 6-23: Normalized vertical velocity contour at two different time instants for
membrane wing at a= 6 on the upper surface.

freestream velocity. A negative velocity component is also observed near the leading

edge. Fig. 6-23b approximately corresponds to a stage of maximal displacement. A

negative velocity near the root indicates that the membrane begins to move down.

Comparison between the Membrane and Rigid Wings

With the same freestream condition, as shown in Table 6-2, the flexible wing

exhibits slightly less lift coefficient than the rigid one at a= 6. The difference in

CL/CD is also small. At higher angle of attack of 15, the membrane wing generates

a lift coefficient about 2'. less than the rigid wing does; however, its CL/CD is

1.5'. more than the rigid one. These differences lie behind the high frequency of the

membrane vibration. Under the external force, the membrane wing changes its shape.

This shape change has two effects. On the one hand, it decreases the lift by reducing

the effective angle of attack of the membrane wing; on the other hand, it increases the

lift by increasing the camber. Our numerical findings show that membrane and rigid

wings exhibit comparable aerodynamic performance before stall limit, which was also

experimentally observed by Waszak et al. [111].

Fig. 6-24 shows the time averaged vertical displacements of the trailing edge at

a=6 and a=15. The displacements are normalized by the maximal camber of the

Table 6-2: Aerodynamic comparison between membrane and rigid wings.

Membrane Wing Rigid Wing
6 CL/CD 7.05 7.06
CL 0.530 0.532
15 CL/CD 3.94 3.88
CL 0.920 0.954

wing. At a =6 the deflection is about 15'. of the maximal camber and increases to

more than 211'. at a=15. Due to the trailing edge deflection the effective angle of

attack of the membrane wing is less than that of the rigid wing. The spanwise angles

of attack between rigid and membrane wings under the same flow condition and with

identical initial geometric configurations are shown in Fig. 6-25. As we can see from

Fig. 6-25a, the rigid wing has an incidence of 6" at the root and it monotonically

increases to 9.5" at the tip; the membrane wing shares the same angles of attack with

the rigid wing in the 3.' of the inner wing; however, the effective angle of attack

toward the tip is less than that of the rigid wing. At the tip, the angle of attack of

the membrane wing lowers by about 0.8". Fig. 6-25b compares the angle of attack at

a=15", and it shows that the effective angle of attack of the membrane wing is about

1 less than that of the rigid wing at the tip. The reduced effective angle of attack

causes the decrease in the lift.

The membrane wing shape change has another effect. The pressure difference

between the upper and lower wing surfaces inflates the membrane wing's camber,

which is shown in Fig. 6-26. Two airfoil shapes at different spanwise positions are

plotted together with the corresponding rigid wing shapes. The camber increase is

more visible in the outer wing than that in the inner wing. This is consistent with

Fig. 6-24 where the maximal displacement occurs near the tip. As expected, the

increase is larger at a=15" than that at a=6.

The wing surface movement affects the overall pressure distribution. Fig. 6

27 compares the time averaged pressure contours between the rigid and membrane


S01 -

Qi 5-_

o Rigid part

0 02 04

w 02



0 05
:z Rigid part

06 08 1 0 02 04

06 08

Figure 6-24: Averaged displacement of the membrane wing trailing edge: (left) a=6;

(right) a=15.

Rigid Wing 1ng
Membrane Wing Membrane Wng
95 185
9 18

65 55
6 15
0 02 04 06 08 1 0 02 04 06 08 1
Spanwise z/Z Spanwise z/Z

Figure 6-25: Comparison of the spanwise angles of attack between the rigid and

membrane wings: (left) a=60; (right) a=15.

wings. As expected, the differences are mostly located at the trailing edge where

large movement is observed. These differences are also reflected in the chordwise

pressure distributions shown in Fig. 6-28. At the root section there is almost no

difference between the rigid and membrane wings since the fuselage is fixed. How-

ever, the difference becomes clear toward the tip; at z/Z=0.37, the difference is

visible. At z/Z=0.69, the time averaged membrane wing shift the lowest pressure

point downstream comparing with the rigid wing result. This shift helps to delay

flow separation.


Spanwise zZ=O 40


0 02 04 06 08 1

Figure 6-26: Time averaged airfoil shapes at
membrane wing at a= 6 and 15.

Pressure Coefficient Contour

-0.299 -,n -0.5



-0.567 -0.


Spanwise (cm)

Figure 6-27: Comparison of c, contours at a=6:
wing; (right) steady rigid wing.

(left) time averaged membrane

02 04 06 08

different spanwise positions for the


Spanwise z/Z=O 60


Time Averaged Membrane
Rigid Wing

o- 0


0.2 0.4 0.6 0.8 1

Time Averaged Membrane
Rigid Wing


0.2 0.4 0.6 0.8




Time Averaged Membrane
Rigid Wing

0 0.2 0.4 0.6 0.8 1

Figure 6-28: C'! i Il -. c, comparison at different spanwise locations between time
averaged membrane wing and steady rigid wing at a=6: (a) z/Z 0; (b) z/Z=0.37;

(c) z/Z 0.69.


In C'! ipter 6 we have performed the aerodynamic study of the wing aerodynamics

of both the rigid and flexible wings. In this chapter we will develop a systematic

approach to enhance the membrane wing performance.

7.1 Introduction

It is desirable not only to understand the pertinent physical phenomena of the

membrane wing dynamics, but to improve its performance based on the understand-

ing. In C'! lpter 6 we have investigated both rigid and membrane wing aerodynamics.

The wing shape adopted is from an existing design. One question then arises natu-

rally: is that the best wing design, or, can we improve its performance. To answer

this question, it is worthwhile to pause and take a look at how the existing wing is de-

signed. The design of the wing is mostly based on observations and a two-dimensional

analysis tool to test the performance of different airfoils. Based on the tests, one air-

foil shape is then chosen to form a three-dimensional wing. The wing shape is then

tested with flight. In this trial and error process, the three-dimensional effects and the

membrane material properties are not considered. A systematic approach is necessary

to improve the wing performance.

There are many existing mature techniques for shape optimization. Could we

simply use one of them to achieve our goal? The answer is negative. As discussed

in Chapter 6, the MAV wing performs under the low Reynolds number condition,

and it exhibits a high frequency vibration even under steady freestream conditions.

Therefore, the optimization of a membrane wing introduces multiple challenges. First,

coupled, time dependent simulations of the interaction between viscous fluid flow

and flexible structure are very expensive. Second, an efficient and automatic grid

regeneration is essential. To deal with these difficulties, we propose the following


Adopt the rigid wing as a surrogate model for membrane wing optimization.

The surrogate approach was emploil,, in [6, 32, 51, 63], in our study the sur-

rogate model is the rigid wing. We shall optimize the rigid wing under a given

flow condition, and evaluate the performance of a flexible wing with output

from the surrogate model.

Automate a moving grid technique by treating the optimization process as a

iu1' img boundary piI, i, This technique facilitates not only the optimiza-

tion process but the computation of fluid and membrane wing interactions.

In the following, we first introduce the optimization problem. This is followed

by an optimization algorithm in which we show how to integrate the analysis tool,

the moving grid technique, and the optimizer Design Optimization Tools (DOT) [12]

into a system to do shape optimization. The optimization process begins with a rigid

wing surrogate optimization, then the flexible wing performance is assessed with the

optimized shape. Optimization results and the aerodynamic characteristics of both

optimized rigid and membrane wings are then highlighted.

7.2 Optimization and Computational Framework

The objective is to improve the lift-to-drag ratio of a membrane wing. The

baseline membrane wing shape selected is an existing design [33, 34] devised from

solutions provided by XFOIL [15] and empirical modification based on flight tests

[34]. XFOIL provides two-dimensional fluid flow solutions based on coupled thin l-iv r

and inviscid flow models. A schematic baseline wing shape is shown in Fig. 6-1 in

C'! lpter 6. The baseline shape has three carbon-fiber-made battens that are labelled

also in Fig. 6-1. The shading area corresponds to the fuselage. The fuselage is made

of rigid carbon fiber prepreg cloth which does not deform during the flight. A flexible

latex membrane covers the top surface which changes its shape under aerodynamic


The baseline wing has a variable spanwise camber. The nominal camber of the

baseline design is 7.5' at the root, and 2' at the tip. The wing has a maximal chord

length of 13.7 cm, a mean chord of 10.5 cm, and a span of 15.2 cm. C is used to

indicate the chord length at different span stations, and Z for the half span. The

angle of attack is defined in reference to the root geometry. At the design point, the

MAV flies with a speed of 10 m/s and an angle of attack of 6. The objective is to

maximize the lift-to-drag ratio at such a flight condition. At such a condition the

root chord Reynolds number is 9x 104.

Six design variables are chosen, they are the y-coordinates at six points on the

surface. Three are located at the root whose positions are approximately at 1(1'

7'. ^ and 77'. of the chord, as seen in Fig. 7-1. The other three are located at

Batten 2 with the same relative chord positions.

The wing surface can be implicitly represented by a function y = f(x, z), where

x, y, and z represent the wing surface coordinates. We indicate the design variable

with a capital letter Yi whose value is the perturbation of the y coordinate at point

i with respect to its baseline value. We interpolate the perturbations from the six

design variables over the wing surface with the thin plate spline (TPS) interpolation

method discussed in C'!i lter 6. Therefore, the interpolated perturbations over the

wing surface is

6y(x,z) G(Y), (7.1)

where Y is a vector indicating the design variable values, G is the interpolation

operator, and by is the interpolated value over the wing surface. The perturbed wing

surface is therefore represented by the summation of the baseline and the interpolated

value in Eq. (7.1), ynew baseline + 6y. One advantage of this perturbation approach

is that we can recover the original shape when the perturbation vector Y is zero. In






0.5- (a)


0 0.2 0.4 0.6 0.8 1


Design point 3


0.5- (b)
Baseline Shape
--- Perturbd Shape


0 0.2 0.4 0.6 0.8 1

Figure 7-1: A cross section of membrane wing along a batten: (a) straight lines
indicating how convexity constraint is applied; (b) effect of y coordinate perturbation
at a point. The design points are located at 1C('. 27'. and 77' of the chord.

our optimization process, both the leading and trailing edges are fixed to maintain

the identical angle of attack between the baseline shape and the optimized shape.

The optimization problem investigated can be formulated as follows.

Maximize CL/CD

Subject to

1: CL > CLbaseline (7.2)
Yi + y o yo 0 Y2 +y2
2: Convexity constraint: >
X1 Xo X2 X0
Constraint 1 requires that the lift coefficient be no less than that of the baseline wing.

Constraint 2 maintains the convexity of the airfoil (see Fig. 7-la) so to eliminate ob-

viously inappropriate shapes. Only one of the six constraints is shown for illustration

purpose. Constraint 3 gives the lower and upper bounds of the design variables whose

values are listed in Table 7-1.

Table 7-1: Design variables bounds and their optimal values at 6.

Initial Lower Upper Case 1: opti- Case 2: optimal
Parameter design bound bound mal values with values with 1000
(cm) (cm) (cm) 600 iterations per iterations per
analysis analysis
X1 0.00 -0.40 0.00 -0.258 -0.259
X2 0.00 -0.40 0.00 -0.400 -0.400
X3 0.00 -0.20 0.20 -0.152 -0.139
X4 0.00 -0.20 0.20 -0.100 -0.095
Xs 0.00 -0.10 0.20 -0.082 -0.068
X6 0.00 -0.10 0.20 0.135 0.121
Analysis 97 108
Cycle 6 6
CL 0.530 0.529 0.529 0.529
CL/CD 7.06 7.55 7.55
Number of analyses with convexity violations 13 22

Full Text


First,Iwouldliketoexpressmysinceregratitudetomyadvisor,Dr.WeiShyy.Heprovidedagoodstudyenvironmentandoeredplentyofopportunitiesformetoexploreandpursuemyresearchinterest.Iwouldliketothankhimforhisenduringenthusiasmineducatingmeasaresearcherandaperson.IhopethatIemulatehisneexampleofasuccessfulinternationalstudent,anintelligentandindustriousscholar,aresponsibleandknowledgeableadvisor,andanoptimisticandhumorousperson.Iowethankstoseveralcolleaguesfortheirwillingnesstosharetheirexperiencesandknowledgewithme.Ihavebenetedmuchfromtheircollaborationinvariousareas:Dr.PeterG.Ifju,microairvehicle;Dr.AndrewKurdila,structuraldynamics;andDr.RaphaelT.Haftka,optimization.Dr.SiddharthThakurgenerouslysharedhistimeandexperience,whichmademyworkmucheasier.Ienjoyedthegoodrelationshipswithmyprofessors.Dr.RenweiMeigavememanyhelpfulsuggestionsandimprovedmyknowledgeofadvanceduidmechanics.Dr.MarkSheplakspentmuchtimegatheringclassmaterialsforme.BoththeUniversityofFloridaandtheMechanicalandAerospaceEngineeringDepartmentprovidedapleasantenvironmenttostudyandlive.Iamproudtohavebeenhere;Ihopethattheywillbeproudtohavemehere.MyfamilymembersinChinahavegivenmetremendoussupportformystudyabroad.Theirtrustandlovehavesupportedmeovertheyearsandwillsupportmeinthefuture.Mywife,LihuiBai,hasbeenwithmeandsupportedmeallthetime.Herpatience,understanding,andencouragementaretheinvaluablewealthofmylife.ButnoacknowledgementcouldpossiblystateallthatIowetoher. iv


page ACKNOWLEDGMENTS ............................. iv ABSTRACT .................................... vii CHAPTER 1INTRODUCTION .............................. 1 1.1MAVandMembraneWingAerodynamics ............. 1 1.2MotivationsandChallenges ..................... 4 1.3Outline ................................. 5 2THEORETICALANDCOMPUTATIONALMODELSFORFLUIDFLOWS .................................... 8 2.1GoverningEquationsandFluidSolver ............... 9 2.1.1GoverningEquationsforFluidFlow ............. 9 2.1.2FluidFlowSolvers ....................... 12 2.1.3GeometricConservationLaw ................. 12 2.2TurbulenceModels .......................... 13 2.2.1ReynoldsEquationsandk-"TurbulenceModel ....... 13 2.2.2BoundaryConditionsforTurbulenceModels ........ 15 2.3Laminar-turbulentTransition .................... 18 3MOVINGGRIDTECHNIQUE ....................... 19 3.1PerturbationMethod ......................... 20 3.2Master/SlaveConcepts ........................ 21 3.3NumericalTests ............................ 26 4ACTIVEFLOWCONTROLWITHANADAPTIVELY-SHAPEDWING 27 4.1Introduction .............................. 27 4.2ExperimentalandComputationalSet-up .............. 28 4.3NumericalStudy ........................... 29 4.3.1GridSensitivityAnalysis ................... 29 4.3.2AssessmentoftheTurbulenceModel ............. 30 4.3.3EectsofGapsbetweenFlaps ................ 30 4.3.4EectsofFlappingAmplitude ................ 32 4.4Conclusions .............................. 34 v


.................. 36 5.1Introduction .............................. 36 5.2StructuralSolver ........................... 38 5.2.1Mooney-RivlinModel ..................... 38 5.2.2PrincipleofVirtualWork ................... 39 5.2.3TriangularMembraneElement ................ 40 5.2.4Green-LagrangeStrainTensor ................ 42 5.2.5InternalForces ......................... 43 5.2.6ExternalForces ........................ 45 5.2.7DynamicEquationsfortheMembrane ............ 46 5.2.8SolutionoftheDynamicEquations ............. 47 5.3NumericalTestofMembraneModel ................. 48 6MEMBRANEWINGAERODYNAMICSFORMAVS .......... 53 6.1Introduction .............................. 53 6.2CouplingbetweentheFluidandStructuralSolvers ........ 55 6.3MAVWingAerodynamicsandStructuralResponse ........ 58 6.3.1ComputationalBackground .................. 59 6.3.2RigidWingAerodynamics .................. 64 6.3.3MembraneWingDynamics .................. 76 7WINGSHAPEOPTIMIZATION ...................... 88 7.1Introduction .............................. 88 7.2OptimizationandComputationalFramework ............ 89 7.3ResultsandDiscussion ........................ 93 7.3.1SensitivityAnalysis ...................... 93 7.3.2PerformanceoftheOptimizedRigidWing ......... 95 7.3.3FlexibleWingOptimization ................. 97 7.4SummaryandConclusions ...................... 102 8CONCLUSIONSANDFUTUREWORK .................. 104 APPENDIX AAPPENDIXSTRAIN-DISPLACEMENTMATRIX ............ 107 REFERENCES ................................... 109 BIOGRAPHICALSKETCH ............................ 118 vi


Microairvehicles(MAVs),withwingspanof15cmorlessandightspeedaround10m/s,havemanyapplicationsinbothcivilianandmilitaryareas.TheReynoldsnumberbasedonthegivenparametersisaround104,whichoftenyieldsinsucientlift-to-dragratio.Furthermore,oneexpectstheunsteadyeecttobenoticeableforsuchightvehicles.Theexiblewinghasbeendemonstratedtoexhibitfavorablecharacteristicssuchaspassiveadaptationtotheightenvironmentanddelayedstall. ThepresentstudyfocusesondevelopingcomputationalandmodelingcapabilitiestobetterunderstandtheMAVaerodynamics.Bothexiblewings,utilizingmem-branematerials,andadaptively-shapedwings,utilizingpiezo-actuatedaps,havebeenstudied.Intheadaptively-shapedwingstudy,weusepiezo-actuatedapstoactivelycontroltheow.Weassesstheimpactsoftheapgeometry,appingam-plitude,andturbulencemodelingontheowstructurewithaparallelexperimentaleort.Themembranewingusesapassivecontrolmechanismtodelaythestallangleandtoprovideasmootherightplatform.Ourstudyfocusesonthemutualinterac-tionsbetweenthemembranewinganditssurroundingviscousow.Wecomparethe vii


Besidestheaerodynamicstudy,wealsoperformshapeoptimizationtoimprovethemembranewingperformance.Sincedirectoptimizationofamembranewingistootimeconsumingtobepractical,weoptimizeasurrogaterigidwingmodelbasedonanintegratedoptimizationalgorithm,whichconsistsofaNavier-Stokessolver,anautomaticgridgenerationtool,andagradient-basedoptimizer.Then,weassessthemembranewingperformancebasedontheoutcomefromthesurrogatemodel.Ournumericalresultsconrmthatthemembranewingexhibitsconsistentimprovementinthelift-to-dragratiowiththesurrogatemodel. viii


1.1 MAVandMembraneWingAerodynamics Microairvehicles(MAVs)withamaximaldimensionof15cmorlessandaightspeedaround10m/sareofinteresttobothmilitaryandcivilianapplications.Equippedwithavideocameraorasensor,thesevehiclescanperformsurveillance,reconnaissance,targeting,andbiochemicalsensingataremoteorotherwisehazardouslocation[ 83 ].Withtherapidprogressmadeinstructuralandmaterialtechnologies,miniaturizationofpowerplants,communication,visualization,andcontroldevices,severalgroupshavedevelopedsuccessfulMAVs[ 24 34 ].Asitssizeisreduced,MAVgainsfavorablescalingcharacteristics,includinglowinertiaandreducedstallingspeed[ 83 ].However,intermsofaerodynamics,control,range,andmaneuverability,manyoutstandingissuesarestillpresent.First,anMAViesinalowReynoldsnumberregime(104105)duetoitslowightspeedandsmalldimension.Suchaightenvironmentisoftenaccompaniedbylaminarboundarylayerseparation,transition,andlowlift-to-dragratio.Second,anMAVwingtypicallyhasalowaspectratio,whichcausesstrongvorticalowstructuresandincreasestheinduceddrag.Third,anMAVissusceptibletorollinginstabilities,whichbecomesevenmoreseriousbecauseoftheexistenceoftipvortices.Fourth,anMAVissensitivetowinguctuations,whichcanbecomparabletoMAV'sightspeedandmakesboththeinstantaneousightReynoldsnumberandangleofattackvarysubstantially[ 83 ]. InthedevelopmentofMAVs,twoapproachesexistintheMAVwingdesign:appingwingsandxedwings.Intheappingwingarea,pioneeringworkwasdonebyWeis-Fogh[ 112 ],andLighthill[ 57 ].RecentexperimentalandnumericalstudiesinthisaspectweregivenbyChasmanandChakravarthy[ 8 ],DeLaurier[ 11 ],Smith 1


[ 90 ],Dickensonetal.[ 13 ],Ellington[ 17 ],JonesandPlatzer[ 43 44 45 ],Katz[ 50 ],Kamakotietal.[ 48 ],LiuandKawachi[ 59 ],VestandKatz[ 106 107 ],Walker[ 109 ],andWang[ 110 ].Shyyetal.[ 83 ]gaveanextensivereviewofappingandxedwingcharacteristics.Templin[ 96 ]presentedtheappingwingspectrumofightanimals. Ourstudyfocusesonaxedwing.Specically,wewillstudyaexiblemembranewingforMAVapplications.Figure 1{1 showsa15cmMAVwithaexiblewingdesignedbyIfjuandcoworkers[ 34 ].TheMAVhasaightspeedbetween20and40km/h.Poweredbyanelectricmotor,itcanyforabout15minutes.Withabuild-invideocameraandatransmitter,thevehiclehasatotalmassaround40g.TheMAVhasarigidfuselageandleadingedgesparwhicharemadeoflaminatematerial.Oneachsideofthewingtherearethreeunidirectionalcarbonberbattens.Thesparandbattensarecoveredwithalatexmembranematerial.Overall,thewinghasaspanof15.2cm,arootchordof13.2cm,ameanchordof10.5cm,andawingareaof160cm2,whichresultinalowaspectratioof1.4.Thewinghasavariablecamberdecreasingfrom7.5%attherootto2%atthetip.Inthisconguration,themaximumcamberislocatedatabout27%chordattheroot. Figure1{1: MAVwithmembranewing.


Flexiblewingshavedierentaerodynamicperformancefromrigidwings.Oneadvantageofexiblemembranewingsisthattheyfacilitatepassiveshapeadaptationandresultinadelayedstall[ 83 111 ].Fig. 1{2 ,adoptedfromWaszaketal.[ 111 ],comparestheliftcurvesversusanglesofattackforrigidandmembranewingswithcongurationssimilartothoseshowninFigure 1{1 .Undermodestanglesofattack,therigidandmembranewingsdemonstrateasimilarliftcurveslopeof2.9,withtherigidwinghavingaslightlyhigherliftcoecient.However,themembranewingshaveamuchhigherstallanglesofattackthantherigidwing,whichisakeyelementinenhancingthestabilityandagilityofMAVs.Therigidstallsat24owhiletheexiblewingshavestallanglesbetween30oand49o,whicharesimilartothatofmuchloweraspectratiorigidwings(AR=0.5to1.0).But,theverylowaspectwingsexhibitlowerliftcurveslopesof1.3to1.7insteadof2.9.Hence,exiblewingsappeartocombinethedesirableperformanceofrigidwingswithhigherandloweraspectratiosbyexhibitingstallbehaviorsimilartothatofrigidwingswithaspectratioof0.5to1.0andtheliftgeneratingcapabilitysimilartorigidwingswithaspectratioof2.0[ 111 ].Anotheradvantageofexiblemembranewingsisthattheycanadapttothewindgustandprovidesmootherightplatforms,whichareveriedbywindtunnelexperiments[ 84 ]. Figure1{2: Liftcoecientversusangleofattack.(FromWaszaketal.[ 111 ])


1.2 MotivationsandChallenges Thoughexperimentshaveprovidedinformationonthemembranewingaero-dynamics,tofullyunderstandthemembraneaerodynamicsandtoappreciatethemembranewing'spassivecontrolmechanism,weneedtoperformadetailednumer-icalstudy.Anumericalstudyofthemembranewingbringstwochallenges:thelowReynoldsnumberandthelowaspectratiocondition,andtheinteractionbetweenthemembranewinganditssurroundingviscousows. First,thelowReynoldsnumberconditionpresentstremendouschallengesonMAVsstudy.Inourstudy,thechordReynoldsnumberoftheMAVisaround9104.IntherangeofReynoldsnumberof104106,complexowphenomenaoftenoccur.Thelaminarboundarylayerseparation,transition,andreattachmentusuallycoexistontheupperwingsurface.Anaccuratepredicationofthetransitionfromlaminartoturbulentstateisimportanttoevaluatethewingperformance.However,atransi-tionmodelforgeneralapplicationisnotyetavailable.Furthermore,lowaspectratiowingisusuallyaccompaniedbycomplicatedvorticalowstructures.Thevorticalowstructuresintroduceunsteadyphenomena.ValuableinformationregardinglowReynoldsnumberandlowaspectratiowingswasreviewedbyLissaman[ 61 ]andTani[ 95 ],aswellasaddressedinseveralproceedings[ 65 66 ].Gad-el-Hak[ 20 ]discussedapproachestoprevent/alleviateowseparationonMAVswithmicroelectronicsme-chanicalsystems(MEMS).LianandShyy[ 53 ]performedadetailednumericalstudyoftheowseparationontheMAVwinganditsimpactsonaerodynamicperformance. Second,amembranewingexhibitsself-initiatedvibrationsevenundersteadystatefreestreamcondition.Waszaketal.[ 111 ]conductedwindtunnelexperimentsandfoundthatthemembranewingsexhibitedvibrationatO(102)Hzinasteadyfreestream.SimilarobservationwasreportednumericallybyLianandShyy[ 53 ].Suchvibrationsandtheassociatedshapedeformationchangethepressuredistributionon


themembranewing,which,inturn,aectsthemembranedynamics.Thisrecipro-calactionleadstoauidandstructureinteractionproblem.Aeroelasticanalysismethodcouplingacomputationaluiddynamics(CFD)solverwithacomputationalstructuraldynamics(CSD)solverprovidesaneectivetoolforthemembranewingstudy.AlthoughbothCFDandCSDhaveachievedgreatsuccessindividually,acou-pledcomputationisstillachallengingtask[ 80 ].Todate,thestudyofamembranewinganduidowinteractionsislimited.JacksonandChristie[ 38 ]adoptedathree-dimensionalpotentialow-basedsolverandamembranewingmodeltoanalyzetheaeroelasticbehaviorofmarinesails.SmithandShyy[ 92 ]andShyyetal.[ 88 ]pre-sentedacomputationalapproachtomodeltheinteractionbetweenatwo-dimensionalexiblemembranewingandsurroundingviscousows.Acombinednumericalandexperimentalanalysisofaowoveratwo-dimensionalexiblesailwasconductedbyLorilluandHureau[ 62 ].Inthecoupleduidandmembranewinginteractionstudy,weneedauidsolvertoaccuratelycapturethetransientphenomenonoftheuidowsurroundingthemembranewing;weneedastructuralsolvertopredictthein-herentnonlinearbehaviorsofthemembranematerial.Tocouplethesetwosolvers,wealsoneedtohandledierenttimecharacteristicsassociatedwitheachsolver.Be-sidesthese,weneedamovinggridtechniquetofacilitatetheuidandstructureinteractionproblem.Aninterfacetechniqueisalsorequiredtoexchangeinformationbetweenthesetwosolvers. 1.3 Outline Thisworkstudyboththemembranewingandadaptively-shapedwing.First,wewilldevelopacomputationalcapabilitytostudytheuidandstructureinteraction.Second,basedonthiscapability,wewillstudythepiezo-actuatedappingwingthemembranewing.Third,withthegainedknowledge,wewilloptimizethemembranewingperformancewithmodernoptimizationtechnique.Theoutlineofthecurrentworkisasthefollowing:


InChapter 2 wepresenttheuidowsolverandphysicalmodelsoftheuiddynamics.TheNavier-Stokesequationsforincompressibleowoncurvilinearcoordi-natesareemployed,alongwithatwo-equationk-"turbulencemodel.Inthischapter,theuidsolverisconsideredunderthemovingboundaryframework.Abriefin-troductionisgiventothegeometryconservationlaw[ 100 ]thatplaysanimportantroleinmovingboundaryproblems.Eventhoughthepresentworkdoesnotemploylaminar-turbulenttransitionmodel,abriefreviewisoered. InChapter 3 wedevotetomovinggridtechnology,whichisusedbyboththecoupleduid-structurestudyandtheoptimization.Arobust,ecientandeectivemovinggridtechniqueisproposed,whichusestransniteinterpolation-basedpertur-bationmethodinsingleblockgrid[ 56 77 ],andemploysthemaster/slaveconcept[ 29 ]foramultiblockgrid. InChapter 4 weinvestigatetheactiveowcontrolwithadaptively-shapedapsforalow-Reynoldsnumberwing.Resultsoftheinteractionbetweenpiezo-actuatedapsandsurroundingowarepresented. InChapter 5 wefocusonthestructuralsolver,morespecically,anonlineardynamicmembranemodel.Detailedderivationofthemembranemodel,withanemphasisonthenonlinearityaswellasthetimeintegrationwillbepresented.Themembranemodelwillbeassessedusingwelldenedtestproblems. InChapter 6 weproposeacoupleduidandstructureinteractionalgorithm.Thecouplingisrealizedthroughaninterfacetechniqueandtimesynchronization.Thenumericalimplementationofgeometryconservationlawisgivenandtested.ThelowReynoldsandlowaspectratiowingaerodynamics,includingthelaminarboundarylayerseparation,tipvortices,andtheunsteadyphenomenonareinvestigated.Com-putationalresultswillbepresentedtoelucidatethemembranewingaerodynamics.


InChapter 7 weoptimizeofaexiblemembranewingwithagradient-basedmethod.Sinceadirectoptimizationofamembranewingrequiressubstantialcom-putationalresourceandtime,arigidwingisusedasasurrogate.Themembranewingperformanceisthenevaluatedbasedonthesurrogatemodel. Inbrief,thepresentworkismotivatedby,rst,thepracticalneedtounderstandtheaerodynamicsofmembraneandadaptively-shapedwings,second,ourinterestsinmicroairvehicles.Originalworkincludeamovinggridalgorithmformulti-blockgrid,adynamicmembranemodel,numericalstudyofinteractionbetweenuidandaexiblewing,andshapeoptimizationofaexiblewing.


Ourstudyinvestigatestheaerodynamicsofamembranewingaswellasadap-tivelyshapedwingsunderlowReynoldsnumberconditions.Duringight,themem-branewingundergoesshapechangeundertheexternalforces;meanwhile,itsshapevariationaectsthepatternandstructureofitssurroundinguidow.Themem-branewingvibrationwasobservedbothexperimentally[ 111 ]andnumerically[ 53 55 ].Intheadaptivelyshapedwingstudy,thepiezo-actuatedapsvibratedwithdierentfrequenciesandamplitudes.Theapvibrationscauseunsteadyphenomenalikevor-texshedding[ 30 56 ].Itisdesirabletohaveacomputationalcapabilitytoaccuratelycapturethetransientbehaviorofuidowandstructuraldynamics. Eventhoughtsteadyowcomputationshaveachievedsignicantsuccessintheirdevelopmentandapplications,littleattentionhasbeengiventothesystematicanaly-sisofunsteadycomputations[ 79 80 ].Severalfactorsareresponsibleforthissituation.First,steadyCFDisstillthemainstreamaerodynamicanalysis,andthetechniquesusedtoenhancetheperformanceofsteadyCFDcomputationsmaynotbesuitabletounsteadycomputations.Forexample,insteadycomputations,avarietyoftech-niqueshavebeenproposedtoeliminatethelowfrequencytransientstoacceleratetheconvergencetosteadystate;however,lowfrequencyphenomenaareintrinsictomanyuidandstructureproblems,andsuchtechniqueswillbeinappropriatetotheunsteadyanalysis[ 80 ].Second,theapplicabilityoftheturbulencemodeltounsteadyowsisnotwellunderstood[ 79 ]. However,numerousapproacheshavebeenproposedforunsteadyowcomputa-tions.Pulliam[ 74 ]incorporatedsubiterationtoimprovetimeaccuracyofconven-tionalimplicitschemes.Jameson[ 39 ]andRumseyetal.[ 79 ]developedadualtime 8


steppingtechniqueinthecontextofmultigridmethodologytoenhancetheeciencyandaccuracyoftraditionalfactoredschemes.Issa[ 35 ],andOliveiraandIssa[ 71 ]pro-posedthePressure-ImplicitwithSplittingofOperators(PISO)methodtoimprovetheconvergence.Forrelevantinformationoftheseandsomeofthemostfrequentlyadoptedapproaches,seeShyyandMittal[ 85 ]. TheSIMPLEmethod[ 72 82 ]andthePISOmethod[ 35 71 98 ]aretwopopularmethodstosolvethethree-dimensionalNavier-Stokesequationsforincompressibleows.Bothmethodsbelongtothepressure-basedapproach[ 82 ]whichdevisesanarticiallyderivedpressure(correction)equationbymanipulatingthemasscontinuityandmomentumequations. Three-dimensionalNavier-Stokesequationsforincompressibleowincurvilinearcoordinatesareadoptedinourstudy.Twomethodsareemployedtosolvethegov-erningequations,oneistheSIMPLECmethod[ 102 ]thatisbelongtotheSIMPLEfamily,andtheotheristhePISOmethod[ 35 ].Inthischapterwebrieyintroducetheseuidsolversunderthemovingboundaryframework.Fordetaileddescriptionsofthesemethods,werefertoRefs.[ 35 72 82 98 102 ].Tosolvethemovingbound-aryproblems,wealsoemphasizetheimportanceoftheboundaryconditionsandthegeometricconservationlaw[ 100 ].Anintroductionofturbulencemodels,withanemphasisonthek-"turbulencemodelispresented.Thoughanaccuratetransitionmodelforgeneralapplicationisbeyondthescopeofourstudy,weoerabriefreviewofthetransitionmodelsattheendofthischapterforcompletenesspurpose. 2.1 GoverningEquationsandFluidSolver 2.1.1 GoverningEquationsforFluidFlow Thethree-dimensionalNavier-Stokesequationsforincompressibleuidowsincurvilinearcoordinates[ 82 ]instrongconservativeformcanbewrittenasthefollow-ing, @+@V @+@W @=0;(2.1)


@[ J(q11u+q12u+q13u)]+@ @[ J(q21u+q22u+q23u)]+@ @[ J(q31u+q32u+q33u)][@ @(f11p)+@ @(f21p)+@ @(f31p)]+G1(;;):J; @[ J(q11v+q12v+q13v)]+@ @[ J(q21v+q22v+q23v)]+@ @[ J(q31v+q32v+q33v)][@ @(f12p)+@ @(f22p)+@ @(f32p)]+G2(;;):J; @[ J(q11w+q12w+q13w)]+@ @[ J(q21w+q22w+q23w)]+@ @[ J(q31w+q32w+q33w)][@ @(f13p)+@ @(f23p)+@ @(f33p)]+G3(;;):J; whereu,v,andwarethevelocitycomponentsinthex-,y-,andz-directions,re-spectively,istheviscosity,G1,G2,andG3arethebodyforcecomponentsperunitvolumeinthex-,y-,andz-directions.Thecoecientsappearingintheequations


areexpressedasfollows: Theexpressionsforfij,i=1,3,j=1,3,canbefoundin[ 82 ].TheJacobiantransfor-mationalmatrixisdenedasfollows: ThedeterminantoftheJacobianmatrixisexpressedas Inthegoverningequations,Uistheuxthroughthecontrolsurfacewhosenormaldirectionis;VandWaretheuxesthroughtheandcontrolsurfaces.Whenthegoverningequationsareconsideredunderamovinggridframework,thegridvelocitiesshouldbeconsideredintheuxcomputations.Incurvilinearcoordinates,theyareexpressedasthefollowing: where,forexample,_xisthegridvelocitycomponentinthex-directionthatcanbeapproximatedbytherstorderEulermethod, _x=xx0 wheretisthetimestepandthesuperscript0referstotheprevioustimelevel.Similarly,thevaluesoftheothertwocomponents_yand_zcanbeestimated.


2.1.2 FluidFlowSolvers TheSIMPLECmethod[ 102 ]isavariantoftheSIMPLEmethodoriginallypro-posedbyPatankar[ 72 ]inCartesiancoordinates,anditsextensiontocurvilinearcoordinateswasdiscussedbyShyy[ 82 ]toaccommodatethecomplexgeometrycom-putation.TheSIMPLECmethodusescoordinatedunder-relaxationsforthemomen-tumandpressurecorrectionstoimprovetheconvergencethatisintrinsicallyslowintheoriginalSIMPLEmethod.TheSIMPLECmethodwasalsousedtosolvemov-ingboundaryproblems[ 53 88 ].ContrarytotheSIMPLEfamilymethod,thePISOmethodisanon-iterativeapproachtohandlethepressure-velocitycouplingwhichsplitstheprocessofsolutionintoaseriesofpredictorandcorrectorsteps.ThePISOmethodisgenerallymoreecientfortransientowcomputations[ 35 98 ]. TosolvetheNavier-Stokesequations,werequireproperboundaryconditionsforthecomputations.Attheinterfacebetweentheuidandthesolidstructure,thefollowingboundaryconditionsareapplied, Eq.( 2.10 )expressesthecompatibilitybetweenthewetuidsurfacexfandthestruc-turesurfacexsattheinterface;Eq.( 2.11 )requiresthattheuidvelocityattheinterfacematchesthesurfacevelocityofthestructure.Theseboundaryconditionsareprovidedbycomputingthedisplacementofthesolidstructure. 2.1.3 GeometricConservationLaw Whenbody-ttedcurvilinearcoordinatesareusedinthecomputation,atrans-formationmapsaphysicalowregion(x,y,z)ontoauniformcomputationalspace(,,)whereconservationlawsarecarriedout.Tofacilitatethecoordinatetrans-formation,theJacobianmatrixJisintroduced.TheJacobianmatrixisdenedas


inEq.( 2.6 ).ThedeterminantofJrepresentsavolumeelementinthetransformedcoordinates.Inamovinggridproblem,thecomputationalgridneedstobeupdatedwithtime;meanwhile,theJacobianJneedstobeupdatedsimultaneously.SpecicproceduresarerequiredtocomputetheeectivevalueofJ;otherwise,errorsariseduetotheinconsistencyinevaluatinggeometricquantitiesandtheaddingofnon-physicalquantities.AspointedbyThomasandLombard[ 100 ],intheupdatingprocess,thefollowinggeometricconservationlaw(GCL)neededtobesatised, dtZVJddd=ZV(rWs)ddd;(2.12) whereVisthevolumeboundedbyaclosesurfaceS,WsdenotesthelocalvelocityoftheboundarysurfaceS.ThomasandLombard[ 100 ]proposedtoevaluateJwhilemaintainingthegeometricconservationlawbysetting=1,V=0fromthecontinuityequation,i.e., Implicationsofthegeometricconservationlawwithuid-structureinteractionprob-lemswerediscussedbyShyyetal.[ 88 ]. 2.2 TurbulenceModels 2.2.1 ReynoldsEquationsandk-"TurbulenceModel Forsimplicity,theReynoldsaveragedmomentumandcontinuityequationsforincompressibleowsarepresentedinCartesiancoordinates whereTijisthestresstensor,repeatedindicesinanytermindicateasummationoverallthreevaluesoftheindex,andUiistheReynoldsaveragedmeanvelocity.The


stresstensorinaturbulentowmaybewrittenas 2(@Ui k= Inmanyows,thesenormalstressescontributelittletothetransportofmeanmo-mentum.Theo-diagonalcomponentsofijareshearstresseswhichplayadominantroleinthemeanmomentumtransferbyturbulentmotion. Eq.( 2.14 )containsnineunknowncomponentsofij(ofwhichsixareindependentofeachother)besidesPandthethreecomponentsofUi.Therefore,weneedmoreequationsforijtoclosethesystemofequations.Aturbulenceclosuremodelservesthatpurpose. Dierentversionsofturbulenceclosuremodelshavebeenproposed.Thesuccessofsuchmodelsislargelydependontheapplicationsofthewallfunctionsthatrelatesurfaceboundaryconditionstouidpointsawayfromtheboundaries.Therebyweavoidtheproblemofdirectmodellingtheviscosityinuence.HereweusethewallfunctiontreatmentbyLaunderandSpalding[ 52 ]asoneexampletoillustratetheclosureproblem.Thegoverningequationsforstandardk-"modelscanexpressedin


theformshownbelow, @xi(t @xi)+2tSijSij";(2.20) @xi(t @xi)+2C"1t" kSijSijC"2"2 Theeddyviscositytisdenedas whereCisadimensionlessconstant.Sijisthemeanrateofstrain. Theseequationshaveveadjustableparameters,C,k,",C"1,andC"2.Here,forbrevity,weonlygivethecoecientsforstandardk-"modelsuggestedbyLaunderandSpalding[ 52 ].Theseparametersarevalidforawideofrangeofturbulentows. Toclosethesystemofequationsforturbulencecomputation,weneedanadditionalBoussinesqrelationshipthatrelatestheReynoldsstresstothemeanrateofstrainandturbulencekineticenergy, uiuj=t(@Ui 3kij=2tSij2 3kij:(2.24) 2.2.2 BoundaryConditionsforTurbulenceModels Theellipticcharacteristicofthekand"modelequationsgivesriseofbound-aryconditions.Usuallyweassumezerogradientattheoutletboundary.Atinletboundary,ifthereisnomeasurementofkand"atdisposal,onewayistomakeanapproximationbasedonthefreestreamvelocityU1,turbulenceintensityTi,andthecharacteristiclengthscaleL.Inourstudyweadoptthefollowingformulastodene


thefreestreamboundaryconditionsforkand",k=3 2(U1Ti)2; TheturbulenceintensityTiisusuallysetbetweentherangeof0.02and0.05.AnotherwayistoadjusttheturbulenceReynoldsnumberRet=U1L=t.ForpracticeweusuallysettheturbulenceReynoldsnumberintherangeof100and500. Thesolidwallboundaryconditionsaremorecomplicated.Nearthewallbound-aries,thelocalReynoldsnumbersarelowandtheviscouseectsareimportant.However,thek-"turbulencemodelthatwementionedaboveisbasedonthehigh-Reynolds-numberassumptioninwhichtheturbulentstressforcesaredominated.Moreover,closetothewall,extremelynegridisusuallyrequiredtocapturethesteepvariationsinowprolestosolvetheNavier-Stokesequationsdirectly.TheCPUtimebasedonsuchanegridisalmostbeyondtheabilityofcurrentfacilities.Inordertosolvethesediculties,twoapproachesarecommonlyemployed,oneisthewallfunctiontreatment[ 52 ],theotheristhelow-Reynolds-numberversionsofthek-"[ 46 73 ]. Thewallfunctiontreatmentispopularinpracticebecauseofitssimplicityandcapabilityforcomplexgeometries.UnderthehighReynoldsnumbercondition,owcanbelargelydividedintotwolayers:theouterlayerandthesurfacelayer.IntheouterlayertheReynoldsstressisdominated;whileinthesurfacelayerboththeReynoldsstressandtheshearstressareimportant.Intheouterlayer,theuideldisinuencedbytheretardingeectofwallthroughthevalueofthewallshearstress,butnotbytheviscosityitself. Inthewallfunctiontreatment,velocityproleshavedierentdescriptionsintheoutandsurfacelayers.Thevelocityintheoutlayerisgovernedbythevelocity-defect


lawexpressedas u=g(y );(2.27) whereistheboundarylayerthickness,Umax-Uisthevelocitydecit,gisanunknownfunctionwiththeorderof1.Thevelocityinthesurfacelayerisgivenbythelawofthewall u=f(y+):(2.28) Herey+isanon-dimensionalnormaldistancefromthewall,anditisrepresentedby :(2.29) Intheaboveexpressionyisthedistancefromthesolidwall,anduisthefrictionvelocityatthewalldenedby Theshearstressatthewall,w,isgivenby @y:(2.31) Thelow-Reynolds-numberk-"modelrequiresarelativelynegridresolutionnearthewallregions.FordetaileddescriptionswerefertoRefs.[ 46 73 87 ]. Tomakeappropriateuseofthetwoversionsofturbulencemodels,westrivetoadjustthevalueofy+oftherstpointnexttothewalltosatisfycertainrequirements.Forthelow-Reynolds-numberapproach,y+shouldbelessthan2foranaccurateresult,whiley+intherangebetween30and500[ 103 ]canproduceaccurateresultsforthewallfunctiontreatment.Anequationfortheatplateboundarylayertheoryisusedinthegridgenerationprocesstogetanestimationofpropercellsizefory+,


Italsoshouldberemember,whilethepositionoftherstcelliscritical,itmusthaveenoughcellsinsidetheboundarylayertoresolvethesteepvariationsofowvariables. ThesuccessoftherecentturbulenceclosureisrestrictedtothesucientlylargeenoughReynoldsnumberturbulenceows.Furtherrenementisneededbeforetheyareusedwithcondencetocalculatenear-wallandlowReynoldsnumberows.Amoredetailedreviewaboutturbulencemodelsfornear-wallandlowReynoldsnumberowsisgivenbyPateletal.[ 73 ]. 2.3 Laminar-turbulentTransition Areliablecomputationoftheboundarylayertransitionremainsoneofthemostchallengingdemandsinaerodynamics.Sincethetransitionalseparationbubblesandtheirlossesplayanimportantroleforthepressureanddragdistributions,anac-curaterepresentationofbothlaminarandturbulentseparatedowsiscriticalwhendragpredictionneedstobeaccuratelypredicted.However,thefundamentaluiddy-namicsproblemoftransitiontoturbulencehasresistedanysimplephenomenologicaldescription.Thereisstillnocomprehensivetheoryoftransition. Therangeofexistingtransitionpredictionmethodsextendsfromsimpleempiri-calrelationshipsthroughthosebasedonparallelandlinearstabilitytheories,liketheeNmethods[ 15 75 ],orlinearornonlinearparabolizedstabilityequations(PSE)[ 31 ],orthelow-dimensionalmodels[ 76 ],todirectnumericalsimulation,likelarge-eddysimulationtechniques.


Forproblemsinvolvinginteractionsbetweenuidowsandmovingstructures,thegeometryofthesolidobjectisnotknownapriori.Inthecourseofcomputationitisnecessarytotrackthegeometry,theeldequations,andtheinterfacialconditiontoensurethatallrequirementsaremet.Tofacilitatethesolutionofsuchmovingbound-aryproblems,themovinggridtechnique(ordynamicgridtechnique)isemployedtoadjustthecomputationalgriddynamicallyalongwiththegeometricupdates.Itisdesirablebutdiculttodevelopanautomatedregriddingproceduretoensurethatthenewgridnotonlymatchesthegeometricchangesbutalsomaintainssatisfactorycharacteristicssuchasnon-overlapping,smooth,andnotexcessivelyskewed.Themovinggridtechniquecanalsobeusedbyshapeoptimization[ 54 ]. Formovinggridproblems,severalalternativeapproacheshavebeenproposed.Schusteretal.[ 81 ]usedanalgebraicshearingprocessintheirstaticaeroelasticanal-ysisofaghteraircraft.Thebasicideaistoassumethatthedisplacementofthesubjectisredistributedalongthegridline,whichconnectstheaeroelasticsurfacetotheouterboundary.Thissimplemethodworksquitewellforsingleblockgridwithmodestdeection;however,whenappliedtomultiblockgridarrangements,thisapproachusuallyneedsextensiveinterventions.Thespringanalogymethodwasrstproposedforunstructuredgrids[ 3 ],andlaterwasextendedtostructuredgridsbyRobinsonetal.[ 78 ].Thismethodregardsallgridsasconnectedbysprings,andeachspringhasthestinessinverselyproportionaltothespacingbetweentargetvertexanditsneighbors.Comparedwithothercurrentlyusedgridregenerationmethods,thespringanalogymethodneedsmorememoryandCPUtime.Directtransnite 19


interpolationbyEriksson[ 18 ]isapopularmethodbecauseofitseciencyforstruc-turedgrid.HartwichandAgrawal[ 29 ]alsoproposedtheMaster/Slaveconcepttoexpeditethegridregenerationprocessandminimizetheuserintervention. 3.1 PerturbationMethod Forsmallamplitudedisplacementononefaceofthecomputationalblock,asimplebutecientone-dimensionalperturbationmethod[ 56 77 ]iscanecientlyregeneratethegrid.Theone-dimensionalperturbationmethodobeysthefollowingformula, wherexirepresentsthelocationofanarbitraryinteriorgridpoint,xsrepresentsthelocationofagridpointonthemovingboundary,andSrepresentsthenormalizedarclengthalongtheradialmeshlinemeasuredfromtheouterdomain,whichisgivenby Tousethisone-dimensionalperturbationmethod,oneneedstoknowthedisplace-mentsofthetwoendingpointsofameshline.ThepositionsoftheinterimpointsaresolelydeterminedbythedisplacementsoftheseendingpointsandtheiroriginalpositionsasshowninEq.( 3.1 ).Thesedisplacementsareconsideredasperturbationsources.Intuitively,theperturbationwillspreadthroughoutthewholedomainwhileexertingmoreeectsonnearerpoints.BasedonEq.( 3.1 ),amorecomplicated3-stageperturbationmethod,analogoustothetransniteinterpolation(TFI)method[ 18 ],wasproposedforthree-dimensionalproblemstreatedbythesingle-blockapproach.UnliketheoriginalTFImethod,thisperturbationmethodconsiderstheoriginalgriddistribution,andpreservestheoriginalgridquality. Takeatwo-dimensionalcaseasanexample,inFig. 3{1 ,arectangleundergoesshapechange.Toregeneratethegridforthechangedgeometry,theperturbationmethodproceedsasa2-stageprocess.Therststageistomatchcornerpoints,i.e.,


tomovethecornerpointAtoitsnewpositionA0,andBtoB0,etc.(SeeFig. 3{1 b).Oncethecornerpointsmatchtheirnewpositions,interimedgesaregeneratedbasedonEq.( 3.1 ).ToadjustthepositionofaninternalpointO,eectsfromthefourcornerpointsneedtobeaccountedforinacoordinatedmanner.Inthisapproach,P0,Q0,S0,andT0,whichconnectO0toitsboundary,willberstadjusted.Thedisplacementsofthefouredgepointsin,e.g.,thex-direction(P,Q,S,T)canbecalculatedwithEq.( 3.1 ),andthemovementofO0isevaluatedasacombinationofthesedisplacements.Specically,onecanadoptthefollowingexpression: O0=abs(I)I+abs(J)J where,I=(1Si)P+SiQ,J=(1Sj)S+SjT,SiandSjarethenormalizedarclengthsalongmeshlinesPQandSTrespectively.Displacementsinotherdirectionscanbecomputedinthesimilarway. Aftertherststage,thefourcornerpointshavematchedtheirnewlocations;however,theinterimedgesdonotnecessarilymatchtheirnalpositions.Thesecondstageistomatchtheinterimedgestotheirnaledgelocations.Thedierencebe-tweentheinterimandnaledgelocationsaretakentobethesourceofperturbation.Inthisstage,sincetheperturbationsinbothdirectionsareindependent,thecontri-butionsfromallindividualdirectionsareadded[SeeFig. 3{1 c].Forthree-dimensionalproblems,a3-stageprocessperturbationmethodwasproposedinRefs.[ 56 77 ]. 3.2 Master/SlaveConcepts Theaforementionedperturbationmethodcanbeimplementedecientlyforsin-gleblockandsmalldisplacementproblems.However,whenamulti-blockgridisinvolved,thismethodneedstobeenhancedbyothertechniques.TheMaster/Slaveconceptactsasausefulapproachtomaintaingridqualityandtopreventpotentialgridcross-over.


Figure3{1: Twostageperturbationmethod:(a)ABCDchangestoA0B0C0D0;(b)rststage:matchingcornerpoints;(c)secondstage:matchingedgepoints.


Althoughunstructuredgridhasbeenusedtosolveuiddynamicsproblemsbydierentresearchers[ 3 26 ],here,ourattentionislimitedtostructuredgridonly.Foramulti-blockstructuredgrid,CFDsoftwareoftenrequirespoint-matchedgridblockinterfaces.Abasicrequirementforgridregenerationinmovingboundaryproblemsistomaintainpoint-matchedinterfaces.Whilethisupdatingprocessiseasilysatisedbythoseblocksurfacesthatcoincidewiththebodysurface,theremainingblocksurfacesneedtobehandledwithduecare.Tousetheperturbationmethod,oneneedstospecifythenewpositionsoftheo-bodysurfacepoints,or,atleast,theverticesofsuchblocksurfaces. Whenanobjectchangesitsshape/position,masterpoints,whicharedenedasthepointslocatedatthebodysurface,rstmove,andthenaectdistributionsoftheslavepoints,theo-bodypoints.Whilethetaskofredistributinginteriorgridscanbechallenging,inpractice,itismorediculttomovetheverticesofeachblockifapoint-to-pointmatchedinterfacewithoutoverlapisrequired.Iftwoabuttingblockshavecongruentinterfaces,suchastheinterfaceofblock2andblock3inFig. 3{2 a,oncethecornerverticesaredetermined,theedgepointsandtheinteriorpointscanbeuniquelysettledbasedonEq.( 3.1 ).However,whentwoneighboringgridblocksdonotshareidenticalinterfaces,suchasblock1andblock2showninFig. 3{2 a,thisprocedurecancausediscontinuityattheinterfacebecausetheregridingproceduresattheinterfacearebasedondierentstencilsfordierentblocks.Slavingverticestotheirmasterpointscanavoidcreatingundesirablegriddiscontinuities. TheMaster/Slavecouplingisbasedonthedistancebetweenano-bodypoint(Slavepoint)andasurfacepoint(Masterpoint).Thedistanceisgivenby: wheresubscriptvrepresentsano-bodyvertex,and,b,abodysurfacepoint.Foreacho-bodyvertex,thenearestbodysurfacepointisdenedasitsmasterpoint.


Figure3{2: Testsonthemovinggridscheme:(a)initialgeometry;(b)sideviewofthegrid;(c)deformedshape;(d)newgriddistribution;(e)griddistributionwithhighstinessparameter.


Tosimplifytheconnectionbetweenthemasterandslavepoints,onecanconvertthethree-dimensionalparameterspace(i,j,k)intoaone-dimensionaldatastructure.Eachgridpointhasanidenticationnumberldenedby wherebo(bn)istheidenticationnumberoftherstpointinblockbnintheone-dimensionalarray,imaxandjmaxarethemaximaldimensionsintheiandjdirections,respectively.Inthiswayaslavepointisassociatedwithitsmasterpointbystoringitsmaster'sidenticationnumberintoitsaddress. Oncetherelationisestablished,thenextstepistodeterminehowtheslavepointmovesaccordingtotheinuencefromitsmasterpoint.Intuitively,anearervertexwillhavemoreeect(displacement)thanthosefaraway.Also,toavoidalowersurfacetocrossoveritsoppositefacewhenthemovementisrelativelylarge,oneneedstoscalethemovement.Asimplebuteectiveformulaisexpressedasfollows: ~xs=xs+(~xmxm);(3.6) wheresubscriptsmandsrepresentthemasterandslavepoints,respectively,thetilde()indicatesthenewposition,andisadecayfunction.Inourcomputations,aGaussiandistributionfunctionisused whereisaparameterwhichaectsthestiness.


3.3 NumericalTests Themovinggridtechniqueconsistsoftwoelements,oneistheperturbationmethod,andtheotheristhemaster/slaveconcept.Theproposedmethodistestedonathree-dimensional,multi-block,andstructuredcomputationalgrid.TheinitialcongurationisshowninFig. 3{2 a,block1sharestheinterfacewithblocks2and3.TheinitialgridisshowninFig. 3{2 b,ithasclustereddistributionsatthetopandbottomsurfaces,anditalsohaspoint-to-pointmatchedinterfaces.Thebottomtwofacesarethemovingboundariesthatexperiencelargedeformation.Thecoordinatesofeachvertexarelabelled.ComparingFig. 3{2 aandFig. 3{2 cwecanseethedisplacementisaslargeashalfofthetotalheight.Thenewcomputationalgridiscomputedwiththemethoddiscussedinthischapter.ThenewgridisshowninFig. 3{2 d.Thatguredoesnotshowanycross-overinthenewgridafterthelargedeformation.Thenewgridnotonlykeepsthepoint-to-pointmatchedinterfaces,butmaintainstheclustereddistributionsatthetopandbottomsurfaces.TheCPUtimeforregeneratingthecomputationalgridislessthan10%oftheCPUtimeusedbytheuidsolver. Theparametercontrolsthestiness,wend=1/64workswellforourprob-lems.Alargercausestheblocktobehavemorelikearigidbody(SeeFig. 3{2 e).MoredetailsaboutthemasterandslaveconceptcanbefoundinRefs.[ 29 56 ].


4.1 Introduction Themainobjectivesforuiddynamistsandaeronauticalengineersthroughoutthehistoryofmodernighthavealwaysbeentoincreaseliftandreducedrag.Microairvehiclesareairbornevehicleswithcharacteristiclengthsunder6inchestravelingatspeedsaround30ft/s.Theyaredesignedtoperformanarrayofdierenttaskssuchasreconnaissance,targetdesignationorbordersurveillance.Toenhancetheperformanceofthesesmallvehicles,improvementsshouldbemadeintherespectsofpowerplants,thrustgenerators,navigationsystems,cameras,datatransmitting,andmoreimportantly,theiraerodynamicperformance.Asfarastheaerodynamicsisconcered,MAVshavetooperateinoneoftheworstightregimesimaginable,oneofthemisthelow-Reynoldsnumberswhichisusuallynorotiouslyaccompaniedwithlaminarboundaryseparation,transition,andlowlift-to-dragratio.FurthermorethesmalldimensionsoftheMAVsgreatlyenhancesthetradeodicultiesmentionedabove.Thesmallsizeinhibitsthedesigner'sfreedomtoachievethedesiredaerody-namicperformance.Thelimitationsinthepowersourcesavailable,atpresenttime,makethisequationevenmorecomplicated. Tosolvetheaerodynamicchallengesoftheseproblems,twomajoreldshavebeenevaluated.Therstapproachistryingtocopynature'sownway;thatofappingwings[ 43 45 59 83 ].Thesecondapproachisusingadynamicalwing.Therearesomedierentimplementationsinthiseld.Onewidelystudiedmethodistheuseofanwingwithaexiblemembraneuppersurface[ 53 55 83 92 ].Anothermethodusesactivedevicestocontroltheowontheuppersurfaceofthewing.Thelattermethodhasgainedmoreattentionwiththesuccessmicroelectromechanicalsystems(MEMS). 27


Inactiveowcontrolarea,apiezo-actuatedapmountedonalowReynoldsnumberwinghavebeenstudiedbothexperimentally[ 19 ]andcomputationally[ 30 56 ].Gad-el-Hak[ 20 ]hasdiscussedthepossibilityofemployingMEMStodelayorpreventseparationonMAVs. TounderstandthelowReynoldsnumberwingcharacteristicsandimprovethewingperformance,wewillperformeacomputationalinvestigationforowssurround-ingadynamicallyshapedwing,atachordReynoldsnumberof78,800,alongwithaparallelexperimentaleort.Weadoptapiezo-actuatedapontheuppersurfaceofaxedwingforactivecontrol.Theactuationfrequencytheyfocusonis500Hz.Theircomputationalframeworkconsistsofamulti-block,movinggridtechniquetohandlethegeometricvariationsintime,usingtwo-equationturbulenceclosures,andapressure-basedowsolver.Thedynamicgridredistributionemploysthetransniteinterpolationschemewithaspringnetworkapproach.Comparingtheexperimentalandcomputationalresultsforpressureandvelocityelds,eectsofthedetailedge-ometriccongurationoftheap,theappingamplitude,turbulencemodeling,andgriddistributionsontheowstructurewillbeassessedinthischapter. 4.2 ExperimentalandComputationalSet-up ThewinghasmovingapsontheuppersurfaceasshowninFig. 4{1 .Theapcanbeactuatedatprescribedfrequencies[ 19 ].Thewingstudiedintheexperimenthasaspanof1ftandachordof5:400.Intheexperimentalset-up,thereisagapof0:01600widebetweentwoaps,eachis200wideand0:56400long.TheReynoldsnumber,basedonthefreestreamvelocityU=27.5ft/s,is78,800.Theowisinitiatedaslaminar;afterthetransitionpoint,theturbulencemodelisinvoked.Thetransitionpointisidentiedtobeatthetipoftheap.Inthecomputation,theturbulentReynoldsnumberbasedonthechordlength,Ret=UL=t,attheinletoftheturbulentowregionissetto500.Theamplitudeandshapeoftheapareprescribedbyanalyzing


theexperimentaldataobtainedwithPIV[ 19 ]Forthe500Hzcase,twoappingamplitudesareconsidered,namely0:00100and0:007300. Figure4{1: Movingapwingexperimentalsetupinthewind-tunnel Inthischapter,anumberoftwo-dimensionalandthree-dimensionalcases,in-cludingdierentgriddensities,withandwithoutthegap,willbestudied.For3-Dcases,aneorthasbeenmadetoconsidertheeectofthegapbetweenapsasshowninFig. 4{2 .TheoriginalLaunder-Spaldingk-"model[ 52 ]withthewallfunction,aswellasalow-Reynoldsnumberversionk-"model[ 46 ]arecompared. Figure4{2: Flapwingwithrenedthreedimensionalbendingmodel. 4.3 NumericalStudy 4.3.1 GridSensitivityAnalysis Anumberof2-Dgridswithdierentdistributionsandsizeshavebeentested.Thenareference2-Dgridisselectedbyensuringthatsatisfactorysolutionqualityisobtainedonamosteconomicgridsize.Then,the3-Dgridisestablishedbyextendingthereference2-Dgridalongthespanwisedirection.Theresultinggrid,forcaseswithgap,has33nodesinthespanwisedirection,andconsistsof20blockswithapproximately533,000gridpoints.Forcaseswithoutgap,thegridsystemhas


20nodesinthespanwisedirection,andconsistsof8blockswithabout253,000gridpoints.Withthemulti-blockstrategyitisconvenienttoemploythetransitionmodel,i.e.,whenatransitionlocationisidentied,newcomputationalblocksaregeneratedsurroundingthetransitionpoint,sothatthelaminarowequationsaresolvedintheblocksbeforethetransitionlocation,andtheReynolds-averagedturbulentowequationsaresolvedintheblocksafterthetransitionpoint.Inourcomputations,thetransitionpointisidentiedfromtheexperiments[ 19 ],anditisatthetipoftheap,therefore,ourcomputationalblocksaresetupbasedonthistransitionpoint. 4.3.2 AssessmentoftheTurbulenceModel IntheearlyworkbyHeetal.[ 30 ],thewallfunctionapproachisusedexclusively.ForsuchalowReynoldsnumber,thelow-Reynoldsnumberversionofthek-"modelseemsmoresuitable[ 56 ].Fig. 4{3 comparesthetime-averagedpressuredistributionsat500Hzbetweentwodierentturbulencemodelsforthecasewithoutgapbetweenaps.Bothexperimentalandcomputationaldataarecollectedatthemiddleoftheapandalongthechordwisedirection.Asfarasthecomputationaltimeisconcerned,itisobservedthatthelow-ReynoldsnumbermodeluseslessCPUtimethanthatofthewallfunctiontreatment. 4.3.3 EectsofGapsbetweenFlaps Theinclusionofagapobviouslyincreasesthecomplexityoftheoweldinthespanwisedirection.ItcanbeseenfromFig. 4{4 thatstreamlinesenterthegapandmovetowardthemiddleoftheap.Aweakrecirculataryowwithavortexstretchingfromthewingsurfacetotheapispresentonbothsidesofthegap,rotatingintheclockwisedirectionontherightandcounter-clockwisedirectionontheleftofthegapwhenviewedfromabove. Thegapalsoaectstheoverallpressuredistributiononthewing.Fig. 4{5 comparesthepressuredistributionsobtainedwithandwithoutgapbetweenaps.Thelow-Reynoldsnumberturbulencemodelisusedinbothcases.Timeaverageis


Figure4{3: Comparisonbetweentime-averagedpressuredistributionsforcaseswithoutgap,at500Hz,obtainedwithdierentturbulencemodels.(Timeaveragedfortheentirecycle,amplitude:0:00100). Figure4{4: Streamlinefor3-Dcasewithgapat500Hz,amplitudeis0:00100


madefortheentirecycleforanamplitudeof0:00100andanactuationfrequencyof500Hz.Sincethepressuretapsarelocatedatthemiddleofaap,theeectofthegapmaynotbeclearlydemonstratedbytheguresshown.Fig. 4{6 comparesthevelocityprolesimmediatelydownstreamfromtheaptrailingedge.Thelow-Reynoldsnumberturbulencemodelisusedagainandthecomparisonismadeat0:016200downstreamfromthetrailingedgeoftheap.Casesinvolvingbothstationaryapsandactuatedaps,at500Hz,areexamined.Consistentwiththatalreadyobserved,thevelocityprolesindicatethattheowenterintothegapfromtheuppersurfaceoftheap,maintainingthesamestreamwisedirectionwhilemovinglaterallytowardtheregionbetweentheapandthewingsurface.Arecirculatingowpatternisformedbetweentheapandthewingsurface. Figure4{5: Comparisonoftime-averagedpressuredistributionsforcaseswithandwithoutgapbetweenaps. 4.3.4 EectsofFlappingAmplitude Dependingontheappingamplitude,thereattachmentpointcanmovebackandforth,andtherecirculatingowcaneitherremainattachedtothetipoftheap,orbeshedaway.Inourcomputationthecurvatureoftheapissetbycollecting


Comparisonofvelocityprolesattwodierentspanwisepositionswithgapimplemented:(a)steadycomputation;(b)f=500Hz. theexperimentaldataandthenapproximatedbyasimpliedmodel.Forthe500Hzcaseinourcomputationtwoamplitudesareconsidered,oneis0:00100andtheotheris0:007300.Forthesmalleramplitudecase,themotionofthetrailingedgeismodeledasthefollowing, wherezisthespanwisedistancefromthestudypointtothegap,fisfrequency,andtistime.Thelargeramplitudecasesimplyscalestheaboveequationbasedonthesameexpression. Fig. 4{7 comparesthetime-averagedpressurecoecientsbetweenthetwoap-pingamplitudes.Forthehigherappingamplitude,thetime-averagedreattachmentpointmovesclosertothetrailingedgeoftheap,becausealargermovementoftheapintroducesmorehighmomentumuidsintotheboundarylayer.Intheexperi-mentweobservevortexsheddingasshowninFig. 4{8 a.Inourcomputationswiththesmallerappingamplitude,thevortexcoreisanchoredattheaptip,exhibitingtwoco-rotatingcoresastheapsmoveupwards,asshowninFig. 4{8 b.Ontheotherhand,withthehigheramplitude,vortexsheddingappearsforcaseswithorwithout


Figure4{7: Pressurecoecientatdierentactuationamplitudewithoutgapbetweenaps.Actuationfrequency:500Hz. thegap.Fig. 4{9 showsselectedinstantaneousplotsofthevortexstructureduringasinglecycleoftheapmovement.Asexpected,therecirculationisstrongerwhentheapmovesupwardandbecomesweakerwhentheapmovesdownward. Instantaneoustrailingedgestreamlinesforthefrequencyof500Hzasapismovingupwardsforamplitudeof0:00100.Left:experimentalresult;right:computationalresult. 4.4 Conclusions Baseonourresultswendthatthelow-Reynoldsnumberformulationofthek-"turbulencemodelproduceslongercirculationzone.Thegapbetweentheapsaectsthelocationofthere-attachmentpoint,andcreatesathree-dimensionalowpattern.Aidedbythelimitedexperimentalinformation,thecomputationalmodelcanoer


Trailingedgestreamlinesforthefrequencyof500Hzinthethirdcyclewithanamplitudeof0:007300 Astothephysicalmechanisms,whileitisexpectedthatthevortexsheddingappearswhentheamplitudeishigh,itisnotclearwhyandwhenthisphenomenonceasestoexistwhentheamplitudeisreduced.Similarly,whileitisobviousthatthegapbetweenapsaectstheowreattachment,itishardtopredictthetrendquantitatively.


5.1 Introduction Thedeformationsofmembranestructuresareoftenlargeandthusinherentlynonlinear.Qualitativedescriptionsofthethinmembranebehaviorsareoftenmorecomplexthanthatofthree-dimensionalbodies.Untilrecently,membraneanalyseswereprimarilyreliedontrialanderror.GreenandAdkins[ 25 ]arethersttoderiveageneralnon-lineartheoryformembranes.Themembranetheoryhastwofundamentalassumptionsbecauseoftheirthinness[ 41 ].First,themembranetensionstinessismuchlargerthanitsbendingstiness,therefore,thestresscoupleeectscanoftenbeneglected.Second,theratioofthemembranethicknesstoitsradiusofthecurvatureissmall,whicheectivelydecouplesthestrain-displacementrelationfromthecurvaturetensor.Practicallyspeaking,thereexisttwoapproachestomodeltheresponseofamembranestructure.Therstapproachistoidealizeaplateorshellstructureasamembraneawayfromitsboundaries,inwhichthestresscouplesareneglectedintheinteriorregion.Byformallyequatingtheplatestinesstozero,thismembranemodelcanbedescribedbytheFoppl-vonKarmantheory[ 40 41 ].Inthesecondapproach,thestructurecannotsustainstresscouplesanywhere.ThesecondapproachisadoptedinthecurrentworkbecauseitisappropriateforMAVcomputations. Two-dimensionalmembranewingstudywithxedleadingandtrailingedgeswereoriginallyconductedbyNelsen[ 68 ],Thwaites[ 101 ],Voelz[ 108 ].Theseworksconsideredinextensiblemembranessurroundingbysteady,two-dimensional,andir-rotationalows.Asaconsequenceoftheinextensibleassumption,whencambersandincidenceanglesaresmall,problemscanbelinearizedandexpressedcompactlyina 36


integralformwhichisreferredasthe"sailequation"byNewman[ 69 ].Thefollowingequationcompletelydenesthelinearizedtheoryofinextensiblemembranewingsinsteadyandinviscidowelds. 1CT wherey(x)denesthemembraneproleasafunctionofthexcoordinate,istheincidenceangle,andCTisthetensioncoecient.BasedontheworkofVoelz[ 108 ],thesailequationhasbeensuccessfullysolvedanalyticallyandnumerically.Newman[ 69 ]gaveacomprehensivereviewoftheearlierworksrelatedtomembranewingaerodynamics. Thedevelopmentofathree-dimensionalmembranemodelintroducesseveralcomplicatingfactorsnotencounteredintwo-dimensionalanalyses.First,unlikethe\sailequation",thetensioninathree-dimensionalmembraneisnolongerasingleconstantvaluebutisatbestdescribedbyastateofbiaxialtensionalongthelinesofprincipalstresses[ 38 ].Second,geometricandmaterialpropertiesvaryalongthespanwisedirection,andtheyhavetobedescribedaccurately.Third,ifoneoftheprincipaltensionsvanishes,themembranewillbecompressedandwrinkled,whichfurthercomplicatestheanalysis.Consideringthecomplexityofthethree-dimensionalproblem,simplicationsaremadeintheearlyworks. OdenandSato[ 70 ]presentedaniteelementmethodfortheanalysisoflargedisplacementsandnitestrainsinelasticmembranesofgeneralshape.Intheirworkthestaticequilibriumofaninatedmembranewasconsidered.Astaticanalysiscanbeconsideredasaspecialcaseofadynamicanalysisinwhichtheinertiaeectisincluded.Worksonthedynamicanalysisofmembranesbefore1991werereviewedbyJenkinsandLeonard[ 41 ];mostrecentworkswerereviewedbyJenkins[ 40 ].Verronetal.[ 105 ]performedacombinednumericalandexperimentalstudyofthedynamicinationofarubber-likemembrane.Membranewrinkleproblemsweretackledin


theworkofDingetal.[ 14 ];theyperformedanumericalstudyofpartialwrinkledmembranes.ThestabilityissuesofauidowpastamembranewerereportedbyThaokarandKumaran[ 99 ]. Inthischapterathree-dimensionalmembranemodelisproposedforthemem-branedynamics.ThemembranematerialconsideredobeysthehyperelasticMooney-Rivlinmodel[ 64 ].WeusetheGreen-Lagrangestraintensorforthedescriptionoflargestrains.Thedynamicresponseofthemembraneisdescribedbyasystemofsecond-orderordinarydierentialequations.Aniteelementmethodwithtriangu-larelementsforspatialdiscretizationisutilized.Fortimeintegration,theimplicitWilson-methodisused. 5.2 StructuralSolver 5.2.1 Mooney-RivlinModel Foraninitiallyisotropicmembrane,GreenandAdkins[ 25 ]showedthatthereexistedastrainenergyfunctionWwhichwasexpressedasthefollowingform whereI1,I2,andI3aretherst,second,andthirdinvariantsoftheGreendeformationtensorC.Ifthematerialisincompressible,namely,I3=1,thenthestrainenergyisafunctionofI1andI2only.Itcanbewrittenas: wherec1andc2aretwomaterialparameters.AmaterialthatobeysEq.( 5.3 )isknownastheMooney-Rivlinmaterial[ 64 ],andthemodeliscalledtheMooney-Rivlinmodel.TheMooney-Rivlinmodelisoneofthemostfrequentlyemployedhyperelasticmodelsbecauseofitsmathematicalsimplicityandrelativelygoodaccuracyforreasonablylargestrains(lessthan150percent)[ 64 ].


Foraninitiallyisotropicmembrane,thegeneralstress-strainrelationshipiswrit-tenas @C;(5.4) whereSisthesecondPiola-Kirchhostresstensor,pisthehydrostaticpressure.Ifthemembraneisincompressible,thenthestress-strainrelationcanbefurthersimpliedasthefollowing whereIis3-by-3identitymatrix.Inthemembranetheorythestresseldisessentiallytreatedastwo-dimensional,andthestressnormaltothedeformedmembranesurfaceisnegligibleincomparisonwiththestressesinthetangentplane[ 70 ].Underthisconsideration,thedeformationandstressmatrices,C(t)andS(t),aregivenby Therefore,thehydrostaticpressureisdeterminedbyconsideringtheconditionthatS33=0, where3isdenedby inwhichh(t)andh0arethemembranethicknessinthedeformedandundeformedcongurations,respectively. 5.2.2 PrincipleofVirtualWork Consideranelasticmembrane,whichinitsinitialstateoccupiesaniteregioninspace.ThepositionoftheregioncanbespeciedbyanorthogonalcoordinatesystemXYZwhichisreferredtoasglobalcoordinates.Inglobalcoordinates,the


generalformoftheprincipleofvirtualworkintheLagrangiandescriptionisgivenby whereV0and@V0are,respectively,thevolumeandtheboundarysurfaceoftheundeformedmembrane,0isthemembranedensity,D(t)isthedisplacementvectoringlobalcoordinates,D(t)isthevirtualofvectorD,thesuperscriptTreferstothetransverseofavectororatensor,D(t)istheaccelerationvector,P(t)isthegeneralizedexternalsurfaceforce,EistheGreen-Lagrangestraintensor,SisthesecondPiola-Kirchhostresstensor,andE:Sisthescalarproductoftwotensors.IntherectangularCartesiancoordinates,thescalarproductisdenedby withsummationontherepeatedindices. 5.2.3 TriangularMembraneElement Theprincipleofvirtualworkcanbewritteninthefollowingdiscretizedform whereisthetotalvirtualwork,andthesubscripterepresentselementaryquan-tities,Neisthenumberofelementscoveringtheentiremembrane.Inthepresentworkthetriangularelementisemployed.AsshowninFig. 5{1 ,atypicaltriangularelement,e,isdenedbynodes,i,j,k,andthreestraight-lineboundaries.Intheglobalcoordinatesthethreenodeshavecoordinates(Xi,Yi,Zi),(Xj,Yj,Zj),and(Xk,Yk,Zk),respectively. Tofacilitatethedevelopmentofacomputationalprocedure,itisconvenienttoconsiderthemembraneelementanditsdisplacementinthelocalcoordinates.AsshowninFig. 5{1 ,inlocalcoordinates,thetriangularelementliesintheplaneof


Schematicglobalandlocalcoordinates. whereXiistheglobalcoordinateofnodei,andTisanorthogonaltransformationmatrixwhichsatisesTT=T1.Similarlythedisplacementshavethefollowingrelation whereDisthevalueofadisplacementinglobalcoordinates,disitsvalueinlocalcoordinates. Iftheelementissucientlysmall,wecanapproximatethedisplacementatanypointwithintheelementwithalinearcombinationofthedisplacementsofnodes,i,j,andk,


whereN3isaninterpolationmatrixdenedas where01,01,and011.deisavectordenedby inwhichdi,dj,anddkarethenodaldisplacementvectorsofnodesi,j,andk,respectively.Eachnodehasthreedegreesoffreedom, 5.2.4 Green-LagrangeStrainTensor Whenadisplacementisintroducedthelocationofapointoriginallygivenbycoordinatesxnowhasnewcoordinates whereuisthedisplacementasshowninFig. 5{1 .FollowingGreenandAdkins[ 25 ],theGreen-Lagrangestraintensorisdenedbytherelation 2@yk Sincethethicknessofthemembraneissmallcomparedtotheothercharacteristicdimensions,thefollowingsimplicationsareallowedtobemade[ 25 70 ] 2@u 2(21);;=1;2;(5.20)


whereisthesameas3inEq.( 5.8 ).Tosimplifythecomputation,wedividetheGreen-Lagrangestraintensorintotwoparts.OnepartisE0,whichtakesintoaccountsmallstrains, @x1 2(@u @y+@v @x)01 2(@u @y+@v @x)@v @y00003777775;(5.21) andtheotherisEL,whichisthenonlinearpartofthestraintensor, @x@u @x+@v @x@v @x+@w @x@w @x@u @x@u @y+@v @x@v @y+@w @x@w @y0@u @x@u @y+@v @x@v @y+@w @x@w @y@u @y@u @y+@v @y@v @y+@w @y@w @y000213777775:(5.22) TondeachcomponentinthestraintensorsE0andEL,itisnecessarytocomputethedisplacementgradientur(x;y).Aftersomeoperationsweobtaintheexpressionofdisplacementgradient, @x@u @y@v @x@v @y@w @x@w @y3777775=26666664uiuj BasedontheexpressionsinEq.( 5.23 )wecancomputethecomponentsE0andEL. 5.2.5 InternalForces Whenlargedeformationsand/ornonlinearmaterialpropertiesareinvolved,theequivalentinternalelasticresistingforcesofthestructurecanbeexpressedasthefollowing, whereisthenonlinearstressvector,"isthestraintensor,andBisthestrain-displacementmatrixwhichisalsoafunctionofthedisplacementforlargedeforma-tion.Recallthatintheprincipleofvirtualworkthetermrelatedtotheinternal


forcesfromthediscretizedformofprincipleofvirtualworkEq.( 5.11 )reads Eq.( 5.6 )tellsusthatthemembranestressisessentiallytwo-dimensional,therefore,weknowfromEq.( 5.10 )thatitnecessarytoconsiderthereducedstraincomponentsinsteadofallthecomponentswhencalculatingtheinternalforcesbasedonEq.( 5.24 ).ThegeneralmembranestrainvectorcanbeconvenientlydecomposedintotwopartsaccordingtothestraindecompositioninEqs.( 5.21 )and( 5.22 ), where"0isthereducedstrainvectorforthelinearpartofthestraintensor,E0,and"LforthenonlinearpartEL.BasedonEqs.( 5.21 )and( 5.22 )thesetwostrainvectorscanbeexpressedasthefollowing, @x@v @y@u @y+@v @x9>>>>>=>>>>>;;"L=8>>>><>>>>:"Lx"LyLxy9>>>>=>>>>;=1 28>>>>><>>>>>:@u @x@u @x+@v @x@v @x+@w @x@w @x@u @y@u @y+@v @y@v @y+@w @y@w @y2(@u @x@u @y+@v @x@v @y+@w @x@w @y)9>>>>>=>>>>>;:(5.27) Oncethevirtualdisplacementsaregivenatanypointwithintheelement,thestrainvariationatanypointcanbedetermined.Therelationshipsbetweenthevirtualdisplacementandstrainvariationare whereB0andBLarethestrain-displacementmatrices,theirexpressionscanbefoundinAppendixA.BasedonEq.( 5.27 )thevirtualworkineachtriangularelementcan


becomputedas whereB=B0+BL,"T="T0+"TL,andFinteistheinternalforceinglobalcoordinates inwhichTisthesametransformationmatrixasshowninEq.( 5.13 ). 5.2.6 ExternalForces Asthemembranechangesitsshapeandpositionunderthepressure,notonlythedirectionofthepressurechangesbutalsotheareaonwhichitactschanges[ 70 105 ].Toaccountforthissituation,weassumethatthesizeofeachniteelementissucientlysmallthatthepressureisuniformovereachelement.Tofurthersimplifythecomputation,weagainuselocalcoordinates.AsshowninFig. 5{2 ,thelocalcoordinatesliesinthedeformedplaneinsteadofintheundeformedplane.InthelocalcoordinatesthepressurePisgivenby, withn=[001]Tistheexternalnormaldirectionoftheelement. ThevirtualworkbytheexternalforceateachelementisZ@V0eDTPdS;


Schematicoflocalcoordinatesforexternalforces. thenetexternalforceateachnodeisobtainedbysimplyrepresentingitbythreeequalforcesateachnode 3001 3001 3T=DTeFexte:(5.33) Hence,inglobalcoordinatesthenodeforceonelementeduetoauniformpressureoverthatelementis 3001 3001 3T:(5.34) 5.2.7 DynamicEquationsfortheMembrane Thersttermoftheleftsideoftheprincipleofvirtualworkequation( 5.11 )canbewrittenas


whereMeisthemassmatrixofanelement, 120s0h02666642IIII2IIII2I377775;(5.36) inwhichIisa3-by-3identitymatrix. Uptonowwehaveobtainedtheexpressionsfortheinternalforce,externalforce,andmassmatrixonanelemente.Itisnecessarytoassembleallthediscreteelementsintoasingleunit.Foranyelementetheprincipleofvirtualworkis Assemblingalltheforcesinoneelementwehavethefollowingequation, Thisequationisforasingleelement,ifwesumallelementsinglobalcoordinates,thenweobtainthesystemofgoverningequationsformembraneunderexternalloads, whereMisapositivedenitemassmatrix,whichremainsconstant,Fintistheinternalforce,andFextistheexternalload. 5.2.8 SolutionoftheDynamicEquations Eitheranimplicitoranexplicitschemecanbeadoptedtofollowthetimehistoryofthesystemofequations.Theexplicitschemeiscomputationallyeconomicforeachtimestepandrequireslessstorage.Itseciencycanbefurtherimprovedbyspecialtechniques,suchaslumpingtechniqueofthemassmatrix.However,itstimestepneedstobesmalltosatisfythestabilityrequirement.Asfarastheimplicitschemeisconcerned,thematrixsystemisupdatedmorethanoncepertimesteptoaccount


forthenonlineareect,whichrequiresmorecomputationspertimestepandmorestorage.Generallyspeaking,implicitalgorithmsareeectiveforstructuraldynamicsproblemsinwhichtheresponseiscontrolledbyasmallnumberoflowfrequencymodes,whileexplicitalgorithmsareecientforwavepropagationproblemsinwhichthecontributionofintermediateandhighfrequencystructuralmodestotheresponseisimportant[ 94 ]. Therearenumerouswaystointegratethestructuraldynamicsequationsintime;amongtheseweutilizetheWilson-method[ 113 ].TheWilson-methodhasacontrollablealgorithmicnumericaldissipation,itisaone-stepmethod,andityieldssecondorderaccuracyintime.Thesecharacteristicsmakeiteasiertocodeandbehavewellinouruid/structureinteractionsystem.SincetheWilson-methodisdissipative,weemployatimestepsmallerthanthatrequiredbytheoverallaccuracycriteriontoavoidexcessivedissipation.Thisapproach,with=1.4,workswellforthepresentproblems. 5.3 NumericalTestofMembraneModel Theinationofasphericalmembraneunderanimposedpressurestepwasan-alyzedexperimentallybyAlexander[ 1 ]andtheoreticallybyBeaty[ 4 ].Recently,anumericalstudywasperformedbyVerronetal.[ 105 ].Thecorrespondingsemi-analyticalsolutionwasreportedbyVerronetal.[ 104 ].Thesestudiesfocusonanincompressiblesphericalhyperelasticmembrane.ThemembranehasaninitialradiusofR0,anundeformedthicknessofh0,andamassdensityof0,respectively.ThemembranematerialisassumedtoobeytheMooney-Rivlinmodel.Therearetwocon-stitutiveparametersinthismodel,namelyc1andc2.Tosimplifyourcomputation,wedeneanon-dimensionalMooney-Rivlinconstantbasedonthesetwoconstitutiveparameters


TheanalyticalgoverningequationforthesphericalmembraneisdetailedbyVerronetal.[ 104 ].Hereweonlygivethedimensionlessform =p2+(1 whereisthecircumferentialprincipleextension,whichisdenedastheratioofthedeformedtotheundeformedradius.Thecorrespondingnon-dimensionaltime,andthenon-dimensionalimposedpressurestepparegivenasfollows: whereP0isthedimensionalpressurestep.Theinitialconditionsofthemembranemayvary;here,wefocusonthecasethatthemembraneisinitiallyinitsequilibriumstate,i.e.,withthefollowinginitialconditions: Onecanobtainthesemi-analyticalsolutiontoEqs.( 5.40 )and( 5.43 )usinganystan-dardsolver.Thedirectnumericalsolutionsofthemembraneequationareobtainedwiththeniteelementmethoddiscussedbefore.Thecomputationalmeshhas816nodesand1568triangularelementsonasemi-sphericalsurface.Inourworkwechooset=5104whichis500timeslargerthantheexplicitschemeusedbyVerronetal.[ 105 ]wheretheyuseanexplicitschemewithatimestepof106forthesameproblem. Allcomputationsareperformedusingthedouble-precision.InFig. 5{3 wecom-paretheniteelementresultswiththesemi-analyticalsolutions.ExceptforonecaseinFig. 5{3 b,allnumericalsolutionsmatchwellwiththesemi-analyticalsolutions.Thelossoflinearityofoscillationsiswellreproduced.Themismatchedcaseshownin


Fig. 5{3 bcorrespondsto=0.1.Forthisvalueofthematerialparameter,aspointedoutbyVerronetal[ 105 ],thereisanunstableequilibriumpointatp0.687.Neartheunstableequilibriumpoint,theniteelementresultscanonlyapproachqualitativelytherealbehavior.Thisisthereasonwhythenumericalresultsarehighlysensitivetotheassignedpressurestep.Asexpected,thelargertheMooney-Rivlinconstantis,thehigherpressurethemembranecansustain. TwoprofoundissuesshouldbeconsideredintheuseoftheWilson-method.First,theWilson-methodisunconditionalstableforlinearproblems.Thetimestepinnonlinearproblemsisoftendictatedbytheerrorsintroducedintheapproximationtreatmentofthenonlinearities;therefore,duecareshouldbetakentochooseanapproximatetimestep.Second,asmalltimestepisdesirablefortheWilson-methodduetoitsexcessivedissipation.AsshowninFig. 5{4 ,alargertimesteptendstocausemorenumericaldamping.InFig. 5{5 weshowtheerrornormoftheimplicitscheme,thegureexhibitsaslopearound1.Therefore,theoverallaccuracyfortimeintegrationisapproximatelyrstorder. Thepresentmembranemodelgivesagoodagreementwiththeanalyticalsolu-tionsformembranedynamicsunderlargedeformation.Itleavesoneissueoffunda-mentalimportancesinceitcannothandlewrinklephenomenonwhichhappenswhenthemembraneiscompressed.Themembranewingissometimescompressed,espe-ciallyattheleadingedge.Therefore,furtherinvestigationisnecessaryinthisarea.


Figure5{3: DynamicinationofaMooneysphericalmembrane:(a)=0;(b)=0:1;(c)=0:25.Solidline:Semi-analyticalsolution;Circle:Numericalsolution.


Theeectoftimesteptothenumericaldampingforp=0.5and=0. Errornormversustimestepfortheimplicitschemeforp=0.5and=0.


6.1 Introduction Membranestructuresarefoundinmanyengineeringpractices.Traditionally,theyareutilizedinparachutesandsailboatstoimproveperformance.Arecentinterestingapplicationisamembranewing-basedMicroAirVehicle(MAV).MAVswithamaximumdimensionof15cmorlessandightspeedaround10m/sprovideexpendableplatformsforsurveillanceanddatacollectioninsituationswherelargervehiclesarenotsuitableorelsedangeroustoreach.However,MAVsaresensitivetotheenvironmentaldisturbances,anditishardtoavoidconsideringtheirsmalldimensions,lightweights,andlowightspeeds.Ourexperienceindicatesthatanaeroelasticmembranewingcanbetteradapttotheatmosphericdisturbances,providesmootherviewingplatform,andmakethevehicleeasiertoy[ 83 ]. Duringight,thesurroundingowexertsaerodynamicforcesonthemembranewingandcausesstructuraldeformationsofthemembranesurface;meanwhile,thewingdeformationsaecttheowpatternandchangetheaerodynamicforcedistri-butions.Themutualinteractionsbetweentheuidowandthestructurecauseauidandstructureinteractionproblem.Fluidandstructureinteractionproblemshaveattractedinterestofmanyresearches.Numerousstudieshavebeenconductedonthesesubjects.Numericalschemesforuidandstructureinteractionproblemsmayvary.However,ingeneral,therearetwokeyelementsinacoupleduidandstructureinteractionsystem.Oneistheuidsolverthatevaluatestheaerodynamicforces,andtheotheristhestructuralsolverthatcomputestheshapedeformation.Besidethese,amovinggridtechniqueisalsoneededtoregeneratethecomputational 53


grid;aninterfacealgorithmisutilizedtoexchangeinformationbetweentheuidandstructuralsolvers. Theinteractionbetweenauidandamembraneusuallycausessubstantialstruc-turaldeformationthatisinherentlynonlinear.Themembranebehaviorisingeneralcomplexandearlierstudiesofinteractionsbetweenuidowsandmembranestruc-turesarelargelybasedonempiricalobservations[ 34 ],two-dimensionalowbasedcomputations[ 36 86 92 93 ],orsimpliedthree-dimensionalanalyses[ 37 38 ].Inourstudywepresentacoupleduidandstructureinteractionsystemtoinvestigatetheinterplaybetweenamembraneandincompressibleviscousow.Inthecoupledsystem,theuiddynamicsisgovernedbythethree-dimensionalNavier-Stokesequa-tionsforincompressibleows,andthestructuraldynamicsisgovernedbyanonlineardynamicmembranemodel.TheNavier-Stokesequationsaresolvedusingapressure-basedmethodonastructuredandmulti-blockgrid.ThemembranematerialfollowsthehyperelasticMooney-Rivlinmodel[ 64 ].Themembranemodelisderivedfromtheprincipleofvirtualwork.ItisspatiallydiscretizedwithatriangularniteelementmethodandtemporallyintegratedwiththeWilson-method[ 113 ].Toaccommodatethemembraneshapechange,theowsolverisfurtheraugmentedwithamovinggridtechniquediscussedinChapter 3 Theinteractionsystemiscoupledthroughtheexchangeoftheaerodynamicforcesandshapedeformationsbetweentheuidandstructuresolvers.Sincetheuidsolverandstructuralsolversdonotsharetheidenticalgridattheinterface,thethinplatespline(TPS)interpolation[ 16 ]isusedtomaptheexternalloadandtheshapedeformationtothestructuralanduidsolvers,respectively.Toavoidphaselagbetweentheuidandstructuralsolvers,subiterationsareperformedtosynchronizetheuidandstructuralsolvers. Inthefollowingwerstintroducetheinterfacealgorithmbetweentheuidandstructuralsolvers,thisisfollowedbyastep-by-stepdescriptionoftheuidand


structureinteractionalgorithm.Thiscouplingisachievedthroughexchanginginfor-mationandsynchronizationbetweentheuidandstructuralsolvers.Intheaerody-namicanalysispart,westudyboththerigidandmembranewings.Intherigidwingstudy,wefocusonthecharacteristicsofthelowReynoldsnumberandlowaspectratiorigidwing,includingboundarylayerseparationandreattachment,tipvortices,andunsteadyphenomenonathighangleofattack.Inthemembranewingstudy,westudytheself-initiatedmembranevibration,itsimpactsonthemembranewingperformance,andwealsocomparethemembraneandrigidwingperformances. 6.2 CouplingbetweentheFluidandStructuralSolvers Whenauidowpassesaexiblestructure,thestructurechangesitsshapeduetotheexternalforcessituatedinow;meanwhile,theshapechangeaectstheuidowaroundthestructureandtheforcedistributionsonthestructure.Thisleadstoauidandstructureinteractionproblem.Usually,theuidproblemandthestructuralproblemhavedierentmathematicalandnumericalproperties.De-pendingontheproblemcomplexity(linearproblemornonlinearproblem,simpleorcomplexgeometry,small-scaleorlargescaleproblem),availablesoftware,andcom-putingplatform,numeroustechniqueshavebeendevelopedtoprovidecoupledCFDandCSDmethods.Oneapproachistotreattheowandstructuralequationsasonemonolithicsystemofequations[ 5 ]andsolvethesystemofequationssimultaneouslyona\monolithic"grid.However,sincemultipletimescalespresentinauidandstructuralsystem,auniedtreatmentmaynotbeecienttohandlethecomputa-tionalstiness.Liuetal.[ 58 ]solvedtheowandstructuraleldssimultaneouslybyafullyimplicitmethod.However,unlikethemonolithicmethod,theirCFDsolverandCSDsolveraresolvedondierentgridsystemsinwhichTheuidandstructuralgoverningequationsareregardedasonesinglesystemoftime-dependentequationsinapseudo-time.


Tobetterhandledierentcharacteristicsintheowandstructuraldomains,mostoftheresearcherssolvedtheuidandstructuralequationsseparatelyandthencou-pledthemthroughasynchronizationprocess.SmithandShyy[ 93 ],Kamakotietal.[ 49 ],GordnierandFithen[ 22 ],KalroandTezduyar[ 47 ],ZwaanandPrananta[ 115 ],andLianetal.[ 53 55 ]coupledtheuidandstructuralsolverswithasubiterationapproach.Theuidsolverandthestructuralsolverfunctionindependentlyontheirowncomputationalgridswiththeirowntimesteps,subiterationsbetweenthesetwosolversareemployedateachtimesteptoavoidthephaselag.Alternatively,Farhatetal.[ 26 28 ]formulatedtheuidandstructureinteractionproblemasathree-eldproblem:theuid,thestructure,andthedynamicmeshthatisrepresentedbyapseudo-structuralsystem.Intheirwork,thearbitraryLagrangianEulerian(ALE)nitevolumeschemeisemployedfortheuiddynamicsequations,aniteelementmethodisappliedforthestructuraldynamicsequations.Partitionedproceduresandstaggeredalgorithmsarethenadoptedforthecoupleduidandstructureinteractionproblemsinthetimedomain. Inuid-structureinteractioncomputations,amovinggridtechniquecanbefruit-fullyemployedtoautomaticallyregeneratethenewCFDgridaccordingforthede-formedconguration.Inaddition,sincetheCFDsolverandCSDsolverdonotnecessarilyshareidenticalgridattheinterface,interpolationsareneededtoallowtheinformationexchange.Smithetal.[ 91 ]evaluatedsuitablemethodstotransferinformationbetweenCFDandCSDgrids.Inourstudyathinplateinterpolationmethod[ 16 ]isadoptedforthispurpose.Thethinplatesplineinterpolationisaglobalinterpolationmethod,itisinvariantwithrespecttorotationsandtranslations,andhenceissuitablefortheinterpolationonmovingsurfaces.Aone-dimensionalversionoftheinterpolationis


whereH(x)isthedisplacementdistributionfunctionoverthedomain,iistheunde-terminedcoecient,Nisthenumberofmonitoredstructuralnodesontheinterface,andxi'saretheirlocations.ItcanbeseenfromEq.( 6.1 )thattheinterpolatedfunctionH(x)hasacontinuousrstorderderivative.Oncethesecoecientsiaredeterminedbasedonthestructuralgridlocations,thedisplacementvectordenedontheCSDgridDismappedtotheCFDgridviaamatrixG whereXisthecorrespondingdisplacementvectorontheCFDgrid.SimilarlythesurfaceforcescanalsobemappedfromtheCFDgridtotheCSDgrid. Uponupdatingtheaerodynamicforcesinthestructuralequationsandprovidingthenewboundaryconditionstotheuidsolver,thetemporaldevelopmentbetweentheaerodynamicandstructuralequationscanbecoordinated.Toconducttime-accuratecomputationsfortheuid-membranedynamics,iterationsbetweentheuidandstructuralsolversshouldbedoneateachtimeinstant.Astep-by-stepcoupledalgorithmisdescribedbelow: 1. GeneratetheCFDandCSDgrids. 2. Setn=0,evaluatepressureP0ontheundeformedwingbasedontheinitialvalues0,u0,and0. 3. n=n+1,andtn=tn1+t. 4. SolvetheNavier-StokesequationsandtransferthepressurePn+1tothestruc-turebasedontheTPSinterpolation. 5. UsethestructuralsolvertoevaluatethedisplacementvectorDn+1. 6. Performsubiterationstosynchronizetheuidandstructuralsolvers,andseti=0. (a) i=i+1;


(b) MapthedisplacementsfromtheCSDgridtotheCSDgridwiththeTPSinterpolationXin+1=GDin+1. (c) UpdatetheCFDgridwiththemovinggridtechnique. (d) UpdatetheJacobiantosatisfythegeometricconservationlaw. (e) Updatetheboundaryconditionsfortheuidsolver. (f) ComputethepressurePin+1eldbasedonthenewCFDgridandboundaryconditions. (g) MaptheaerodynamicloadfromtheCFDgridtotheCSDgridwiththeTPSinterpolation. (h) CalculatethestructuraldisplacementDi+1n+1. (i) Checkforconvergenceforthesubiterationprocess. If`No',returntostepa),otherwise,continue. 7. Returntostep3toprocessnexttimestep. Itshouldbenotedthatafullyconvergedsolutionintheprocessishardtoachieve.Thenumberofsubiterationisbasedonthecomputationalresource,problemnonlinearity,andthetimestep.Inourstudywendtwoorthreesubiterationsareenoughtoalleviatethephaselagerrors. 6.3 MAVWingAerodynamicsandStructuralResponse AsshowninFig. 1{1 themembranewingconsistsofaleadingedgesparandthreechordwisebattensoneachsideofthehalfwing.Themembranematerialisbondedtothesparandbattens.Whilethemembraneisexibleanddoesnotsustainbendingmoment,thecarbonberisrigidwhichprovidessupportforthemembraneandcansustainbendingmoment.Tofullymodelthemembranewing,ideally,oneneedstomodelboththerigidbattensandtheexiblemembrane.Twoprincipaldicultiesarisebydoingthat.First,themembranedoesnotbearbendingmoment,ithasthreedegreeoffreedomsateachnode;thebatten,ifmodeledasabeam,hasvedegreesoffreedomateachnode.Aspecialtreatmentisrequiredtotreatthe


interfacepointwhichbelongstodierentregimesandhasdierentdegreesoffreedom.Second,themembraneismuchthinnerthanthebatten,themassmatrixassociatedwiththesetwodierentmaterialsmakestheassembledmassmatrixill-conditioned.Giventhesefactors,wetreatthebattenasaspecialmembranematerialwithalargerdensitywhiletheotherpropertiesarethesame;thedensityratiobetweenthebattenandthemembraneissettobethree. 6.3.1 ComputationalBackground Inourstudyonlythemembranewingisconsideredwhileignoringthefuselageandpropellers.AschematicgeometryofthewingisshowninFig. 6{1 .Theshadedareashowncorrespondstothefuselagewhichisxed.Thephysicalproblemunderconsiderationisviscousowoverawinginanunboundeddomain.Basedonthefreestreamvelocityof10m/s,therootchordReynoldsnumberis9104.Computa-tionsarerstperformedtoassessthesensitivityofthecomputedresultstothewinglocationandthetypesofboundaryconditionsimposedtotheouterowboundaries.Basically,theboundariesshouldbesetsucientlyfarawayfromthemembranesothattheydonothavemuchimpactonthecomputedaerodynamiccoecients.BasedonourtestsaproperwingpositionischosenandshowninFig. 6{2 .Theupstreaminletboundaryis7cfromtheleadingedge,herecisthechordlengthattheroot.AtsurfaceABCDazerogradientboundaryconditionisimposedforthevelocitycompo-nents;alltheotherboundariesareassignedasaDirichlettypeboundarycondition.Sincethefreestreamvelocityisparalleltothechordwisedirectionandnopropellerismodelled,onlyhalfofthedomainiscomputed,andasymmetricmappingcanbeappliedtotheotherhalfdomain. Gridsensitivitytestsarealsoperformed.Threegridsystemshavebeensystem-aticallychosen,rangingfrom1:8105nodesatthecoarseleveland2:3106nodesatthenelevel.Forthecoarsegridthereare41pointsinthechordwisedirectionand31pointsinthespanwisedirectiononthewingsurface.Thegriddistribution


Figure6{1: Schematicwingwiththreebattensoneachhalfwing. Figure6{2: Computationalset-upandboundaryconditions.


atthecoarselevelaroundthewingisshowninFig. 6{3 .Theresultsfromthenegridareusedasthereference.ThecomparisonamongthesethreegridsystemsissummarizedinTable 6{1 .Thecoarsegridsolutionunder-predictsthelift-to-dragratio(CL=CD)by0.5%,andthedierenceinformdragislessthan0.7%.Asimilarconclusioncanbedrawnforthelift,L.Asexpected,theskinfrictiondrag,DFissensitivetothegriddensity;thecoarsemeshoverestimatesDFbyalmost10%.Anintermediategrid,whichhas6:7105nodesandmorethandoublesthesurfacepointsonthewingsurfacecomparingtothecoarsegrid,givesanimprovedpredictiononskinfrictiondrag.FromTable 6{1 itcanbeseenthatskinfrictiondraghasmuchlesscontributiontothetotaldragthanformdrag.Thisisalsovalidathighanglesofattackwheretheformdragisthedominantfactorinthetotaldrag. Figure6{3: Schematicwinggeometryandstructuredgridwith1.8105pointsand4131onhalfwingsurface. Wefurthercomparethemainowfeatures.Flowseparationsareobservedwithallthesecomputationalgridsatthelowerwingsurface.Weextractthechordwisevelocityproleatthequarterchordoftherootregion.Theboundarylayervelocity


Table6{1: Gridrenementeectsonthecomputedaerodynamiccoecients. GRID 5.10E-1 9.70E-3 7.06 Intermediate6.7105 5.08E-1 9.22E-3 7.16 Fine2.3106 5.04E-1 8.80E-3 7.10 proleisshowninFig. 6{4 .Thenegridhasthelargestreversevelocity,thecoarsegridhasthesmallest,andtheintermediategridhasthethickestseparationzone.Theseresultsindicatethatthegridresolutionnoticeablyaectsthesharpnessofthevelocityprole,aswellastheextentoftheowseparation.Nevertheless,thecoarsegridpredictsthemainfeaturesoftheowwhichisqualitativelyconsistentwiththoseonthenegrid. Figure6{4: Spanwisevelocityprolesatthebottomsurfaceatthequarterchordoftherootregionforrigidwingat=6o. Comparisonsbetweentheintermediateandcoarsemeshesatotheranglesofattackarealsoconducted.ThecomputedaerodynamiccoecientsareshowninFig. 6{5 .Boththeliftanddragcoecientsshowgoodagreements,withalittlebitlargerdierenceathigheranglesofattack.Afewmorewordsshouldbeputforthecomputationalresultsathighanglesofattack.Generallyspeaking,vortexshedding


andboundarylayerseparationoccuratlargeanglesofattack,andtheyintroduceunsteadinesstotheaerodynamicperformanceandthereforeaecttheaerodynamiccoecients.Cummingsetal.[ 10 ]reportedthatatlargeanglesofattacktheunsteadycomputationspredictedmuchlowerliftcoecientsthanthesteadycomputations.Weperformsteadycomputationswhentheincidenceislessthan15owhileconductingbothsteadyandunsteadycomputationswhentheangleofattackislargerthanthat.Athighanglesofattack,trulysteadysolutionscannotbeobtained.The\steady"solutionsareobtainedwithsteadycomputationsoncecertainorderisreducedintheresidual.ThePISOmethoddiscussedinChapter 2 isusedfortheunsteadycom-putation.Wesetthetimesteptobe5105s.However,forourcomputations,thedierencebetweenthesteadyandthetime-averagedliftanddragcoecients,arefoundtobelessthan1%overallthestudiedcases.Fig. 6{5 showsthecomputedaerodynamiccoecientsfortherigidwingwithbothcoarseandintermediatemeshes.Atlowanglesofattack,thecoarsegridpredictsslightlylargerliftthantheinterme-diategrid;thetrendisreversedathighanglesofattack.However,bothpredictthatthestalloccursapproximatelyat51o. Aerodynamiccoecientsversustheangleofattack:(left)liftcoecient;(right)dragcoecient.


Foraproblemascomplicatedasthepresentone,trulygridindependentresultsarediculttoattain.Consideringtheavailablecomputingresourceandtheout-comeoftheaforementionedinvestigation,weusetheintermediategridfortherigidwingsimulationwhileusingthecoarsegridtoillustratethekeyissuesintheuidandexiblestructureinteractionsystem,todemonstratethesalientfeaturesofmem-branewingcomputations,andtoillustratetherelatedaerodynamicandstructuralcharacteristics. 6.3.2 RigidWingAerodynamics VorticalFlowStructuresoftheRigidWing Tipvorticesexistonanitewingduetothepressuredierencebetweentheupperandlowerwingsurfaces.Thetipvortexestablishesacirculatorymotionthatpresentsoverthewingsurfacesandexertsgreatinuenceonwingaerodynamics.Forexample,tipvorticesincreasesthedrag.Thetotaldragcoecientonanitewingatsubsonicspeedscanbewrittenas[ 2 ] whereCD;Pistheformdragcoecientduetothepressure,CD;Fisthefrictiondragcoecientduetoviscosity,eisthespaneciencyfactorwhichislessthanone,ARistheaspectratio,andCD;i=C2L 6.3 )demon-stratesthatinduceddragvariesasthesquareoftheliftcoecient.Furthermore,Eq.( 6.3 )illustratesthatasARisdecreased,induceddragincreases.TheMAVwinginourstudyhasalowaspectratioof1.4;therefore,itisimportanttoinvestigatethetipvortexeectsonthewingaerodynamics.Besidestheireectsonthedrag,tipvorticeshavetwo-foldedimpactsonthelift:


2 ]. 67 ]. Fig. 6{6 showstipvorticesaroundthewingsurfaceatdierentchordwisepo-sitionstogetherwiththestreamlinesattheangleofattackof39o[ 53 ].Thehigherpressurefromthelowersurfacedrivestheowtowardtheupperwingsurfacewherethethepressureislower.Thetwovorticesonbothsidesofthewingdemonstrateclockwiseandanti-clockwiserotations;itseemsthattheowisleadingfromthelowersurfacetotheuppersurface.Fig. 6{7 showsthelowpressurezoneduetothevorticalow.Thepressuredropfurtherstrengthenstheswirlbyattractingmoreu-idstowardthevortexcore;meanwhile,thepressuredecreasescorrespondinglyinthevortexcore.Thelowpressureregioncreatedbythevortexgeneratesadditionallift.Towarddownstream,thepressurerecoverstoitsambientvalue.Thentheswirlingmotionweakensandthediameterofthevortexcoreincreases.Thevortexcorelossesitscoherentstructureatdownstream. Fig. 6{8 visualizestheevolutionofthevorticalowstructurewiththeangleofattack.Thepressuredistributionontheuppersurfaceisalsopresentedinthesamegure.At=6otipvorticesareclearlyvisibleeventhoughtheyonlycoverasmallareaandareofmodeststrength.Theowisattachedtotheuppersurfaceandfollowsthechordwisedirection.Alowpressureregionnearthetip,whichisindicatedwithgreencolor,isobserved.Vorticesstrengthenwiththeincreaseofangleofattack.At=27o,asshowninFig. 6{8 ,tipvorticesdevelopastrongswirlingmotionandentrainthesurroundingowtothevortexcore.Thelowpressureareaincreasesastheangleofattackbecomeshigher. ThepressuredropinthevortexcoreisdemonstratedinFig. 6{9 awherethespanwisepressurecoecientontheupperwingsurfaceatx=c=0.4isplotted.At


Figure6{6: Streamlineandvorticesforrigidwingat=39o.Thevorticalstructuresareshownonselectedplanes. Figure6{7: Pressuredistributionaroundtherigidwinginthecrosssectionswithstreamlinesatangleofattackof39o.




6{8 ).Theliftbythevorticesisnotenoughtomaintaintheincreaseoftheliftcoecientwiththeangleofattack,andstalloccurs.Eventhoughttipvorticesaectthepressuredistributionatthelowersurface,theeectedregionismainlylocatednearthetip.Mostregionsinthelowersurfacearenotaectedbythetipvortices.FromFig. 6{9 bitcanbeseenthatthepressureisalmostuniforminthespanwisedirectionevenathighanglesofattack.Fig. 6{9 isillustrativeinregardtothepressuredistributionsversusthevorticalstructures,however,theyarenotindicativeofthetotallevelofthepressureforce.OneshouldalsonotethatFig. 6{9 indicatesthatthemomentexperiencessubstantialvariationsastheangleofattackchanges. TheLaminarBoundaryLayerSeparationonRigidWing Thelaminarboundarylayer,especiallyunderlowReynoldsnumberconditions,ispronetoseparateunderanadversepressuregradient.Thisseparationisusuallyfollowedbyaquicktransitionfromthelaminartoturbulentow.Iftheadversepressuregradientismodest,theresultingturbulentow,byobtainingenergythroughentrainment,mayreattachtothesurfacetoformalaminarseparationbubbleandthenformsanattachedturbulentboundarylayer.


(a)cpatuppersurface (b)cpatlowersurface Figure6{9: Spanwisepressurecoecientdistributionsatx/c=0.4forrigidwingatdierentanglesofattack. EversinceitsrstobservationbyJones[ 42 ],manyexperimentalinvestigationsonthelaminarseparationbubblehavebeenconducted,asreviewedbyYoungandHorton[ 114 ].Typically,alaminarseparationbubblehasthefollowingcharacteris-tics.Justdownstreamoftheseparationpoint,thereisa\dead-air"region,wheretherecalculatingvelocityissignicantlylessthanthefreestreamvelocityandtheowisvirtuallystationary.Theseparatedowformsafree-shearlayer,whichishighlyunstable.Theseparatedshearlayerquicklyundergoestransitionfromlaminartoturbulentow.Theturbulentowmayreattachtothesurfacebehindavorticalstructureandformsaturbulentboundarylayer.Hence,aseparationbubbleforms.ThedynamicsofalaminarseparationbubbledependsontheReynoldsnumber,thepressuredistribution,thegeometry,thesurfaceroughness,andthefreestreamtur-bulenceintensity.AroughruletodeterminethereattachmentpointisgivenbyCarmichael[ 7 ]thatsaystheReynoldsnumberbasedonthefreestreamvelocityandthedistancefromtheseparationpointtothereattachmentpointisapproximately5104.Ontheonehand,iftheReynoldsnumberislessthan5104,anairfoilwillexperienceseparationwithoutreattachment;ontheotherhand,alongseparation


bubblewilloccuriftheReynoldsnumberisslightlyhigherthan5104.Inthefol-lowing,detailedowstructuresaroundarigidwingarepresented.Bycontrastingthedetaileduidowaroundtherigidandmembranewings,onecanbetterevaluatetheMAVtechnologies.Thebubblestructuresatdierentanglesofattackaredemon-stratedinFig. 6{10 .Representativeinstantaneousstreamlinesattherootsectionontherigidwingareplotted.Atthismoment,thespanwisevelocitycomponentisignoredtomakeaclearrepresentationoftheseseparatedstructures.Atalowangleofattackof6o,theowprimarilyattachestothesurfaceandatinyseparationbubbleisobservedonthelowerwingsurfaceneartheleadingedge.Beginningat45%chordfromtheleadingedge,aweakrecirculationzoneisseenontheupperwingsurface,whichproducesareattachmentlengthofx/c=0.5. Withtheincreaseofincidence,theseparationpointmovesforwardtowardtheleadingedge.Attheangleofattackof15o,asshowninFig. 6{10 ,theseparationoc-cursat39%ofthechordfromtheleadingedge,whichisfollowedbyalong\dead-air"zonebeforeitreachesamaximalreversevelocityof0.26U.Inthe\dead-air"zone,thestationaryreverseowdoesnothavemucheectonthepressuredistribution,whichisprimarydeterminedbythewingcurvature.Theshearlayerreattachmenthappensat89%ofthechord.Basedonthefreestreamvelocityandthedistancebetweentheseparationandreattachmentpoints,theReynoldsnumberisapproximately4.5104whichisslightlysmallerthantheReynoldsnumbersuggestedbyCarmichael.Thereattachmentpointcorrespondstoarapidincreaseinthesurfacepressure,itishighlyunstable,andmovesforwardandbackward. ThestreamwisevelocitycontoursatdierentanglesofattackaredemonstratedinFig. 6{11 .Thevelocityisnormalizedbythefreestreamvelocity.At=6o,themaximumreversevelocityislessthan0.5%ofthefreestreamvelocity.Undersuchacondition,experimentshowsthataconsiderableportionoftheshearlayerremainslaminar[ 21 ].




However,withtheincreaseofangleofattack,themagnitudeofthereversevelocityincreasesaswell.At=27o,thereverseowcomponentoftheseparationbubblereachesamaximalvelocityof0.49U.Undersuchasituation,CromptonandBarrett[ 9 ]showedthattheshearlayerisparticularlyenergeticsincemostoftheshearlayeristurbulent.Eventhoughtheowontheuppersurfaceneartheroothasseparatedfromalargeportionofthesurface;however,theliftstillincreaseswiththeangleofattack.Tworeasonsmightliebehindthis.First,tipvorticesgeneratesuctionareasontheupperwingsurfacewhichprovideadditionallift.Second,eventhoughowseparatesneartheroot,itstillattachestotheotherregionofthewingsurface.Whentheangleofattackisincreasedfurther,at=51o,massiveseparationappears,andstalloccurs.Meanwhile,asseenfromFig. 6{12 ,vortexsheddingisobservedneartheroot.Atinyseparationbubbleemergesattheleadingedgerst,itthenobtainsenergyfromtheshearlayerandincreasesitssize.Thisbubblethenmergeswithanotherbubbledownstreamtoformalargerbubble.Thelargerbubbleisnotstable,itbreaksintotwosmallbubblesnearthemaximalcamberposition,andvortexsheddingbegins. Forlowaspectratiowings,tipvorticeshaveconsiderablecontributionstothelift.Thisscenarioisquitesimilartothatofdeltawings[ 23 ].Wehaveobservedthatthelowaspectratiowingsuerslessfromtheseparation.AsseenfromFig. 6{5 a,theliftkeepsalinearrelationshipwiththeangleofattackevenatveryhighanglesofattack.ExperimentsbyTorresandMueller[ 97 ]onlowaspectratiowingshavesimilarndings.






UnsteadyPhenomenonatHighAnglesofAttack Bothvortexsheddingandboundarylayerseparationoccuratlargeanglesofattack,andtheyintroduceunsteadinesstotheaerodynamicperformance.Asmen-tionedbefore,weperformsteadycomputationswhentheincidenceislessthan15owhileconductingbothsteadyandunsteadycomputationswhentheangleofattackislargerthanthat.Interestingly,thedierencebetweenthesteadystatecomputa-tionsandthetime-averagedunsteadycomputationsaresmallevenatlargeanglesofattackinwhichunsteadyphenomenonsuchasvortexsheddingareprominent.Thedierencesintheliftcoecientandthelift-to-dragratioarefoundtobelessthan1%(Fig. 6{13 ).Noticeabledierenceappearsonlywhentheangleofattackbecomessucientlyhigh. ComparisonoftherigidwingCLbetweenthesteadyandunsteadycomputations. Thepressuredistributionsarealsocomparedattherootsection,whereowseparationusuallyappearsrstduetothelargecamberinthepresentMAVwingdesign[ 34 ].Fig. 6{14 showsthatat=6othetimeaveragedpressurecoecient


Figure6{14: Thecpcomparisonontherigidwingattherootbetweenthesteadyandunsteadycomputations:(left)=6o;(right)=15o. matchescloselywiththesteadystateresult.ThisndingisconsistentwithFig. 6{10 andFig. 6{11 whereaveryweakrecirculationisseenatthisangleofattack.However,at=15o,cleardierenceisshownafterx/c=0.6whichapproximatelycorrespondstothelocationofthemaximalreversevelocityinFig. 6{11 .Thetimeaveragedvalueyieldsasmoothpressuredistribution;thevariationinthesteadyresultisapparentfromtherecirculationzone.Velocitydistributionsshowthesametrendasthepressuredistributions.At=6o,asseenfromFig. 6{15 ,thereisalmostnodierenceinthechordwisevelocitydistribution.However,at=15o,nodierenceneartheleadingedgewherenoseparationoccurs,butcleardierenceisshownaftertheseparationbubble(Fig. 6{16 ). 6.3.3 MembraneWingDynamics TimeSynchronization Beforepresentingthemembranewingresults,somekeyissuesincoupledmembrane-uidinteractionsarerstaddressed.Specically,thetimesynchronizationbetweentheuidandstructuralsolutions,andthegeometricconservationlawinregardtothemovinggridmethodareemphasized.


Figure6{15: Comparisonofchordwisevelocityprolesontherigidwingbetweenthesteadyandtime-averagedunsteadycomputationsattheangleofattackof6o.Velocityprolesareattwochordwiselocations,andbothattherootsection. Figure6{16: Comparisonofthechordwisevelocityprolesontherigidwingbe-tweenthesteadyandtime-averagedunsteadycomputationsatangleofattackof15o.Velocityprolesatattwochordwiselocations,andbothattherootsection.


Thetimesynchronizationinourstudymeansthatiterationsshouldbeenforcedbetweentheuidandstructuralsolversateachtimesteptoavoidphaselagerrors.Weshallsimulatetheuidandstructureinteractionbyintegratingtheuidandthestructuralsolverswithsubiteration.Todemonstratetheimportanceofthetimesynchronizationbetweentheuidandstructuralsolvers,WecomparetheresultswithandwithoutsynchronizationinFig. 6{17 thatdepictsthedisplacementhistoryofatrailingedgenode.Thecomputationwithoutsynchronizationfailstocontinuewhilethatwithsynchronizationkeepsongoing.Thedisplacementhistorywithoutsynchronizationmatcheswellwiththatwithsynchronizationattheveryearlystage.However,thelaggingerrorsaccumulatewithtimeandeventuallyresultinamuchlargerdisplacementthanthatwithsynchronization.GordnierandVisbal[ 23 ]havealsoreportedtheimportanceofthetimesynchronization.ThevelocityandpressurehistoriesofthesamenodewithandwithouttimesynchronizationarecomparedinFig. 6{18 andFig. 6{19 ,respectively.Beforethecomputationdiverges,thevelocitywithoutsynchronizationismorethantwotimeslargerthanthatwithsynchronization.Thelargervelocityofthemembranenodemeansthatthemembranenodeobtainingmoreenergyfromthesurroundingowthanthatwithsynchronization.ThehighvelocityisbelievedtobeassociatedwiththesuddenpressuredropshowninFig. 6{19 ,whichcausesthewrinkleofthemembraneandthefailureofthestructuralsolver. TheGeometricConservationLaw Thegeometricconservationlaw[ 100 ]isanotherimportantfactorwhenthemov-inggridtechniqueisemployed.Fig. 6{20 ashowsthetimehistoryofCL/CDforthemembranewingat=6o.Withoutsatisfyingthegeometricconservationlaw,irreg-ularspikesareobservedinthecourseofcomputations.However,thecomputationsarebetterregulatedoncethegeometricconservationlawisenforcedasinFig. 6{20 b.Farhatetal.[ 27 ]arguedthatthegeometricconservationlawwasanecessarycon-ditiontomaintainthestabilityofaschemeutilizedinmovingboundaryproblems.


Figure6{17: Eectsoftimesynchronizationonthedisplacementofatrailingedgepointformembranewingat=6o. Figure6{18: Eectsoftimesynchronizationonverticalvelocityofatrailingedgepointformembranewingat=6o.


Figure6{19: Eectsoftimesynchronizationonpressuredistributiononatrailingedgepointformembranewingat=6o. Wedonotcomeacrossinstabilityinourstudy.Twopossibilitiesliebehindthat:oneisthesmalltimestepweusefortheuidsolver,theanotheristherelativelysmalldisplacementinthemembranewing.Nevertheless,thecomputationsbehaveerraticallywhenthegeometricconservationlawisnotenforced. Self-initiatedMembraneVibration Earlierworksinmembranewings[ 88 92 ]werebasedonconsiderationsthatthemembranewasmasslessandtherewasnotimedependentmovementinthesteadyfreestream.Experiments[ 111 ]showedthatthemembraneexperiencedhighfrequencyvibration.Byadoptingadynamicmembranemodel,ourcomputationsshowthatthemembranewingalsoexertshighfrequencyvibrationsinthesteadyfreestream.Fig. 6{21 demonstratesthedisplacementhistoryofthetrailingedgewithtime.Themaximaldisplacementduringthatperiodoccursnearthetip.Thedisplacementisnormalizedbythemaximalcamberthatisabout0.9cm.Toinvestigatethevibrationfrequency,wechooseapointbetweenthebatten1andbatten2onthetrailingedge


Figure6{20: EectofthegeometricconservationlawonCL/CDformembranewingcomputationat=6o:(left)notsatisfyingthegeometricconservationlaw;(right)satisfyingthegeometricconservationlaw asanexample.ItsdeectionhistoryisshowninFig. 6{17 ,andtheestimateddom-inatedfrequencyisaround120Hz(Fig. 6{22 ).Analysesatotherpointsalsoshowthefrequenciesaround120Hzthatiscomparabletotheexperimentalfrequencyof140Hz.[ 111 ]ConsideringthattheexperimentalconditionsandthedetailedwingcongurationbetweenourstudyandthosereportedbyWaszaketal.[ 111 ]arenotidentical,thequalitativeagreementbetweencomputationsandexperimentsinthisregardisdeemedsatisfactory.Liu[ 60 ]reportedthatunderatypicalwinggustsitua-tiontheenergywasmainlylocatedinthelowfrequencyrangeofseveralHzorlower.Themembraneuctuationisnotexpectedtocausesensitiveresponsetothevehiclesincethemembraneuctuatesinatimescalemuchfasterthaneitherthevehiclecontrolscaleortheexpectedwinggusttimescale. Whenthemembranesurfacevibrates,thevelocityonthewingsurfaceisnolongerzero.Aconsiderablevelocitycomponentisobservedatthewingsurface.TheverticalvelocitycontoursattwodierenttimeinstantsareplottedinFig. 6{23 .Fig. 6{23 acorrespondstoastagewhenthemembranewingmovesupundertheaerodynamicforces.Atthetrailingedge,theverticalvelocityisabout2%ofthe


Figure6{21: Timehistoryoftrailingedgedisplacementformembranewingat=6o.Thecamberattherootis0.90cm. Figure6{22: Spectrumanalysisofthetrailingedgepointvibrationformembranewingat=6o.


Figure6{23: Normalizedverticalvelocitycontourattwodierenttimeinstantsformembranewingat=6oontheuppersurface. freestreamvelocity.Anegativevelocitycomponentisalsoobservedneartheleadingedge.Fig. 6{23 bapproximatelycorrespondstoastageofmaximaldisplacement.Anegativevelocityneartherootindicatesthatthemembranebeginstomovedown. ComparisonbetweentheMembraneandRigidWings Withthesamefreestreamcondition,asshowninTable 6{2 ,theexiblewingexhibitsslightlylessliftcoecientthantherigidoneat=6o.ThedierenceinCL/CDisalsosmall.Athigherangleofattackof15o,themembranewinggeneratesaliftcoecientabout2%lessthantherigidwingdoes;however,itsCL/CDis1.5%morethantherigidone.Thesedierencesliebehindthehighfrequencyofthemembranevibration.Undertheexternalforce,themembranewingchangesitsshape.Thisshapechangehastwoeects.Ontheonehand,itdecreasestheliftbyreducingtheeectiveangleofattackofthemembranewing;ontheotherhand,itincreasestheliftbyincreasingthecamber.Ournumericalndingsshowthatmembraneandrigidwingsexhibitcomparableaerodynamicperformancebeforestalllimit,whichwasalsoexperimentallyobservedbyWaszaketal.[ 111 ]. Fig. 6{24 showsthetimeaveragedverticaldisplacementsofthetrailingedgeat=6oand=15o.Thedisplacementsarenormalizedbythemaximalcamberofthe


Table6{2: Aerodynamiccomparisonbetweenmembraneandrigidwings. MembraneWing RigidWing 6o 7.06 0.532 15o 3.88 0.954 wing.At=6othedeectionisabout15%ofthemaximalcamberandincreasestomorethan20%at=15o.Duetothetrailingedgedeectiontheeectiveangleofattackofthemembranewingislessthanthatoftherigidwing.ThespanwiseanglesofattackbetweenrigidandmembranewingsunderthesameowconditionandwithidenticalinitialgeometriccongurationsareshowninFig. 6{25 .AswecanseefromFig. 6{25 a,therigidwinghasanincidenceof6oattherootanditmonotonicallyincreasesto9:5oatthetip;themembranewingsharesthesameanglesofattackwiththerigidwinginthe36%oftheinnerwing;however,theeectiveangleofattacktowardthetipislessthanthatoftherigidwing.Atthetip,theangleofattackofthemembranewinglowersbyabout0.8o.Fig. 6{25 bcomparestheangleofattackat=15o,anditshowsthattheeectiveangleofattackofthemembranewingisabout1olessthanthatoftherigidwingatthetip.Thereducedeectiveangleofattackcausesthedecreaseinthelift. Themembranewingshapechangehasanothereect.Thepressuredierencebetweentheupperandlowerwingsurfacesinatesthemembranewing'scamber,whichisshowninFig. 6{26 .Twoairfoilshapesatdierentspanwisepositionsareplottedtogetherwiththecorrespondingrigidwingshapes.Thecamberincreaseismorevisibleintheouterwingthanthatintheinnerwing.ThisisconsistentwithFig. 6{24 wherethemaximaldisplacementoccursnearthetip.Asexpected,theincreaseislargerat=15othanthatat=6o. Thewingsurfacemovementaectstheoverallpressuredistribution.Fig. 6{27 comparesthetimeaveragedpressurecontoursbetweentherigidandmembrane


Figure6{24: Averageddisplacementofthemembranewingtrailingedge:(left)=6o;(right)=15o. Figure6{25: Comparisonofthespanwiseanglesofattackbetweentherigidandmembranewings:(left)=6o;(right)=15o. wings.Asexpected,thedierencesaremostlylocatedatthetrailingedgewherelargemovementisobserved.ThesedierencesarealsoreectedinthechordwisepressuredistributionsshowninFig. 6{28 .Attherootsectionthereisalmostnodierencebetweentherigidandmembranewingssincethefuselageisxed.How-ever,thedierencebecomescleartowardthetip;atz/Z=0.37,thedierenceisvisible.Atz/Z=0.69,thetimeaveragedmembranewingshiftthelowestpressurepointdownstreamcomparingwiththerigidwingresult.Thisshifthelpstodelayowseparation.


Figure6{26: Timeaveragedairfoilshapesatdierentspanwisepositionsforthemembranewingat=6oand15o. Figure6{27: Comparisonofcpcontoursat=6o:(left)timeaveragedmembranewing;(right)steadyrigidwing.


(a) (b) (c) Figure6{28: Chordwisecpcomparisonatdierentspanwiselocationsbetweentimeaveragedmembranewingandsteadyrigidwingat=6o:(a)z/Z=0;(b)z/Z=0.37;(c)z/Z=0.69.


InChapter 6 wehaveperformedtheaerodynamicstudyofthewingaerodynamicsofboththerigidandexiblewings.Inthischapterwewilldevelopasystematicapproachtoenhancethemembranewingperformance. 7.1 Introduction Itisdesirablenotonlytounderstandthepertinentphysicalphenomenaofthemembranewingdynamics,buttoimproveitsperformancebasedontheunderstand-ing.InChapter 6 wehaveinvestigatedbothrigidandmembranewingaerodynamics.Thewingshapeadoptedisfromanexistingdesign.Onequestionthenarisesnatu-rally:isthatthebestwingdesign,or,canweimproveitsperformance.Toanswerthisquestion,itisworthwhiletopauseandtakealookathowtheexistingwingisde-signed.Thedesignofthewingismostlybasedonobservationsandatwo-dimensionalanalysistooltotesttheperformanceofdierentairfoils.Basedonthetests,oneair-foilshapeisthenchosentoformathree-dimensionalwing.Thewingshapeisthentestedwithight.Inthistrialanderrorprocess,thethree-dimensionaleectsandthemembranematerialpropertiesarenotconsidered.Asystematicapproachisnecessarytoimprovethewingperformance. Therearemanyexistingmaturetechniquesforshapeoptimization.Couldwesimplyuseoneofthemtoachieveourgoal?Theanswerisnegative.AsdiscussedinChapter 6 ,theMAVwingperformsunderthelowReynoldsnumbercondition,anditexhibitsahighfrequencyvibrationevenundersteadyfreestreamconditions.Therefore,theoptimizationofamembranewingintroducesmultiplechallenges.First,coupled,timedependentsimulationsoftheinteractionbetweenviscousuidowandexiblestructureareveryexpensive.Second,anecientandautomaticgrid 88


regenerationisessential.Todealwiththesediculties,weproposethefollowingstrategies: 6 32 51 63 ],inourstudythesur-rogatemodelistherigidwing.Weshalloptimizetherigidwingunderagivenowcondition,andevaluatetheperformanceofaexiblewingwithoutputfromthesurrogatemodel. Inthefollowing,werstintroducetheoptimizationproblem.Thisisfollowedbyanoptimizationalgorithminwhichweshowhowtointegratetheanalysistool,themovinggridtechnique,andtheoptimizerDesignOptimizationTools(DOT)[ 12 ]intoasystemtodoshapeoptimization.Theoptimizationprocessbeginswitharigidwingsurrogateoptimization,thentheexiblewingperformanceisassessedwiththeoptimizedshape.Optimizationresultsandtheaerodynamiccharacteristicsofbothoptimizedrigidandmembranewingsarethenhighlighted. 7.2 OptimizationandComputationalFramework Theobjectiveistoimprovethelift-to-dragratioofamembranewing.Thebaselinemembranewingshapeselectedisanexistingdesign[ 33 34 ]devisedfromsolutionsprovidedbyXFOIL[ 15 ]andempiricalmodicationbasedonighttests[ 34 ].XFOILprovidestwo-dimensionaluidowsolutionsbasedoncoupledthinlayerandinviscidowmodels.AschematicbaselinewingshapeisshowninFig. 6{1 inChapter 6 .Thebaselineshapehasthreecarbon-ber-madebattensthatarelabelledalsoinFig. 6{1 .Theshadingareacorrespondstothefuselage.Thefuselageismadeofrigidcarbonberprepregclothwhichdoesnotdeformduringtheight.Aexible


latexmembranecoversthetopsurfacewhichchangesitsshapeunderaerodynamicforces. Thebaselinewinghasavariablespanwisecamber.Thenominalcamberofthebaselinedesignis7.5%attheroot,and2%atthetip.Thewinghasamaximalchordlengthof13.7cm,ameanchordof10.5cm,andaspanof15.2cm.Cisusedtoindicatethechordlengthatdierentspanstations,andZforthehalfspan.Theangleofattackisdenedinreferencetotherootgeometry.Atthedesignpoint,theMAVieswithaspeedof10m/sandanangleofattackof6o.Theobjectiveistomaximizethelift-to-dragratioatsuchaightcondition.AtsuchaconditiontherootchordReynoldsnumberis9104. Sixdesignvariablesarechosen,theyarethey-coordinatesatsixpointsonthesurface.Threearelocatedattherootwhosepositionsareapproximatelyat10%,27%,and77%ofthechord,asseeninFig. 7{1 .TheotherthreearelocatedatBatten2withthesamerelativechordpositions. Thewingsurfacecanbeimplicitlyrepresentedbyafunctiony=f(x;z),wherex,y,andzrepresentthewingsurfacecoordinates.WeindicatethedesignvariablewithacapitalletterYiwhosevalueistheperturbationoftheycoordinateatpointiwithrespecttoitsbaselinevalue.Weinterpolatetheperturbationsfromthesixdesignvariablesoverthewingsurfacewiththethinplatespline(TPS)interpolationmethoddiscussedinChapter 6 .Therefore,theinterpolatedperturbationsoverthewingsurfaceis whereYisavectorindicatingthedesignvariablevalues,Gistheinterpolationoperator,andyistheinterpolatedvalueoverthewingsurface.TheperturbedwingsurfaceisthereforerepresentedbythesummationofthebaselineandtheinterpolatedvalueinEq.( 7.1 ),ynew=ybaseline+y.OneadvantageofthisperturbationapproachisthatwecanrecovertheoriginalshapewhentheperturbationvectorYiszero.In


(a) (b) Figure7{1: Acrosssectionofmembranewingalongabatten:(a)straightlinesindicatinghowconvexityconstraintisapplied;(b)eectofycoordinateperturbationatapoint.Thedesignpointsarelocatedat10%,27%,and77%ofthechord.

PAGE 100

ouroptimizationprocess,boththeleadingandtrailingedgesarexedtomaintaintheidenticalangleofattackbetweenthebaselineshapeandtheoptimizedshape. Theoptimizationprobleminvestigatedcanbeformulatedasfollows. MaximizeCL=CDSubjectto1:CLCLbaseline2:Convexityconstraint:Y1+y1y0 Constraint1requiresthattheliftcoecientbenolessthanthatofthebaselinewing.Constraint2maintainstheconvexityoftheairfoil(seeFig. 7{1 a)sotoeliminateob-viouslyinappropriateshapes.Onlyoneofthesixconstraintsisshownforillustrationpurpose.Constraint3givesthelowerandupperboundsofthedesignvariableswhosevaluesarelistedinTable 7{1 Table7{1: Designvariablesboundsandtheiroptimalvaluesat6o. Parameter Initialdesign(cm) Lowerbound(cm) Upperbound(cm) Case1:opti-malvalueswith600iterationsperanalysis Case2:optimalvalueswith1000iterationsperanalysis -0.40 0.00 -0.258 -0.259 -0.40 0.00 -0.400 -0.400 -0.20 0.20 -0.152 -0.139 -0.20 0.20 -0.100 -0.095 -0.10 0.20 -0.082 -0.068 -0.10 0.20 0.135 0.121 Analysis 97 108 Cycle 6 6 0.529 0.529 0.529 7.55 7.55 Numberofanalyseswithconvexityviolations 13 22

PAGE 101

TheoptimizationprocedurebeginswiththebaselinedesignbycallingtheNavier-Stokesequationssolverrst.Theoptimizerthenperturbsoneormoredesignvari-ablesforapotentialdecliningdirection.Afterthat,theperturbationswillbeinter-polatedwiththeTPSmethodshowninEq.( 7.1 ).AnewcomputationalgridisthengeneratedbasedonthenewcongurationwiththemovinggridtechniquediscussedinChapter 3 .TheNavier-Stokesequationssolveriscalledtoevaluatetheobjec-tivefunctionbasedonthenewcomputationalgridagain.Adesigncycleincludestheevaluationofderivativesoftheobjectivefunctionandconstraints(herebynitedierences)andaonedimensionalsearchinadirectionselectedbasedonthederiva-tives.ThesearchterminatesafterafewcycleswhentheoptimizerDOTconvergencecriteriaaresatised.Foranobjectivefunctionwithndesignvariables,atleastnanalysesareneededtoevaluatethegradientoftheobjectivefunctionwithoutcon-straints.Thenumberofanalysisincreasesforconstraintoptimization.WiththeTPSinterpolation,theinterpolatedsurfacehasasecond-ordercontinuity.And,asshowninFig. 7{1 b,theinterpolatedsurfacepassesthroughallthecontrolpoints. Toevaluatetheaerodynamicsperformanceofthemembranewing,weusethecoupledmodeltoperformauidandstructureinteractioncomputation.Thiscou-pledmodelconsistsofthreemodules:aNavier-Stokesequationssolver,amembranestructuralsolver,andaninterfacingtechnique.Theuidsolverisapressure-basedmethodforincompressibleowthatwementionedinChapter 2 ;thestructuralsolverisaniteelementapproachforthedynamicresponseofamembranematerialwhichobeystheMooney-Rivlinhyperelasticmodel[ 64 ]. 7.3 ResultsandDiscussion 7.3.1 SensitivityAnalysis Whenemployingagradient-basedoptimizationmethod,akeyquestionisthesensitivityoftheaerodynamiccoecientsinresponsetotheminutegeometricchange.ToanswerthisquestionwerstrecallthegridsensitivitystudyinChapter 6 .Inthat

PAGE 102

studyweconcludethatthefrictiondragissensitivetothegriddistribution.However,theliftandtheformdragonlyhavemodestvariationswiththegridrenement.Sincethefrictiondragisonlyasmallportionofthetotaldrag,thecoarsemeshisusedhereagainforshapeoptimization.Inthefollowing,weinvestigatehowtheuidsolver'sconvergencecriterionaectstheoptimizationresults. Dierentconvergencecriteriaareprobedtoassesstheirimpactsontheopti-mizationoutcome.TheimplicitSIMPLECmethodisusedfortheNavier-Stokesequations.Twoiterationnumbersarechosentointheoptimizationprocess:case1uses600iterationsperanalysis,whilecase2uses1000iterationsperanalysis.Inbothcasestheoverallresidual,denedasthesumoftheabsolutevaluesofthemassandmomentumuximbalanceofallcomputationalcell[ 88 ],isconsistentlyreducedbyanorderoffour,withCase2slightlymore.Bothcasesyieldacceptablesteadystatesolutions,withminordierencesintheresults.Intheoptimizationprocess,itisobservedthataftersixdesigncyclesbothcasesresearchthesamelocaloptimalalthoughtheirconvergencepathsaredierent.Eventhoughtheoptimalvaluesarethesame,thereissmalldierenceinthegeometryatthetip.Weplottheconver-gencehistoriesinFig. 7{2 .TheoptimizationresultsaretabulatedinTable 7{1 .InthefollowingweexclusivelyfocustheoutcomefromCase1asoursurrogatetostudythemembranewingperformance. Inourcomputations,weimposetheconvexityconstraintstoeliminateobviouslyunacceptableshapes.Inthesearchprocess,occasionally,violatingtheconstraintscouldacceleratetheconvergencetowardtheoptimum.However,frequentviolationswilldegradetheeciencyoftheoptimizationprocess.AsshowninTable 7{1 ,thereare13and22analyseswhichviolatethegeometryconstraintsforCase1andCase2,respectively.

PAGE 103

Figure7{2: Optimizationdesignhistorywithdierentnumberofiterationsperanal-ysis. 7.3.2 PerformanceoftheOptimizedRigidWing Thebaselineshapehasthelargestcamberattheroot,about7.5%ofthechordlength;itdecreasesto2%atthetip.Theanalysisofthebaselinewingshowsthatthelargecamberneartherootcancauseowrecirculation.Theoptimizationreducesthecamberneartherootwhileincreasingthecambernearthetip,asseeninFig. 7{3 Theoptimizedshapehasamoreuniformcamber,graduallydecreasingfrom4.8%attherootto4%atthetip.Thismodicationresultsintheeliminationoftherecirculationintherootregion.Fig. 7{4 showsthestreamwisevelocityprolesofbothbaselineandoptimizedrigidwingsneartherootonthepressureside.Itdemonstratesthatthereverseowassociatedwiththebaselinewingshapehasnowbeeneliminatedintheoptimizedshape.Fig. 7{5 illustratesthespanwisepressuredistributionsofbothwings.Attheroot,thebaselinewingyieldsacross-overpressuredistributionattheleadingedge,whichdecreasesthetotallift.Theoptimizedshapeeliminatesthecross-over.Inaddition,theadversepressuregradientonthelower

PAGE 104

(a) (b) (c) Figure7{3: Comparisonbetweenoptimizedshapeandbaselineshapeatdierentspanpositions.Thecamberdecreasesfrom4.8%attherootto4%atthetip:(a)z/Z=0;(b)z/Z=0.4;(c)z/Z=0.8.

PAGE 105

Figure7{4: Comparisonofstreamwisevelocityprolesattherootbetweenbaselineandoptimizedshapes. surfaceissmootherfortheoptimizedshapethanforthebaseline;thishelpsexplainwhytheoptimizedwingdoesnotexperienceowseparationinthatregion. ThespanwiseL/DdistributionisdepictedinFig. 7{6 .Theaerodynamicim-provementislargelyrealizedinthe70%oftheinnerwing.Withareducedcamberintherootregion,thedragcoecientdecreasesnoticeablythere,whichleadstoahigheroveralllift-to-dragratio. Theoptimizationissetattheangleofattackof6o.Toprobetheperformanceoftheoptimizedshapeatotheranglesofattack,theaerodynamiccharacteristicsfordierentanglesofattackareshowninFig. 7{7 .Overall,theimprovementismoresubstantialatthemodestangles.ThespanwiseL/DdistributionsatdierentanglesareshowninFig. 7{8 .Consistently,theimprovementislargelyrealizedfromloweringtheformdrag. 7.3.3 FlexibleWingOptimization Withthesurrogateoptimizationoutput,themembranewingperformanceisevaluatedbasedonthecoupleduid-structureinteractionmodelinChapter 6 .To

PAGE 106

(a) (b) (c) Figure7{5: Comparisonofpressurecoecientsonbaselineandoptimizedwingsatdierentspanwiselocationsat=6o:(a)z/Z=0;(b)z/Z=0.4;(c)z/Z=0.8.

PAGE 107

Figure7{6: Comparisonofspanwiseaerodynamiccoecientsdistributionbetweenthebaselineandoptimizedwingsat=6o.

PAGE 108

Figure7{7: Comparisonsofthebaselineandoptimizedrigidwingsatdierentanglesofattack.

PAGE 109

(a) (b) Figure7{8: Comparisonofspanwiselift-to-dragratiobetweenthebaselineandopti-mizedrigidwingsatdierentanglesofattack.(a)=3o;(b)=9o.

PAGE 110

investigatethespanwiseforcedistributionatdierenttimeinstants,theoverallaero-dynamicsofthemembranewings,forbothbaselineandoptimizedshapes,issum-marizedinFig. 7{9 .Consistentwiththerigidwingcase,showninFig. 7{7 ,theoptimizedshapeimprovestheL/Dconsistentlywithmoresubstantialimprovementatmodestanglesofattack.Furthermore,betweentherigidandmembranewings,themembranewingvarieslessintheL/Dversustheangleofattack,indicatingthatitcanoeramoresteadyightbehaviorthantherigidwing.Thisnalndingisconsistentwiththatobservedinighttests[ 34 ]. 7.4 SummaryandConclusions Wehaveperformedanoptimizationofaexiblemembranewingusingarigidwingassurrogate.Ananalysisofthenaldesignhasveriedthattheexiblewingcanbeimprovedbyoptimizationofarigidwingwiththesameun-deformedgeometry. Comparedtothebaseline,theoptimizationleadstolowercamberneartherootandhighercambernearthetip,whilestillleavingthecamberslightlyhigherattherootthanatthetip.Thelift-to-dragratioisimprovedoverarangeofanglesofattack.However,atlargeanglesofattack,theimprovementwiththeoptimizedshapediminished,theL/Dofthemembranewingvarieslesswithchangesintheangleofattack,indicatingthatitcanoeramoresteadyightbehaviorthattherigidwing.Theimprovementislargelylocatedwithin70%oftheinnerwing.Theliftoftheoptimizedwingremainsthesameasthebaselinewing,eventhoughitscamberreducesattheroot.Theimprovementinaerodynamicsoftheoptimizedwingislargelyrealizedvisreducingtheformdrag.Forbothrigidandmembranewings,theoptimizationimprovestheL/Dconsistently.

PAGE 111

Figure7{9: Theoptimizedmembranewingperformanceatdierentanglesofattack.

PAGE 112

WehaveinvestigatedMAVwingaerodynamicsandshapeoptimizationwithapplicationstoMAVs.Eveninthesteadyfreestream,themembranewingexhibitsself-initiatedvibration.Toaccuratelysimulatetheinteractionbetweentheexiblemembranestructureanditssurroundingviscousow,wehaveperformedcoupleduidandstructurecomputations.Intermsofcomputationalrequirementsforthecoupleduidandstructureinteractions,wehaveenforcedthetimesynchronizationbetweentheuidandstructuralsolverstoreducethephaselagerror.Wehavealsoconrmedtheimportanceofthegeometricconservationinthecontextofmovinggridproblems. InregardtoMAVaerodynamics,tipvorticesplayanimportantroleforthelowaspectratiowing.Ontheonehand,tipvorticeslowertheliftbyreducingtheeectiveangle;ontheotherhand,tipvorticesprovideadditionalliftbyforminglowpressureregionsontheupperwingsurface.FlowontheupperwingsurfaceispronetoseparateunderlowReynoldsnumbercondition.Eventhough,MAVsmaintainaveryhighhighstallangleofattack.Tipvorticesandowseparationcauseunsteadybehaviors.Vortexsheddingoccursathighanglesofattack;however,thevortexsheddingdoesnothavemuchimpactontheaerodynamiccoecients. WehaveconductedCFD-basedoptimizationinamethodicalmanner.wehaveoptimizedamembranewingusingtherigidwingasthesurrogatemodel.Comparingwiththebaseline,theoptimizedwinghasalowercamberneartherootandahighercambernearthetip.However,theoptimizedwingstillleavesthecamberslightlyhigherattherootthanatthetip.Eventhoughtheoptimizedwinghaslowercamber 104

PAGE 113

thanthebaselineattheroot,itmaintainsthesameliftasthebaseline.Theim-provementinthelift-to-dragratioislargelyrealizedviareducingtheformdragandabetterpressuredistribution.Ouroptimizationhasorientedatanangleofattackof6o;testsoverarangeofanglesofattackshowdierentlevelsofimprovement.Atlargeanglesofattack,theimprovementwiththeoptimizedshapediminishes.Atlowanglesofattack,theimprovementismostlyachievedwithin70%oftheinnerwing. Muchprogresshasbeenmadeinthisexcitingresearcharea.Withthecontinuousimprovementinsimulation,materials,fabrication,andmeasurementtechniques,aswellasthefastdevelopmentinmicro-systems,signicantpotentialexiststofurtheradvancethemicroairvehicleasavehicleplatformforperformingnumerousmissions,andasatechnologyplatformforevaluatingsystem.ThefutureworkinthisareawillbeacombinednumericalandexperimentalstudyonthelowReynoldsnumberandlowaspectratioMAVaerodynamics.Theprincipalfeaturesoffuturestudyare ThesearethekeyelementsforanintegratedstudyoftheMAVaerodynamics.Thefulllmentofthisworkhasbothpracticalandtheoreticalsignicance.Practically,thisworkwillimproveourunderstandingontheMAVaerodynamics,enhancetheMAVperformance;theoretically,itwillshedlightonthepassivecontroltheory,andprovidedeepinsightforthetransitionmechanism. Thewell-developedandvalidatedcomputationalalgorithmpresentedherecanbeutilizedforthisresearch.ThecurrentnumericalworksinMAVsexclusivelyconsiderthemembranewingitself.Thefuselageandthepropellerarenotconsidered.Inthefuturework,wewillextendourformerworktoconsiderthefuselageandpropellereects.

PAGE 114

ExperimentsbyWasazketal.[ 111 ]testedanoldMAVmodel,whichisquitedierentfromthecurrentonewestudyhere.ThecurrentMAVmodelfeaturesastreamlizedbodyandanbuilt-inelectricmotor.Thenewmodelissuperiortotheoldoneintheaerodynamicperformance,endurance,control,andmaneuverability.Therefore,acombinedexperimentalandnumericalstudyonthecurrentMAVmodelisneeded. TheMAViesattherelativelylowReynoldsnumberconditions.Complicatedowstructuresoccurattheupperwingsurface,includinglaminarboundarylayerseparation,laminartoturbulentowtransition,andturbulentowreattachment.Sincethereisnocomprehensivetheorytoaccountforallthese,itisimportanttoinvestigatethesephenomenonwithbothexperimentsandnumericalsimulations.Fornumericalcomputations,duetothelowReynoldsnumbercondition,computationswithrenedresolutionsareneededtoaccuratelypredictlaminar-turbulenttransition,detailedowstructures,andtheresultingaerodynamiccoecients. Astothedevelopmentofoptimizationtechniques,oneshouldexplorenewmethodtooptimizethewingshape.Insteadofconsideringasinglepointsolution,weneedtotakeotheranglesofattackinourconsiderationandconsideramulti-objectiveoptimizationproblem.

PAGE 115

Thecomponentsofthenonlinearpartofthestrain-displacementmatrix,BL,canbewrittenasthefollowing, 107

PAGE 117

[1] Alexander,H.,\TheTensileInstabilityofInitiallySphericalBalloons,"Interna-tionalJournalofEngineeringScience,Vol.9,1971,pp.151-162. [2] Anderson,J.D.,Jr.,IntroductiontoFlight,3rded,McGraw-Hill,1989. [3] Batina,J.,\UnsteadyEulerAirfoilSolutionsUsingUnstructuredDynamicMeshes,"AIAAPaper1989-0115. [4] Beaty,M.F.,\TopicsinFiniteElasticity:HyperelasticityofRubber,Elastomers,andBiologicalTissues-withExamples,"AppliedMechanicReview,Vol.40,1987,pp.1699-1734. [5] Bendiksen,O.O.,\ANewApproachtoComputationalAeroelasticity,"1991,AIAAPaper91-0939-CP. [6] Burgee,S.,Giunta,A.A.,Narducci,R.,Watson,L.T.,Grossman,B.,andHaftka,R.T.,\ACoarseGrainedParallelVariable-ComplexityMultidisciplinaryOptimizationParadigm,"TheInternationalJournalofSupercomputerApplica-tionsandHighPerformanceComputing,Vol.10,1996,pp.269-299. [7] Carmichael,B.H.,\LowReynoldsNumberAirfoilSurvey,"NASACR-165803,1981. [8] Chasman,D.,andChakravarthy,S.,\ComputationalandExperimentalStudiesofAsymmetricPitch/PlungeFlapping-TheSecretofBiologicalFlyers,"AIAAPaper2001-0859. [9] Crompton,M.J.,andBarrett,R.V.,\InvestigationoftheSeparationBubbleFormedBehindtheSharpLeadingEdgeofaFlatPlateatIncidence,"JournalofAerospaceEngineering,Vol.214,2000,pp.157-176. [10] Cummings,R.M.,Morton,S.A.,Siegel,S.G.,andBosscher,S.,\NumericalPredictionandWindTunnelExperimentforaPitchingUnmannedCombatAirVehicle,"AIAAPaper2003-0417. [11] DeLaurier,J.D.,\AnAerodynamicModelforFlappingWingFlight,"Aeronau-ticalJournal,1993,pp.125-130. [12] [13] Dickinson,M.H.,Lehmann,F.,andSane,S.P.,\WingRotationandtheAero-dynamicBasisofInsectFlight,"Science,Vol.284,1999,pp.1954-1960. 109

PAGE 118

[14] Ding,H.,Yang,B.,Lou,M.,andFang,H.,\NewNumericalMethodforTwo-DimensionalPartiallyWrinkledMembranes,"AIAAJournal,Vol.41,2003,pp.125-132. [15] Drela,M.,andGiles,M.B.,\Viscous-InviscidAnalysisofTransonicandLowReynoldsNumberAirfoils,"AIAAJournal,Vol.25,1987,pp.1347-1355. [16] Duchon,J.P.,\SplinesMinimizingRotation-InvariantSemi-NormsinSobolevSpaces,"ConstructiveTheoryofFunctionsofSeveralVariables,Oberwolfach1976,editedbySchempp,W.andZeller,K.,Springer-Verlag,Berlin,1977,pp.85-100. [17] Ellington,C.P.,\TheAerodynamicsofHoveringInsectFlight,I.TheQuasi-SteadyAnalysis,"Philos.Trans.R.Soc.LondonSer.A,1984,Vol.35,pp.1-15. [18] Eriksson,L.E.,\GenerationofBoundary-ConrmingGridsAroundWing-BodyCongurationsUsingTransniteInterpolation,"AIAAJournal,Vol.20,1982,pp.1313-1320. [19] Fuentes,C.,He,X.,Carroll,B.,Lian,Y.,andShyy,W.\LowReynoldsNumberFlowsAroundandAirfoilwithaMovableFlap-Part1:Experiments,"AIAAPaper2000-2239. [20] Gad-el-Hak,M.,\Micro-Air-Vehicles:CanTheybeControlledBetter,",tex-titJournalofAircraft,Vol.38,2001,pp.419-429. [21] Gault,D.E.,\AnInvestigationatLowSpeedoftheFlowOveraSimulatedFlatPlateatSmallAnglesofAttackUsingPitotStaticandHot-wireProbes,"1957,NACATN-3876. [22] Gordnier,R.E.,andFithen,R.,\CouplingofaNonlinearFiniteElementStruc-turalMethodwithaNavier-StokesSolver,"ComputersandStructures,Vol.81,2003,pp.75-89. [23] Gordnier,R.E.,andVisbal,M.R.,\DevelopmentofaThree-DimensionalVis-cousAeroelasticSolverforNonlinearPanelFlutter,"JournalofFluidsandStructures,Vol.16,2002,pp.497-527. [24] Grasmeyer,J.M.,andKeennon,M.T.,\DevelopmentoftheBlackWidowMicroAirVehicle,"AIAAPaper2001-0127. [25] Green,A.E.,andAdkins,J.E.,LargeElasticDeformations,TheClarendonPress,Oxford,1960. [26] Farhat,C.,Geuzaine,P.,andBrown,G.,\ApplicationofaThree-eldNonlinearFluid-structureFormulationtothePredictionoftheAeroelasticParametersofanF-16Fighter,"ComputersandFluids,Vol.32,pp.3-29.

PAGE 119

[27] Farhat,C.,Geuzaine,P.,andGrandmont,C.,\TheDiscreteGeometricConser-vationLawandtheNonlinearStabilityofALESchemesfortheSolutionofFlowProblemsonMovingGrids,"JournalofComputationalPhysics,Vol.174,2001,pp.669-694. [28] Farhat,C.,andLesoinne,M.,\TwoEcientStaggeredAlgorithmsfortheSe-rialandParallelSolutionofThree-DimensionalNonlinearTransientAeroelasticProblems,"ComputerMethodsinAppliedMechanicsandEngineering,Vol.182,2000,pp.499-515. [29] Hartwich,P.M.,andAgrawal,S.,\MethodforPerturbingMultiblockPatchedGridsinAeroelasticandDesignOptimizationApplications,"AIAAPaper1997-2038. [30] He,X.,Fuentes,C.,Shyy,W.,Lian,Y.,andCarroll,B.,\ComputationofTransitionalFlowsaroundanAirfoilwithaMovableFlap,"AIAAPaper2000-2240. [31] Herbert,T.,\ParabolizedStabilityEquations,"Annu.Rev.Fluid.Mech.,Vol.29,1997,pp.245-283. [32] Hutchison,M.G.,Unger,E.R.,Mason,W.H.,Grossman,B.,andHaftka,R.T.,\Variable-ComplexityAerodynamicOptimizationofaHigh-SpeedCivilTransportWing,"JournalofAircraft,Vol.31,1994,pp.110-116. [33] Ifju,P.G.,Ettinger,S.,Jenkins,D.,andMartinez,L.,\CompositeMaterialsforMicroAirVehicles,"SAMPEJournal,Vol.37,2001,pp.7-12. [34] Ifju,P.,Jenkins,D.,Ettinger,S.,Lian,Y.,Shyy,W.,andWaszak,R.M.,\Flexible-Wing-BasedMicroAirVehicles,"AIAAPaper2002-0705. [35] Issa,R.I.,\SolutionoftheImplicitlyDiscretisedFluidFlowEquationsbyOperator-Splitting,"JournalofComputationalPhysics,Vol.62,1985,pp.40-65. [36] Jackson,P.S.,\ASimpleModelforElasticTwo-DimensionalElasticSails,"AIAAJournal,Vol.21,1983,pp.153-155. [37] Jackson,P.S.,\TheoryforConicalMembraneWingsofHighAspectRatio,"AIAAJournal,Vol.39,2001,pp.781-786. [38] Jackson,P.S.,andChristie,G.W.,\NumericalAnalysisofThree-DimensionalElasticMembraneWings,"AIAAJournal,Vol.25,1987,pp.676-682. [39] Jameson,A.,\TimeDependentCalculationsUsingMultigrid,withApplicationstoUnsteadyFlowspastAirfoilsandWings,"AIAAPaper1991-1596. [40] Jenkins,C.H.,\NonlinearDynamicResponseofMembranes:StateoftheArt{Update,"AppliedMechanicReview,Vol.49,1996,pp.S41-S48.

PAGE 120

[41] Jenkins,C.H.,andLeonard,J.W.,\NonlinearDynamicResponseofMem-branes:StateoftheArt,"AppliedMechanicReview,Vol.44,1991,pp.319-328. [42] Jones,B.M.,\Stalling,"J.R.Aero.Soc.Vol.38,1938,pp.747-770. [43] Jones,K.D.,andPlatzer,M.F.,\AnExperimentalandNumericalInvestigationofFlapping-WingPropulsion,"AIAAPaper1999-0995. [44] Jones,K.D.,andPlatzer,M.F.,\Flapping-WingPropulsionforaMicroAirVehicle,"AIAAPaper2000-0897. [45] Jones,K.D.,andPlatzer,M.F.,\ExperimentalInvestigationoftheAerody-namicCharacteristicsofFlapping-WingMicroAirVehicles,"AIAAPaper2003-0418. [46] Jones,W.P.,andLaunder,B.E.,\ThePredictionofLaminarizationwithaTwo-EquationModelofTurbulence,"Int.J.HeatMassTrans.,Vol.15,1972,pp.301-314. [47] Kalro,V.,andTezduyar,T.E.,\AParallel3DComputationalMethodforFluid-StructureInteractionsinParachuteSystems,"ComputerMethodsinAppliedMechanicsandEngineering,Vol.190,2000,pp.321-332. [48] Kamakoti,R.,Berg,M.,Ljungqvist,D.,andShyy,W.,\AComputationalStudyforBiologicalFlappingWingFlight,"TransactionofAeronauticalandAstronau-ticalSocietyoftheRepublicofChina,KeynotePaper,Vol.32,2000,pp.265-279. [49] Kamakoti,R.,Lian,Y.,Regisford,S.,Kurdila,A.,andShyy,W.,\Compu-tationalAeroelasticityUsingaPressure-basedSolver,"ComputerModelinginEngineering&Sciences,Vol.3,2002,pp.773-790. [50] Katz,J.,Low-SpeedAerodynamics:FromWingTheorytoPanelMethods,Mc-Graw-Hill,SanFrancisco,CA,1979. [51] Knill,D.L.,Giunta,A.A.,Baker,C.A.,Grossman,B.,Mason,W.H.,Haftka,R.T.,andWatson,L.T.,\ResponseSurfaceMethodsCombiningLinearandEulerAerodynamicsforSupersonicTransportDesign,"JournalofAircraft,Vol.36,1999,pp.75-86. [52] Launder,B.E.,andSpalding,D.B.,\TheNumericalComputationofTurbulentFlows,"Comp.Meth.Appl.MechEngg.,3,1974,pp.269-289. [53] Lian,Y.,andShyy,W.,\Three-DimensionalFluid-StructureInteractionsofaMembraneWingforMicroAirVehicleApplications,"AIAAPaper2003-1726. [54] Lian,Y.,Shyy,W.,andHaftka,R.,\ShapeOptimizationofaMembraneWingforMicroAirVehicles,"AIAAPaper2003-0106.AcceptedforpublicationatAIAAJournal.

PAGE 121

[55] Lian,Y.,Shyy,W.,Ifju,P.,andVerron,E.,\AComputationalModelforCou-pledMembrane-FluidDynamics,"AIAAPaper2002-2972.Acceptedforpubli-cationatAIAAJournal. [56] Lian,Y.,Steen,J.,Trygg-Wilander,M.,andShyy,W.,\LowReynoldsNumberTurbulentFlowsaroundaDynamicallyShapedAirfoil,"ComputersandFluids,Vol.32,2003,pp.287-303. [57] Lighthill,M.J.,\HydrodynamicsofAquaticAnimalPropulsion,"Ann.Rev.FluidMech.,1969,pp.413-445. [58] Liu,F.,Zhu,Y.,Tsai,H.M.,andWong,A.S.F.,\CalculationofWingFlutterbyaCoupledFluid-StructureMethod,"JournalofAircraft,Vol.38,2001,pp.334-342. [59] Liu,H.,andKawachi,K.,\ANumericalStudyofInsectFlight,"Vol.J.Comput.Phys.,146,1998,pp.124-156. [60] Liu,H.T.,\UnsteadyAerodynamicsofaWortmannWingatLowReynoldsNumber,"JournalofAircraft,1992,Vol.29,pp.532-539. [61] Lissaman,P.B.S.,\Low-Reynolds-numberAirfoils,"Annu.Rev.FluidMech.Vol.15,1983,pp.223-239. [62] Lorillu,O.,andHureau,J.,\NumericalandExperimentalAnalysisofTwo-DimensionalSeparatedFlowsoveraFlexibleSail,"J.Fluid.Mech.,Vol.466,pp.319-341. [63] Mason,B.H.,Haftka,R.T.,Johnson,E.R.,andFarley,G.L.,\VariableCom-plexityDesignofCompositeFuselageFramesbyResponseSurfaceTechniques,"ThinWallStructures,Vol.32,1998,pp.235-261. [64] Mooney,M.,\ATheoryofLargeElasticDeformation,"JournalofAppliedPhysics,Vol.11,1940,pp.582-592. [65] Mueller,T.J.(ed.),ProceedingsoftheConferenceonLowReynoldsNumberAerodynamics,Univ.ofNotreDame,NotreDame,IN,1989. [66] Mueller,T.J.(ed.),ProceedingsoftheConferenceonFixed,FlappingandRotaryWingVehiclesatVeryLowReynoldsNumbers,Univ.ofNotreDame,NotreDame,IN,2000. [67] Mueller,T.J.,andDeLaurierJ.D.,\AerodynamicsofSmallVehicles,"AnnualReviewofFluidMechanics,Vol.35,2003,pp.89-111. [68] Nelsen,J.N.,\TheoryofFlexibleAerodynamicsSurfaces,"JournalofAppliedMechanics,Vol.30,1963,pp.435-442.

PAGE 122

[69] Newman,B.G.,\AerodynamicsTheoryforMembraneandSails,"ProgressinAerospaceSciences,Vol.24,1987,pp.1-27. [70] Oden,J.T.,andSato,T.,\FiniteStrainsandDisplacementsofElasticMem-branesbytheFiniteElementMethod,"InternationalJournalofSolidsandStructures,3,1967,pp.471-488. [71] Oliveira,P.J.,andIssa,R.I.,\AnImplicitPISOAlgorithmfortheComputationofBuoyancy-DrivenFlows,"NumericalHeatTransfer,PartB,Vol.40,2001,pp.473-493. [72] Patankar,S.V.,andSpalding,D.B.,\ACalculationProcedureforHeat,MassandMomentumTransferinThree-DimensionalParabolicFlows,"Int.J.HeatMassTransf.,Vol.15,1972,pp.1787-1806. [73] Patel,V.C.,Rodi,W.,andScheuerer,G.,\TurbulenceModelsforNear-WallandLowReynoldsNumberFlows,"AIAAJournal,23,1985,pp.1308-1319. [74] Pulliam,T.,\TimeAccuracyandtheUseofImplicitMethods,"1993,AIAAPaper93-3360-CP. [75] Reed,H.L.,andSaric,W.S.,\LinearStabilityTheoryAppliedtoBoundaryLayers,"Annu.Rev.Fluid.Mech.,Vol.28,1996,pp.389-428. [76] Rempfer,D.,\Low-DimensionalModelingandNumericalSimulationofTransi-tioninSimpleShearFlows,"AnnualReviewofFluidMechanics,Vol.35,2003,pp.229-265. [77] Reuther,J.,Alonso,J.J.,Rimlinger,M.J.,andJameson,A.,\AerodynamicShapeOptimizationofSupersonicAircraftCongurationviaanAdjointFormu-lationonDistributedMemoryParallelComputers,"Comput.Fluids,28,1999,pp.675-700. [78] Robinson,B.A.,Batina,J.T.,andYang,H.T.Y.,\AeroelasticAnalysisofWingsUsingtheEulerEquationswithaDeformingMesh,"JournalofAircraft,Vol.28,1991,pp.781-788. [79] Rumsey,C.L.,Sanetrik,M.D.,Biedron,R.T.,Melson,N.D.,andParlette,E.B.,\EciencyandAccuracyofTime-AccurateTurbulentNavier-StokesCom-putations,"ComputersandFluids,Vol.25,1996,pp.217-236. [80] Schuster,D.M.,Liu,D.D.,andHuttsell,L.J.,\ComputationalAeroelasticity:Success,Progress,Challenge,"AIAAPaper2003-1725. [81] Schuster,D.,Vadyak,J.,andAtta,E.,\StaticAeroelasticAnalysisofFighterAircraftUsingaThree-DimensionalNavier-StokesAlgorithm,"J.Aircraft,Vol.27,1990,pp.820-825.

PAGE 123

[82] Shyy,W.,ComputationalModelingforFluidFlowandInterfacialTransport,Elsevier,Amsterdam,TheNetherlands,1994. [83] Shyy,W.,Berg,M.,andLjungqvist,D.,\FlappingandFlexibleWingsforBi-ologicalandMicroVehicles,"ProcessinAerospaceSciences,Vol.35,1999,pp.455-506. [84] Shyy,W.,Jenkins,D.A.,andSmith,R.W.,\StudyofAdaptiveShapeAirfoilsatLowReynoldsNumberinOscillatoryFlow,"AIAAJournal,Vol.35,1997,pp.1545-1548. [85] Shyy,W.,andMittal,R.,\SolutionMethodsfortheIncompressibleNavier-StokesEquations,"inR.Johnson(ed.)HandbookofFluidDynamics,CRCPress,BocaRaton,FL,1998,pp.31-1to31-33. [86] Shyy,W.,andSmith,R.W.,\ComputationofLaminarFlowandFlexibleStruc-tureInteraction,"inM.HafezandK.Oshima(eds.)ComputationalFluidDy-namicsReview,JohnWiley&Sons,1995,pp.777-796. [87] Shyy,W.,Thakur,S.S.,Ouyang,H.,Liu,J.,andBlosch,E.,ComputationalTechniquesforComplexTransportPhenomena,CambridgeUniversityPress,1997. [88] Shyy,W.,Udaykumar,H.S.,Rao,M.M.,andSmith,R.W.,ComputationalFluidDynamicswithMovingBoundaries,1996,Taylor&Francis,Washington,DC. [89] Smith,M.J.,\SimulatingMothWingAerodynamics:TowardstheDevelopmentofFlapping-WingTechnology,"AIAAJournal,Vol.34,1996,pp.1348-1355. [90] Smith,M.J.\ComputationalConsiderationsofanEuler/Navier-StokesAeroe-lasticMethodforaHoveringRotor,"JournalofAircraft,Vol.33,1996,pp.429-434. [91] Smith,M.J.,Hodges,D.H.,andCesnik,C.E.S.,\EvaluationofComputationalAlgorithmsSuitableforFluid-StructureInteractions,"JournalofAircraft,Vol.37,2000,pp.282-294. [92] Smith,R.W.,andShyy,W.,\ComputationalModelofFlexibleMembraneWingsinSteadyLaminarFlow,"AIAAJournal,Vol.33,1995,pp.1769-1777. [93] Smith,R.W.,andShyy,W.,\IncrementalPotentialFlowBasedMembraneWingElement,"AIAAJournal,Vol.35,1997,pp.782-788. [94] Subbaraj,K.,andDokainish,M.A.,\ASurveyofDirectTime-IntegrationMeth-odsinComputationalStructuralDynamics-II.ImplicitMethods,"ComputersandStructures,Vol.32,1989,pp.1387-1401.

PAGE 124

[95] Tani,I.,\Low-SpeedFlowsInvolvingBubbleSeparations,"ProgressinAeronau-ticalScience,Vol.5,editedbyD.Kuchemannandl.H.G.Sterne,Pergamon,NewYork,1964,pp.70-103. [96] Templin,R.J.,\TheSpectrumofAnimalFlight:InsectstoPterosaurs,"ProgressinAerospaceSciences,Vol.36,2000,pp.393-436. [97] Torres,G.E.,andMueller,T.J.,\AerodynamicCharacteristicsofLowAspectRatioWingsatLowReynoldsNumber."InFixedandFlappingWingAerody-namicsforMicroAirVehicleApplications,ed.TJMueller,Vol.195,2001,pp.191-213.Reston,VA:AIAA. [98] Thakur,S.,Wright,J.,andShyy,W.,STREAM:AComputationalFluidDynam-icsandHeatTransferNavier-StokesSolver.TheoryandApplications,StreamlineNumerics,Inc.,andComputationalThermo-FluidsLaboratory,DepartmentofMechanicalandAerospaceEngineeringTechnicalReport,Gainesville,Florida,2002. [99] Thaokar,R.M.,andKumaran,V.,\StabilityofFluidFlowPastaMembrane",J.FluidMech.,Vol.472,2002,pp.29-50. [100] Thomas,P.D.,andLombard,C.K.,\GeometricConservationLawandItsApplicationtoFlowComputationsonMovingGrids,"AIAAJournal,Vol.17,1979,pp.1030-1037. [101] Thwaites,B.,\AerodynamicsTheoryofSails-PartI,"ProceedingsoftheRoyalSociety,A261,1961,pp.402-442. [102] VanDoormaal,J.P.,andRaithby,G.D.,\EnhancementsoftheSIMPLEMethodforPredictingIncompressibleFluidFlows,"NumericalHeatTransfer,Vol.67,1985,pp.147-163. [103] Versteeg,H.K.,andMalalasekera,W.,AnIntroductiontoComputationalFluidDynamics,LongmanScientic&Technical,1995. [104] Verron,E.,Khayat,R.E.,Derdouri,A.,andPeseux,B.,\DynamicInationofHyperelasticSphericalMembranes,"JournalofRheology,Vol.43,1999,pp.1083-1097. [105] Verron,E.,Marckmann,G.,andPeseux,B.,\DynamicInationofNon-linearElasticandViscoelasticRubber-likeMembranes,"Int.J.Numer.Mech.Engng.,Vol.50,2001,pp.1233-1251. [106] Vest,M.S.,andKata,J.,\UnsteadyAerodynamicModelofFlappingWings,"AIAAJournal,Vol.34,1996,pp.1435-1440. [107] Vest,M.S.,andKatz,J.,\AerodynamicStudyofaFlappingWingMicro-UAV,"AIAA1999-0994.

PAGE 125

[108] Voelz,K.,\ProlundLuftriebeinesSegels,"ZAMM,Vol.30,1950,pp.301-317. [109] Walker,J.A.\FunctionalMorphologyandVirtualModels:PhysicalCon-straintsontheDesignofOscillatingWing,Fins,Legs,andFeetatIntermediateReynoldsNumbers,"Integ.andComp.Biol.,Vol.42,2002,pp.232-242. [110] Wang,Z.J.,\VortexSheddingandFrequencySelectioninFlappingFlight,"JournalofFluidMechanics,Vol.410,2000,pp.323-341. [111] Waszak,R.M.,Jenkins,N.L.,andIfju,P.,\StabilityandControlPropertiesofanAeroelasticFixedWingMicroAerialVehicle,"AIAAPaper2001-4005. [112] Weis-Fogh,T.,\QuickEstimatesofFlightFitnessinHoveringAnimals,Includ-ingNovelMechanismsforLiftProduction,"JournalofExperimentalBiology,Vol.59,1973,pp.169-230. [113] Wilson,E.L.,Farhoomand,I.,andBathe,K,J.,\NonlinearDynamicAnalysisofComplexStructures,"EarthquakeEngineeringandStructuralDynamics,Vol.1,1973,pp.241-252. [114] Young,A.D.,andHorton,H.P.,\SomeResultsofInvestigationofSeparationBubbles,"AGARDCP,4,1966,pp.779-811. [115] Zwaan,R.J.,andPrananta,B.B.,\Fluid/structureInteractioninNumericalAeroelasticSimulation,"InternationalJournalofNon-LinearMechanics,Vol.37,2002,pp.987-1002.

PAGE 126

YongshengLianwasborninthesmallvillageofZiBo,ShandongProvince,inthePeople'sRepublicofChina.Heearnedhisbachelor'sdegreefromShandongUniversity,majoringinmathematics.BeforehejoinedtheUniversityofFlorida,heearnedtwomaster'sdegrees(bothinmathematics):onefromtheChineseAcademyofSciences;andtheotherfromtheHongKongUniversityofScienceandTechnology.In1999hejoinedtheUniversityofFloridaasagraduatestudent.HewastherecipientoftheAlumniFellowshipattheUniversityofFlorida.Yongshengenjoyssportsevents.HeisafanofGatorsports. 118