A viscous flow based membrane wing model

MISSING IMAGE

Material Information

Title:
A viscous flow based membrane wing model
Physical Description:
vi, 95 leaves : ill. ; 29 cm.
Language:
English
Creator:
Smith, Richard W., 1956-
Publication Date:

Subjects

Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1994.
Bibliography:
Includes bibliographical references (leaves 92-94).
Statement of Responsibility:
by Richard W. Smith.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 002045436
notis - AKN3360
oclc - 33393748
System ID:
AA00003639:00001

Full Text









A VISCOUS FLOW BASED MEMBRANE WING MODEL


By

RICHARD W. SMITH

















A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


1994













ACKNOWLEDGEMENTS

I would first like to thank my advisor, Professor Wei Shyy, for allowing me to follow

my own research interests and for providing a productive working environment for all his

students. Without his support and tolerant attitude this thesis could not have been written.

It is with the same sense of gratitude and good fortune that I acknowledge here the support

and encouragement my family has continuously provided over the years. I would also like

to thank all my colleagues in the AeMES Department at the University of Florida. In

particular, Dr. Jeffrey Wright and Dr. Siddharth Thakur were very generous with their time.

Their willingness to share their knowledge and experience made my work much easier.

Finally, I would like to thank Eglin Air Force Base and the Waterways Experiment Station

Supercomputing Center for the use of their Cray computers.












TABLE OF CONTENTS


page

ACKNOWLEDGEMENTS ................................ ii

ABSTRACT ............................................ vi

CHAPTERS

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

1.1 Overview of the Sail Analysis Problem .................. 1
1.2 Review of the Literature ............................. 2
1.2.1 Two-Dimensional Potential Flow Based Models ....... 2
1.2.2 Two-Dimensional Viscous Flow Based Models ........ 4
1.2.3 Three-Dimensional Potential Flow Based Models ...... 5
1.3 Focus of the Present Work ............................ 6

2 GOVERNING EQUATIONS ............................... 9

2.1 Membrane Equilibrium ............................. 9
2.2 Fluid Dynamic Conservation Laws ..................... 11
2.3 Nondimensionalization of the Governing Equations ........ 12

3 NUMERICAL METHOD ................................. 17

3.1 Membrane Equilibrium ............................ 17
3.2 Fluid Dynamic Conservation Laws .................... 17
3.3 Consistent Implementation of the Continuity Equation ..... 26
3.4 The Aeroelastic Computational Procedure ............... 27
3.5 A Potential Flow Model for Thin Wings ................ 29

4 ELEMENTARY TEST CASES ............................ 31

4.1 Rotated Channel Flow ............................... 31
4.2 Uniform Flow Using a Moving Grid .................... 33
4.3 Elastic Membrane Under a Uniform Pressure Load ........ 34

5 MEMBRANE WINGS IN STEADY FLOW .................. 37

5.1 Rigid Wing in a Viscous Fluid ......................... 37
5.1.1 Effect of Outer Boundary Location ................. 37
5.1.2 Effect of Grid Refinement ........................ 39








5.2 Flexible Membrane Wing in a Viscous Fluid ............. 45
5.2.1 Elastic Membrane Case .......................... 46
5.2.2 Constant Tension Membrane Case .................. 50
5.2.3 Inextensible Membrane Case ...................... 51

6 MEMBRANE WINGS IN UNSTEADY FLOW ................ 56

6.1 Constant Tension Membrane Case ...................... 58
6.2 Elastic Membrane Case .............................. 59
6.3 Inextensible Membrane Case .......................... 66

7 A THREE-DIMENSIONAL MEMBRANE WING MODEL ..... 70

7.1 The Elastic Membrane Problem in Three-Dimensions ...... 70
7.2 The Aerodynamic Problem in Three-Dimensions .......... 77
7.3 The Aeroelastic Problem ............................. 81
7.4 Three-Dimensional Test Cases ........................ 82
7.5 Application to an Elliptical Planform Elastic Wing ........ 86

8 SUMMARY AND CONCLUSIONS ......................... 90

REFERENCES .......................................... 92

BIOGRAPHICAL SKETCH............................... 95














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

A VISCOUS FLOW BASED MEMBRANE WING MODEL

By

Richard W. Smith

December 1994

Chairman: Wei Shyy
Major Department: Aerospace Engineering, Mechanics and Engineering Science

In this work a numerical model simulating the aeroelastic characteristics of a flexible

membrane wing is presented. The use of the Navier-Stokes equations as the fluid dynamic

model in the present membrane wing theory is a substantial departure from previous work

on sail aerodynamics which has, almost universally, adopted a potential based description

of the flow field. The two-dimensional, viscous aeroelastic problem is nondimensionalized

and a set of six basic dimensionless parameters is derived which govern the physical

problem. An additional parameter, the frequency ratio, is proposed as a meaningful

parameter for characterizing the harmonically driven unsteady aeroelastic problem. A

numerical procedure is developed for the solution of the coupled aeroelastic problem and

is shown to yield results that are in agreement with available analytic solutions for several

appropriate limiting cases. The numerical procedure is also shown to satisfy certain

identities as dictated by the fundamental fluid dynamic conservation laws. The role of

viscosity in membrane wing aerodynamics is investigated using the numerical model for

both steady and unsteady flows. These investigations are facilitated by distinguishing three

distinct classes of problems which are associated with limiting cases of the dimensionless

parameter set. The aerodynamic characteristics at Reynolds numbers between 103 and 104







are shown to differ substantially from those predicted by a potential based membrane wing

theory. The role of viscosity is shown to be preeminent in the harmonically forced unsteady

flow about a membrane wing. In this case in the influence of viscosity is enhanced since the

acceleration and deceleration of the freestream velocity strongly influences the separation

and reattachment of the flow. The periodic appearance and collapse of recirculation zones,

along with an attendant adjustment in the membrane configuration, result in an aeroelastic

response which may not be characterized as a simple harmonic response at the freestream

forcing frequency. A three-dimensional potential based membrane wing element is also

derived and is shown to yield results that, for several limiting cases, are in agreement with

available analytic solutions.













CHAPTER 1
INTRODUCTION

1.1 Overview of the Sail Analysis Problem

Marchaj, in his second book on the science of sailing (Marchaj 1979), begins the

section on sail design with the following comments.



Despite the fact that mathematics, computers and wind tunnel testing are
playing an increasing part in the designing of sails, sailmaking as well as sail tuning
are still strongholds of art based on a hit-or-miss technique rather than on science..
.. After all, unlike the aeroplane wing, which can be regarded as a rigid structure whose
shape is unaffected by variation in incidence and speed, sail shape is a function of both;
in which the shape of the sail affects the pressure distribution and vice versa, in a rather
unpredictable manner. (p. 500)


These comments by Marchaj suggest two things concerning the analysis of marine

sails: first, the behavior and performance of a sail is governed by the aeroelastic interaction

of a fluid dynamic field and a deformable surface; and secondly, as a result of this interaction,

as well as the presence of other complexities, analytic and computational methods which

have enjoyed considerable success in the design of aircraft have not as yet proven useful to

sail designers. It is likely that the second observation concerning the usefulness of analytic

methods to sail designers will not always hold true.

As an illustration of the potential usefulness of computational methods in the analysis

of sails consider the membrane wing of Fig. 1, which is shown operating near a boundary

such as the surface of the sea. Since the sail is located very near the boundary it is immersed

in the shear layer adjacent to the free surface. The kinematic inflow conditions to the sail are

consequently nonuniform and may also be unsteady due to the presence of wind gusts.

Furthermore, since the sail is not stationary but rather has some forward velocity, the relative







inflow velocity far upstream of the sail will vary in both magnitude and direction along the

sail span as well as possibly varying in time. Computational methods provide a method of

simulating such a complex flow environment which would otherwise be nearly impossible

to reproduce experimentally.




membrane wing




inflow velocity /








free surface V


Figure 1.1. Schematic of a membrane wing of finite span operating near a free surface.


1.2 Review of the Literature


1.2.1 Two-Dimensional Potential Flow Based Models

The vast majority of the published works related to membrane wing aerodynamics

have made several simplifying approximations concerning both the elastic characteristics of

the membrane itself as well as the nature of the surrounding flow field. Perhaps the most

significant of these simplifying assumptions is that the fluid dynamics can be adequately

described by a potential based model of the flow field. In addition to the almost universally

adopted potential flow assumption, the additional approximations associated with thin

airfoil theory small camber and incidence angle are also often made and the membrane

itself is generally considered to be inextensible. A comprehensive review of the work







published prior to 1987 related to membrane wing aerodynamics is given by Newman

(1987).

The analysis of membrane wings begins with the historical works of Voelz (1950),

Thwaites (1961) and Nielsen (1963). These works considered the steady, two dimensional,

irrotational flow over an inextensible membrane with slack. As a consequence of the

inextensible assumption and the additional assumption of small camber and incidence angle

the membrane wing boundary value problem is linearized and may be expressed compactly

in nondimensional integral equation form as
1 d2(y/a)
SCT d2 = d(y/a)
1 d 2(5 x)d (1.1)
2 27t( x) dx
0

where y(x) defines the membrane profile as a function the x coordinate, a is the flow

incidence angle and CT is the tension coefficient. Equation (1.1) has been referred to as the

'Thwaites sail equation' by Chambers (1966) and simply as the 'sail equation' by Newman

(1987) and Greenhalgh et al. (1984). This equation, together with a dimensionless geometric

parameter e which specifies the excess length of the membrane, completely defines the

linearized theory of an inextensible membrane wing in an steady, inviscid flow field.

Different analytic and numerical procedures have been applied to the sail equation

in order to determine the membrane shape, aerodynamic properties, and membrane tension

in terms of the angle of attack and excess length. In particular, Thwaites (1961) obtained

eigensolutions of the sail equation which are associated with a wing at the ideal (singularity

free) angle of incidence. Nielsen (1963) obtained solutions to the same equation using a

Fourier series approach which is valid for wings at angles of incidence other than the ideal

angle. Other more recent but similar works are those by Chambers (1966), Vanden-Broeck

and Keller (1981), Greenhalgh et al. (1984) and Sugimoto and Sato (1988).

Various extensions of the linear theory have appeared in the literature over the years.

Vanden-Broeck (1982) as well as Murai and Maruyama (1980) developed nonlinear theories







valid for large camber and incidence angle. The effect of elasticity has been included in the

membrane wing theories of Jackson (1983) and Sneyd (1984) and the effects of membrane

porosity have been investigated by Murata and Tanaka (1989). In a paper by de Matteis and

de Socoi (1986) experimentally determined separation points were used to modify the lifting

potential flow problem in an attempt to model flow separation near the trailing edge. The

effect of elasticity and porosity were also considered in this work.

Comparisons of the various potential flow based membrane wing theories with

experimental data have been reported by several authors including Greenhalgh et al. (1984),

Sugimoto and Sato (1988) and Newman and Low (1984). In general, there has been

considerable discrepancy between the measurements made by the different authors which

have all been in the turbulent flow regime at Reynolds numbers between 105 and 106

(Jackson and Fiddes 1993). As a result of the discrepancies in the reported data primarily

due to differences in Reynolds number and experimental procedure the agreement between

the potential based membrane theories and the data has been mixed. In particular, the

measured lift is in fair agreement with the predicted value when the excess length ratio is less

than .01 and the angle of attack of attack is less than 5. However, even for this restricted

range of values the measured tension is significantly less than predicted by theory.

Furthermore, for larger excess lengths and incidence angles, both lift and tension are poorly

predicted by the theory. Flow visualization studies indicate the main reason for the

disagreement is the existence of a thick boundary layer or region of separated flow on the

membrane, typically near the trailing edge. It has been noted by several authors that the

presence of viscous effects such as thick boundary layers and separation regions will

overshadow any implications associated with the linearizing approximations made by thin

wing theory (Newman 1987).

1.2.2 Two-Dimensional Viscous Flow Based Models

Although the importance of viscous effects on membrane wing aerodynamics has

been recognized for quite some time (Nielsen 1963, Newman and Low 1984), very few







works have been published which address the issue. Jackson and Fiddes (1993) coupled a

boundary layer calculation to a potential flow model which was then combined with an

inextensible membrane model. The first use of the Navier-Stokes equations as the fluid

dynamic model in a membrane wing theory appears to be the work of Smith and Shyy (1993).

In this work an incremental, material coordinate formulation of the elastic membrane

problem was coupled with a body-fitted, finite volume formulation of the Navier-Stokes

equations. Results from the viscous flow based membrane wing model were compared with

a potential flow based membrane wing theory. Although the results were in qualitative

agreement with the flow visualization studies given by Newman and Low (1984) a

meaningful quantitative comparison with the experimental data was not possible due to the

assumption of laminar flow in the model.

1.2.3 Three-Dimensional Potential Flow Based Models

The development of a three-dimensional aeroelastic membrane wing model

introduces several complicating factors in the analysis which do not have to be considered

in the two-dimensional analysis. In contrast to the two-dimensional case where the

inextensible approximation is often valid, the effects of elasticity are preeminent in

three-dimensional wings where spanwise twist and trailing edge deflection are strongly

influenced by the elastic properties of the membrane. Also, even within the context of

inviscid flow, the tension in a three-dimensional membrane is no longer a single constant

value but is best described by a state of biaxial tension directed along lines of principal stress

(Jackson 1985). Furthermore, if one of the principal tensions vanishes the membrane will

compress and wrinkle which further complicates the analysis. Finally, the effect of viscosity,

which is known to affect strongly the aerodynamics of two-dimensional membrane wings

under many conditions, must also play an important role in the behavior of finite span

membrane wings.

Several works have appeared in the literature which have considered various

idealizations of the general three-dimensional membrane wing problem. The







two-dimensional theory described in the previous section was combined in a stripwise

fashion with lifting line theory by Nielsen (1963), Sneyd (1984) and Ormiston (1971). Murai

and Maruyama (1982) developed a model for a rectangular planform wing assuming the

tension in the chordwise direction to be zero. Holla et al. (1984) also investigated a

rectangular membrane wing with fixed edges assuming a state of uniform biaxial tension to

exist in the membrane. Kroo (1986) developed a numerical method for arbitrary planforms

and edge conditions which was based on a strain energy minimization formulation of the

membrane problem and a vortex lattice treatment of the potential flow problem. Kroo (1986)

used this model to investigate the stability of hang gliders. A summary of moderate aspect

ratio membrane wing theories is given by Newman (1987).

The numerical model developed by Jackson (1985) and Jackson and Christie (1987)

is probably the most general three-dimensional model presented to date that has been applied

to the aeroelastic analysis of marine sails. Jackson's model combines a vortex lattice

treatment of the lifting, potential flow problem with a Lagrangian, finite element based

formulation of the elastic membrane problem. An iterative procedure is then used to solve

the coupled aeroelastic boundary value problem. Results for rectangular and triangular

planform sails with several different edge boundary conditions are presented in the 1985 and

1987 papers. The issue of membrane wrinkling is also discussed in these papers.

1.3 Focus of the Present Work

The primary focus of the present work is to develop a computational model of a

two-dimensional membrane wing which fully accounts for the effect of viscosity in the fluid.

The incorporation of a viscous fluid dynamic model in the analysis distinguishes the present

work from previous analytic and computational work on sail aerodynamics in which the

fluid dynamics were characterized by a potential based description. Many authors have noted

that flow separation is often observed in flow visualization studies with membrane wings and

several have gone on to say that the existence of flow separation may be the primary cause







for the disagreement between the experimental data and existing membrane wing theories

(Newman and Low 1984, Newman 1987, Sugimoto and Sato 1988).

In Chapter 2 the governing fluid dynamic conservation laws are presented along with

the spatial coordinate based equations of equilibrium for a flexible, linearly elastic

membrane. The boundary conditions on fluid pressure, membrane geometry and fluid

velocity which couple the membrane equilibrium equations to the fluid dynamic

conservation laws are also presented. The governing equations and boundary conditions are

then nondimensionalized to determine the dimensionless parameters that characterize the

unsteady, viscous, aeroelastic membrane wing problem.

In Chapter 3 the discrete form of the governing equations is given. The

Navier-Stokes equations are first transformed from Cartesian to general curvilinear

coordinates and a pressure based control volume formulation is adopted as the basic point

of departure for the solution algorithm. Special attention is given during the development

of the numerical method to the consistent transformation of the components of the fluid

velocity vector from Cartesian to curvilinear coordinates. Special care is also given to the

numerical treatment of the time derivative term in the transformed continuity equation to

insure a consistent, source free implementation of mass conservation. The membrane

equilibrium equations are implemented in a straightforward way using conventional finite

difference approximations.

The numerical method developed in Chapter 3 is applied in Chapter 4 to several

elementary test problems for which analytic solutions are known. The test cases presented

for the Navier-Stokes algorithm are chosen to illustrate the need for special care in certain

aspects of the implementation. The test case given for the membrane equilibrium algorithm

is included simply as a check on the accuracy of the computation.

In Chapter 5 the computational procedure is applied to three separate classes of

aeroelastostatic membrane wing problems. These three classes of problems arise naturally

from the nondimensionalization of the governing equations and are referred to as the







constant tension case, the elastic case, and the inextensible case. Results from the viscous

flow based model are compared with results using a potential flow model of the fluid

dynamics. A very limited comparison of the computational results with selected

experimental data is also made in this section. A comprehensive comparison with

experimental data is precluded in the present work by the assumption of laminar flow since

the experimental data were obtained at Reynolds numbers in the turbulent regime.

In Chapter 6 the numerical method is applied to membrane wings subjected to a

harmonically varying freestream velocity. The three classes of wings are investigated and

the time dependent aerodynamic properties of the wing are compared with the steady state

values under comparable conditions.

In Chapter 7 membrane wings of finite span are considered and an incremental,

velocity stepping, finite element based formulation is presented. The fluid dynamics are

approximated by a potential based description of the flow, and a material coordinate based

description of the membrane kinematics is adopted. A three-dimensional aeroelastostatic

membrane wing element is derived and tested on several limiting cases of wing geometry,

pressure distribution and material stiffness before the algorithm is applied to an elliptical

planform, edge constrained membrane wing of moderate aspect ratio.

Finally, in Chapter 8, the present work is summarized and several observations are

made concerning the role of viscosity in the steady and unsteady flow about a membrane

wing. Suggestions for future work are made in this chapter with particular emphasis on

incorporating the effect of turbulence in the membrane wing model. The importance of

comparing the computational predictions with experimental data is also discussed in this

chapter.














CHAPTER 2
GOVERNING EQUATIONS

2.1 Membrane Equilibrium

In this section the general equilibrium equations are presented for a two-dimensional

elastic membrane subjected to both normal and shearing stresses. The membrane is assumed

to be massless and the equilibrium conditions are stated in terms of the instantaneous spatial

Cartesian coordinates xi and the body-fitted curvilinear coordinates i. The basic

formulation is essentially identical to many previously published works such as de Matteis

and de Socoi (1986) and Sneyd (1984). Index notation with the summation convention is

used throughout the following chapters.







low pointE
p point W point P t r '

X2 LL T
T membrane free body




c

Figure 2.1. End constrained elastic membrane


Figure 2.1 illustrates an elastic membrane restrained at the leading and trailing edges

subjected to both fluid dynamic pressure and shear stress, p and T respectively. Imposing

equilibrium in the normal and tangential directions requires








d2X2 I + = (2.1)


dT
d (2.2)


where T is the membrane tension. The net pressure and shear stress acting on a segment of

the membrane are given respectively by

Ap = p- p+ (2.3)

= T- + z+ (2.4)

where the superscript indicates the value at the upper and lower surface of the membrane as

shown in the figure. If the membrane material is assumed to be linearly elastic the nominal

membrane tension 7 may be written in terms of the nominal membrane strain 5 as

T = (S + E)h (2.5)

where S is the membrane prestress, E is the elastic modulus, and h is the membrane

thickness. The nominal membrane strain is given by

6- = L L (2.6)
Lo

where LO is the prestrained length of the membrane and L is the length of the membrane after

deformation which may be expressed in terms of the spatial coordinates xi as

c2
L = 1 + dxl (2.7)
0

where c is the chord length.

At the leading and trailing edges of the membrane the following boundary

conditions are imposed on Eq. (2.1)


x2 = 0 at x, = 0,c


(2.8)







2.2 Fluid Dynamic Conservation Laws

The governing conservation laws for unsteady, laminar, incompressible flow are the
Navier-Stokes equations which may be written in two-dimensional Cartesian coordinates as

a a ap (2.9)
( vl) + t- ( eviv) = (/ Vlj) x (2.9)
SaxJ

a a a ap
S) ( ?vz )+ + ( Vj2 ( V a x2 (2.10)


a + a (evj) = 0 (2.11)

where vi are the Cartesian components of the fluid velocity vector, Q is the fluid density, jt
is the fluid viscosity, t is time, and p is the fluid pressure.
When new independent variables 1\ and k2 are introduced Eqs. (2.9) through (2.11)
change according to the general transformation 1= 1 (xi, x2, t) and k2 = k2 (XI, X2, t).
Equations (2.9) through (2.11) can be rewritten as follows where the subscript i,j indicates
the partial derivative of the i Cartesian component of velocity or position with respect to the
general curvilinear coordinate j.

-( JQvi ) + -a(eoVj,) (_ ql, q21,2 )) +

S-( q3v1,2 21,1 )) (2,2 P) + 2,l P) (2.12)


a( Jv2 ) + -( Vjv2) l( ( qlv2,1 q2v2,2 )) +

( ( q32,2 q2v2,1 )) + (1,2 P) (xl,1 P) (2.13)


( J ) + ( PV) = 0 (2.14)

where the contravariant velocity components are given by

V1 = (vI x-) x2,2 (2 X2)X1,2 (2.15)


V2 = (v2 x2) X1,1 ( xl) x2,1


(2.16)




12


and

ql = (xl,)2 + (x2,2)2 (2.17)

q2 = Xl,lx,2 + x2,1x2,2 (2.18)

q3 = (X,1)2 + (2,1)2 (2.19)

and the Jacobian of the transformation is defined as

J = Xl,lX2,2 Xl,2X2,1 (2.20)

The following kinematic boundary conditions are imposed on the Cartesian

components of the fluid velocity vector at the membrane surface

vi = xi (2.21)

and, in the present work, the fluid velocity far upstream of the membrane is chosen to be a

harmonic function of time given by

v = v, cosa(1 + sinwt ) (2.22)


V2 = v. sina(1 + fsinot ) (2.23)

where 0o is the circular frequency of the freestream velocity and 3 is a parameter defining

the magnitude of the oscillation.

2.3 Nondimensionalization of the Governing Equations

The aeroelastic boundary value problem can be written in nondimensional form after

introducing the following dimensionless variables

v = (2.24)

2 = (2.25)


= tA (2.26)
t = t(O)


S-X1
Xi = 7


(2.27)








x2 =

A Al
Ap = -A

A T
OSh


A T
T =Eh


(2.28)


(2.29)


(2.30)


(2.31)


where either Eq. (2.30) or Eq. (2.31) is used to nondimensionalize the membrane tension

depending on whether the tension is dominated by pretension or by elastic strain.

Substituting these variables into Eq. (2.1) leads to the following dimensionless equilibrium

equation


d2X2 + _2


(2.32)


= T


when membrane tension is dominated by elastic strain, with 11 defined to be


7 Eh
171 ~ [e00


(2.33)


d2x21 + 2\ 2d
dx2 dx 1


)- 1 P


when membrane tension is dominated by pretension, with 12 defined to be


H2 = (- Sh (2.35)
2 p2QV2 C

The use of the cube root in the definition off 1 in Eq. (2.33) is suggested by the exact solution

of Eq. (2.1) given by Seide (1977).

The boundary conditions for the membrane equilibrium equation may also be written

in nondimensional form as


(2.34)








x2 =0 at x1 = 0,1 (2.36)

where the hat has been dropped for convenience from the dimensionless variables in Eqs.

(2.32), (2.34), and (2.36).

The physical significance of the aeroelastic parameters 171 and 172 is that the

nondimensional deformation of an initially flat, elastic membrane is inversely proportional

to l1 in the absence of pretension. Alternatively, the deformation of a membrane is inversely

proportional to H2 in the presence of large initial pretension. Consequently, the steady state,

inviscid, aeroelastic response of an initially flat membrane wing at a specified angle of attack

is controlled exclusively by 17 in the limit of vanishing pretension and exclusively by 172

in the limit of vanishing material stiffness.

Substituting the same dimensionless variables into Eqs. (2.9) through (2.11) leads

to the following nondimensional form of the incompressible Navier-Stokes equations


St (IV ( Viv ) J ) ax (2.37)


Sta ( VjV2 ) aL) (2.38)


a( v ) = 0 (2.39)


with boundary conditions at the membrane surface

vi = St xi (2.40)

and far upstream from the membrane

v, = cosa(l + fsint ) (2.41)

V2 = sina(1 + #sint ) (2.42)

where again the hat has been dropped for convenience from the dimensionless variables. The
dimensionless parameters appearing in Eqs. (2.37) and (2.38), the Reynolds number and the
Strouhal number are defined as








VgQC
Rn = wC (2.43)


St = (2.44)

If the membrane is not initially taut the geometry of the wing may be characterized

by an additional dimensionless quantity, the excess length e, defined as follows.

Lo c (2.45)
c
C

In order to relate the membrane excess length to a more familiar airfoil property,

Greenhalgh et al. (1984) offered the following approximate formula, based on a parabolic

arc airfoil, which relates E to the wing camber by

X2 61 (2.46)


where the left hand side of Eq. (2.46) is taken to be the value of the x2 coordinate at the

midchord point. The set of dimensionless parameters given above 17, 112, St, Rn, E, and

a completely characterize the physical problem considered in the present work.

In addition to the basic dimensionless parameter set other physically significant

dimensionless parameters may be derived. In particular, a frequency ratio parameter, Q, may

be introduced by drawing an analogy between a one degree of freedom spring/mass system

and the transverse motion of the membrane in response to a harmonically varying freestream

velocity. If the ratio of the system forcing frequency to the system natural frequency is

defined as

S= (2.47)

the following parameters may be defined and substituted for the Strouhal number in the basic

parameter set.


1 = 3 (2.48)
111


2 = 2 (2.49)
2







Again, as in the case of steady flow, Q2 is the appropriate dimensionless frequency when the

membrane is pretensioned, whereas Q2 is the appropriate dimensionless frequency when the

membrane develops tension elastically.

Finally, the aerodynamic lift, drag and moment about the quarter-chord as well as the

membrane tension and fluid dynamic pressure are nondimensionalized in the customary way

as stated below.
L
CL = j2 (2.50)

D
CD = 1 (2.51)


Mc/4
CM = 1, M1 (2.52)


CT = v (2.53)

P -PO
cp 2 (2.54)

In Eqs. (2.50) through (2.52), the lift, drag, and quarter-chord moment are obtained

by integration of the aerodynamic pressure and shear stress over the membrane chord and

resolving the net force vector into components parallel and perpendicular to the freestream

velocity vector.













CHAPTER 3
NUMERICAL METHOD
3.1 Membrane Equilibrium

A discrete form of the elastic membrane boundary value problem can be obtained

at a finite number of points on the fixed interval [0,c] by replacing the derivatives in Eq. (2.1)

and the integral in Eq. (2.7) with appropriate finite difference and finite sum approximations.

Applying central difference approximations to Eq. (2.1) leads to a three point difference

kernel centered around point P with neighboring points E and W as shown in Fig. 2.1. At

each point P an equation of the general form

Apx2, = AEX2E + AwX2, + Sx2 (3.1)

is then obtained where the A's are coefficients associated with the finite difference

approximation, and Sx2 is a source term containing all terms in Eq. (2.1) that cannot be

expressed as a linear combination of the x2 coordinates. Consequently, all of the nonlinearity

in the boundary value problem is contained in the source term. The resulting set of finite

difference equations may then be solved using a line iteration method with underrelaxation.

As a measure of the degree to which the discrete form of Eq. (2.1) has been satisfied the

residual of the set of membrane equilibrium equations is defined as follows.

R, = I- AX2, + AEX2E + A 2w + Sx2 c (3.2)
all P

3.2 Fluid Dynamic Conservation Laws

A pressure based numerical procedure originally proposed by Patankar and Spalding

(1972) for Cartesian coordinates was chosen for computing the laminar incompressible flow

surrounding a two-dimensional wing of vanishing thickness in an unbounded domain. The

details of the basic pressure correction algorithm are given in Patankar (1980) with the








extension of the procedure to general curvilinear coordinates given in Shyy (1994). The

present implementation of the algorithm follows the work of Braaten and Shyy (1986) with

the statement of the conservation laws in terms of a time dependent grid taken from

Thompson et al. (1985).









P



S rV V22
X1



xl

Figure 3.1. Staggered grid arrangement in the physical domain




A staggered grid arrangement is adopted for the discretization of Eqs. (2.9) through

(2.11) as shown in Fig. 3.1. Specifically, a nonorthogonal body-fitted grid system is

generated numerically at the positions marked by triangles. Both the Cartesian and

contravariant velocity components, vi and Vi, as well as the Cartesian components of the grid

velocity vector, xi, are located at the midpoint of the control volume faces as shown in the

figure. The discrete value of pressure is located at the arithmetic center of the four adjacent

grid points. The finite volume form of the conservation laws may be obtained by integrating

Eqs. (2.9) through (2.11) using a first order, fully implicit time integration scheme over

appropriately staggered control volumes leading to a five point difference kernel centered

around the point P with neighbor points N, S, E, and W. By arbitrarily taking the grid spacing








in the transformed domain to be unity and the time step to be A t, the momentum and mass

conservation laws take on the following discrete control volume form



( JQv1 JOQgv )
At

'QVlV1 (ql1 2V1,2) + 2,2P)e (QV1vi -(ql1,1 q2v1,2) + 2,2)w +

V2 (- q21,1 + q3v1,2) x2,1P)n (QV2V (-21, + q31,2) x2,P)s

= 0 (3.3)



( Jov2 JOv0 )
j7 +
At

QVlv2 (qlv2,1 q2V2,2) + Xl,2P)e (QVlV2 -(ql2,1 q2V2,2) + Xl,2P) +

QV2v2 2-(- qv2,1 + q32,2) X1,1P)n (V22 (- q2V2,1 + q32,2) X1,P)s

= 0 (3.4)



(JQAt + (GVI)e (eV1), + (QV2)n (V2)s = 0 (3.5)



where the superscript refers to the previous time level. The Cartesian components of the grid

velocity vector appearing in Eqs. (2.15), (2.16) and (2.21) are approximated by the first order

backward time difference given below.

xi x9
x. =
i At (3.6)


In order to clarify some of the details of the numerical implementation consider the

control volume used to derive the discrete form of the vj momentum equation as shown in

Fig. 3.2.













: 4
|1
!
o1
W,,. .0. w l--DO- e "A v-
!I
\''wt ''!
!I


-D VIs


Figure 3.2. Control volume used for the vj momentum equation




If the convective momentum flux is approximated using a first order upwind scheme

and the diffusive flux is approximated by a central difference expression the discrete form

of the vj momentum equation may be written as


ApVP = AEVIE + Awvl, + ANVlY + ASVi + SvI


(3.7)


where, after explicitly imposing mass conservation, the coefficients in Eq. (3.7) are given

by


AE = -QVIe, 0] + iqie



A = [Q ilw o] + qilw



AN= [- eV2il, ]+' q3ln



As = [ V2s 0] + -q31s



Ap = AE+A +AN + A + JQ


(3.8)



(3.9)



(3.10)



(3.11)



(3.12)


with the source term given by


-y VIN









Sv, = q21,21e + q2V1,2w- q2v n + q2V1,


+ x2,2Plw x2,2PIe + X2,1Pln X2,IPls + dJ (3.13)

where the various terms may be evaluated at the control volume faces using appropriate

interpolations. In the present work the second order upwind scheme of Shyy, Thakur and

Wright (1992) is adopted. As a result of the higher order interpolation scheme used for the

convection terms, additional source terms will appear in Eq. (3.13). The explicit form of

these additional terms may be found in that reference.




2N2
------





V I P
S
-^-4--t


V2s


Figure 3.3. Control volume used for the v2 momentum equation

Similarly, the discrete form of the v2 momentum equation may be derived by

considering the control volume shown in Figure 3.3. The discrete form of the equation may

be written as

ApV2p = AEV2E + Awv2w + ANV2, + ASV2, + Sv, (3.14)

with the coefficients given by








AE = [-QVl, o] + q Iel


A = Q[evw o] +-+qllw


AN = [- V2n ] + 0 q In


As = [ V2s 0] + -Jq3s

AE+A+A+A
Ap = AE + Aw + AN + As + A"


(3.15)


(3.16)


(3.17)


(3.18)


(3.19)


and the source term given by


S, = --q2V2,21e + j2 ,21 J- q2,1ln + q2,1s

J+ ,2Pe ,2 + QV
+ Xl,2pIe Xi,2PIw X1,p"n + XiPIs + A---


(3.20)


Similarly, the discrete form of the pressure correction equation may be derived by

considering the control volume shown in Fig. 3.4.


+ 4


P w
*


P E


4 4 C-


I ,1*
Ps

Figure 3.4. Control volume used for the p' equation









The discrete form of the pressure correction equation may be written as

Ap p'p = AE P'E + AW P'w + AN P'N + As P's + Sp' (3.21)

with the coefficients given by


AE = +x 2] (3.22)


Ex2 x2 "
AW = Q 2 + x22 (3.23)



A = 2p A1 (3.24)

A [x21,1 x2,11

A [ A lJ A (3.25)


Ap = AE + Aw + AN + As (3.26)

and the source term given by

S, = Q(V. V, + V) + J (3.27)
At(3.27)

The sequence of events in the solution procedure is as follows. The momentum

equations are first solved using an ADI method with an initial guess for the Cartesian

velocities and the pressure field. When solving these equations, the contravariant velocities

are first computed directly from the guessed Cartesian components. The new velocity field

will generally not satisfy the mass conservation law since the pressure field used in the

solution process was only approximately correct. Consequently the pressure field must be

corrected to more closely satisfy the continuity equation.

In the basic pressure correction procedure of Patankar (1980), the correct pressure

field p is obtained from








p = P* + p (3.28)

where p* is an approximate pressure field and p' is called the pressure correction.

Corresponding contravariant velocity corrections may be introduced in a similar way

following the procedure described in Shyy (1994) and Braaten and Shyy (1986).

V1 = V + V\ (3.29)


V2 = V + V2 (3.30)

To derive the pressure correction equation we first subtract the momentum equations with

the approximate pressure and velocity fields from the same momentum equations using the

correct flow field variables. The resulting momentum correction equations are then assumed

to be adequately approximated by the following truncated form

V1p Alp x2,2 P,1 + X2,1 P,2 (3.31)


v2p A2P 2x x1,2 P, Xl,1 P,2 (3.32)

Using these approximate correction formulas, the Cartesian velocity components may be

corrected by the following relationships

S (- x2,2 P'l + X2,l P2) (3.33)
v1 = Vl + A


(x,2 P,' xl P,) (3.34)
v2 = v2 + A
A2P

Subsequently, the contravariant component correction formulas may be obtained by

substituting Eqs. (3.33) and (3.34) into Eqs (2.15) and (2.16). After dropping terms that are

not representable on a five point finite-difference molecule the following simplified

correction equations for the contravariant velocity components are obtained








/(x2,2)2 1,2)2 ,
V, = V A (, + A2P Pl' (3.35)

((x +2 (x,)2)

V2 (= 21) + (1 P,2 (3.36)

These equations are then substituted into Eq. (3.5) to obtain Eq. (3.21).
The following residuals are defined to assess the degree of convergence of the
discrete form of the fluid dynamic conservation laws.



AviP + Alnb nb + Sv,
Rv = allnb (3.37)
Ov2 C
all P 0v0c


A2Pv2P + A2nb nb + S
Rv = all nb (3.38)
all P 0vic


R-as 1 Vc (3.39)
all P
In the present procedure the pressure corrections are used to update the contravariant
components using Eqs. (3.35) and (3.36). Provided the pressure correction equation is solved
to convergence at each iteration, the resulting Vi and V2 fields satisfy mass conservation
exactly. After satisfying continuity in terms of the contravariant components, it is necessary
to recover the Cartesian velocity components before returning to the momentum equations
with the updated velocity and pressure fields. It is of preeminent importance to recover
consistently the Cartesian velocity components from the updated contravariant components
since, from Eq. (3.5), mass conservation is stated explicitly in terms of the contravariant
components in generalized curvilinear coordinates. Care must be taken to ensure a consistent
transformation between Cartesian velocity components and contravariant components. If








the transformation is done naively, an inconsistency will arise in the discrete form of the

conservation laws. This issue is briefly discussed in the following paragraph.

The obvious way to compute the Cartesian components from the contravariant

components is to invert analytically the transformation given by Eqs. (2.15) and (2.16). This

inversion gives the following formulas for the Cartesian components in terms of the

contravariant velocity components

X1,1 X1,2
Vl = V, + V2- (3.40)


2 =V2 + VI 2 (3.41)

In order to estimate the contravariant velocity components, it is necessary to interpolate for

the value of v2 when computing VI in Eq. (2.15) and similarly interpolate for the value of

vi when computing V2 in Eq (2.16). A corresponding interpolation is needed when

computing vi and v2 in Eqs. (3.40) and (3.41). The need for these interpolations causes an

inconsistency in mass conservation to arise in the procedure. An efficient method for

consistently recovering the Cartesian components from the updated contravariant

components based on the so called D'yakonov iteration procedure is described in Braaten

and Shyy (1986) and Shyy (1994). In the present numerical procedure this method for

computing the Cartesian velocity components is adopted.

3.3 Consistent Implementation of the Continuity Equation

Again, because of the need for interpolation in computing the contravariant velocity

components appearing in Eqs. (2.15-2.16) a possible inconsistency arises in the

implementation of the discrete form of the continuity equation when the grid is time

dependent. An inconsistent numerical implementation of the continuity equation would lead

to the generation of artificial mass sources in the flow calculation. As suggested by

Thompson et al. (1985) the following equation relating the time rate of change of the







Jacobian to the Cartesian components of the grid velocity is used to update the Jacobian at

the implicit time level.


t+ a1(- X2,2 + 2X1,2 )+ 2X + X-2,1 ) = 0 (3.42)

This identity may be derived from Eq. (2.14) by requiring mass to be conserved for a constant

density, uniform velocity field under a coordinate transformation that is time dependent.

Integrating Eq. (3.42) using a first order, fully implicit time integration scheme over

the same control volume used for mass conservation leads to Eq. (3.43) which is the discrete

form of the identity given above. Adopting this equation as the updating formula for the

Jacobian guarantees that the basic requirement of geometric conservation (Shyy and Vu,

1991, and Thomas and Lombard, 1979) is satisfied in the discrete form of the conservation

laws when the grid is time dependent.

(J J0)
At + (- X1X2,2 + X2XI,2)e (- XX2,2 + 2l,2)w

+ (- x2x1,1 + XX2,1)n (- X2xl,1 + XIX2,1)s = 0 (3.43)


3.4 The Aeroelastic Computational Procedure

The primary objective of the present work is to determine the equilibrium

configuration and associated aerodynamic characteristics of an elastic membrane wing as a

function of time. Consequently, the present aeroelastic problem consists of finding

membrane configurations and aerodynamic surface pressures and shear stresses which

simultaneously satisfy Eq. (2.1) and Eqs. (2.9) through (2.11), respectively. An iterative

procedure is used to solve the coupled boundary value problem by computing the elastic and

aerodynamic problems cyclically until a solution is obtained at each time step that satisfies

the governing equations to a predetermined convergence criterion.

Since a fully implicit time integration scheme is used to advance the solution from

one time level to the next, all grid dependent terms appearing in Eqs. (3.3) through (3.5) must

be reevaluated each time the membrane profile and field grid points are adjusted during the







aeroelastic iteration procedure. The Cartesian components of the grid velocity must also be

reevaluated since the grid velocity enters the solution not only through the definition of the

contravariant velocity components defined by Eqs. (2.15) and (2.16), but also through the

kinematic boundary condition given in Eq. (2.21). Consequently, the metrics of the

transformation as well as the components of the grid velocity vector are recomputed directly

from the updated grid coordinates each time the membrane equilibrium equation is solved.

An exception to this strategy is made when the Jacobian of the transformation at the implicit

time level appearing in Eq. (3.5) is evaluated. This grid dependent quantity is not computed

directly from the updated grid coordinates, but is computed using the identity given by Eq.

(3.43).

During the course of the aeroelastic iteration procedure it is necessary to update the

grid in the physical domain so that it remains body-fitted as the wing deforms. The

maintenance of the body-fitted grid is achieved by updating the grid points in the vicinity

of the membrane during the course of iteration according to the following formula

x = -1) + g(f2)(X~ x i-)) (3.44)

where g(2) is a general function that decays with distance away from the membrane surface.

The superscripts appearing in Eq. (3.44) indicate the level of the aeroelastic iteration

procedure. Presently, g(2) is chosen as an exponential function defined as


g(2) = exp 1l (3.45)

where cl is a constant that depends on the grid resolution and the Reynolds number. A similar

strategy using an exponentially decaying function for updating field grid points while

maintaining a body-fitted curvilinear grid has been used by Boschitsch and Quackenbush

(1993).








3.5 A Potential Flow Model for Thin Wings

In this section a method is presented for computing the surface pressure distribution

resulting from the irrotational, incompressible flow over thin, two-dimensional wings of

arbitrary shape. The reason for including in the present work a numerical procedure for the

solution of lifting potential flows is to highlight the difference in aerodynamic quantities

predicted by a viscous flow based membrane wing model when compared to previous work

on membrane wing aerodynamics based on a potential description of the flow field.

The fundamental assumption concerning the flow field around the wing is that it is

irrotational, and consequently the velocity field may be derived from a scalar potential.

Implicit in this assumption is that the flow is inviscid and free of vorticity far upstream. The

use of point vortex singularities to model the lifting potential flow around thin wings has its

historical origin and justification in work by James (1972). A description of the method in

its modem form may be found in Katz and Plotkin (1991).



%+ n

p-



point vortex panel



Figure 3.5 Potential flow wing model

Figure 3.5 shows a membrane airfoil that has been discretized into a number of linear

segments or panels. A typical segment is composed of a point vortex and a control point

where the flow tangency condition is enforced. If the point vortices and control points are

located at the local quarter and three-quarter chord points of the panel, respectively, the

Kutta condition is implicitly satisfied and no additional boundary condition is needed at the

trailing edge (Katz and Plotkin, 1991, James 1972). By modeling the wing as an assemblage







of vortex singularity segments of unknown strength and enforcing the zero normal relative

velocity condition at each control point, designated as point i, the following set of linear

algebraic equations may be formed and solved for the vortex strength F at point j.

Aij ij = (- v. n )i (3.46)

In Eq. (3.46) Ai is a matrix of vortex influence coefficients, defined as the normal velocity

induced at control point i by a unit strength vortex at point j, v is the freestream velocity

vector, and n is the unit normal vector at the control point.

Once the vortex strengths have been determined from Eq. (3.46) the perturbation

velocities on the upper and lower surfaces, v and v-, respectively, are given by


v 21 [- n" (3.47)



v = [ n (3.48)


where I is the length of the segment or panel and nj and n2 are the Cartesian components of

the unit normal vector.

With the velocity field determined on the upper and lower surface of the wing the

pressure on the wing surfaces, p+ and p, may be computed directly from the Bernoulli

equation as follows


p* = ( v2 v* ) (3.49)

The total fluid velocity on the upper and lower wing surfaces is the sum of the perturbation

velocity and the free stream velocity and is given by

v = v" + v* (3.50)

It should be pointed out that the leading edge suction force arising from the pressure

singularity at the leading edge of the infinitely thin membrane is not included in the above

potential flow model.















CHAPTER 4
ELEMENTARY TEST CASES

4.1 Rotated Channel Flow

Before applying the present method to an aeroelastic problem, the performance,

accuracy and convergence properties of the numerical formulations presented in Chapter 3

are investigated by applying the membrane equilibrium and fluid dynamic formulations

independently to certain relevant test problems. The first test case to be presented is that of

steady, developing parallel flow in a channel that is rotated 450 to the Cartesian axes. Once

the flow becomes fully developed the convection terms vanish in the Navier-Stokes

equations and an exact solution to the conservation laws may be obtained (Schlichting 1979).


le+02

le+01

le+00

le-01

le-02

le-03

le-04


100 200 300
iteration


400 500


Figure 4.1. Pressure contours (a), and convergence path (b), for a 450 rotated channel at
Rn=40. The results shown are for a 41 x 41 uniform grid. D'yakonov iteration was used
in the computation.

Figure 4.1 (a) shows the constant pressure contours computed using the numerical

algorithm described in the previous chapter with the D'yakonov iterative procedure used to








recover the Cartesian velocity components from the contravariant components. A uniform

velocity field is specified at the channel inlet and a 0th order extrapolation of the Cartesian

velocity components is used at the channel exit to enforce a gradient free boundary condition

in the streamwise direction. As may be seen from the figure the streamwise pressure gradient

becomes a constant after some initial adjustment of the velocity and pressure fields near the

inlet. A comparison of the computed streamwise pressure gradient with the analytic solution

for fully developed parallel flow is shown in Fig. 4.3. It may be seen from the figure that once

the flow becomes fully developed the computed pressure gradient is indistinguishable from

the analytic solution. The residuals of the discrete form of the conservation laws, defined

by Eqs. (3.37) through (3.39), are shown in Fig. 4.1 (b). The momentum equations converge

to a terminal level of 5 x 10-4 and the continuity equation converges to a terminal level of

5 x 10-5. These levels of residuals are consistent with the single precision floating point

accuracy of the arithmetic used in the calculation.





le+02



le+00

le-01 I R
R
l e-02 r R

le-03

Sle-04 i

0 100 200 300 400 500
0 (a) (b) iteration



Figure 4.2. Pressure contours (a), and convergence path (b), for a 450 rotated channel at
Rn = 40. The results shown are for a 41 x 41 uniform grid. D'yakonov iteration was
not used in the computation.








Figure 4.2 shows the results for the same problem but now the D'yakonov procedure

is not used and the Cartesian velocity components are computed directly from Eqs. (3.40)

and (3.41). Interestingly, the pressure contours look qualitatively correct but the momentum

equation residuals have converged to a terminal level which is order one. The reason for the

dissatisfaction of the momentum equations is evident in Fig. 4.3 where it may be seen in that

the pressure gradient is computed incorrectly when the D'yakonov iterative procedure is not

used. The convergence path shown in Fig. 4.2 (b) clearly illustrates the inconsistency that

arises in the transformation due to the noncollocated velocity components. Fortunately, the

convergence path shown in Fig. 4.1 (b) shows that the inconsistency can be resolved using

the D'yakonov procedure during the course of the computation.


5.0


4.0


3.0
0.
U
2.0


1.0


0.0
0 10 20
1,


30 40 50


Figure 4.3. Centerline pressure coefficient for the 450 rotated channel at Rn = 40.


4.2 Uniform Flow Using a Moving Grid

The second elementary test case to be considered is the computation of a uniform

flow on a grid that moves arbitrarily in space as a function of time. This is an important test








case since any formulation that is to be applied to unsteady flow problems using time

dependent, body-fitted coordinates must admit a uniform flow field as a solution identically

when the grid coordinates are an arbitrary function of time.

Figure 4.4 shows the grid and pressure contours after one time step for a uniform flow

field inclined at 450 to the Cartesian coordinate axes. The initial grid at time t is shown in

Fig. 4.4 (a) and the grid at time t + At is shown in Fig. 4.4 (b). The boundary conditions

imposed on the solution at all computational boundaries are Dirichlet conditions on the

Cartesian velocity components. The pressure contours computed using Eq. (3.43) to update

the Jacobian of the transformation at the t + At time level are shown in Fig. 4.4 (c). The

computed pressure gradient in the flow is approximately 10-7 everywhere in the domain and

the contours shown in the figure result from roundoff error at the limit of single precision

floating point arithmetic. Contrastingly, the pressure contours computed using Eq. (2.20) to

calculate the Jacobian at the t + At time level are shown in Fig. 4.4 (d). The pressure gradient

in this case is order one, which clearly demonstrates the inconsistency that arises in the

computation by simply taking a snapshot of the grid at the implicit time level and computing

the Jacobian using Eq. (2.20). Adopting Eq. (3.43) as the updating formula for the Jacobian

guarantees the basic conservation laws are satisfied for unsteady flows using time dependent,

body-fitted coordinates.

4.3 Elastic Membrane Under a Uniform Pressure Load

The last elementary test problem to be investigated is the deformation of an elastic

membrane under a prescribed uniform pressure load. Figure 4.5 (a) shows the computed

equilibrium shape of an initially flat membrane under uniform pressure loading. For this test

case the pretension is chosen to be zero, consequently the elastic response of the membrane

is completely described by the nondimensional elastic parameter 1lj. For the value of 7II

shown in the figure the numerical algorithm converged to machine accuracy in less than 100

iterations.



























(a)



















pressure gradient 0(10-7)


(d) pressure gradient 0(1)


Figure 4.4. Comparison of Jacobian updating strategies for the computation of a 450
uniform flow on a time dependent grid. Initial and t=t+A t grids are shown in (a) and
(b), respectively. Computed pressure contours using Eq. (4.43) are shown in (c), and
using Eq. (2.20) are shown in (d).








As the inextensible limit is approached and H1 tends to infinity the geometric

nonlinearity intrinsic in the problem becomes more pronounced and the elastic boundary

value problem becomes algorithmically stiff. Consequently, the use of underrelaxation

becomes essential to ensure convergence of the numerical scheme. This is significant since

most of the practical applications for membrane wings, as well as the available experimental

data, are for flexible but nearly inextensible membranes. In the aeroelastic membrane wing

cases presented in Chapters 5 and 6, relaxation factors as small as 10-4 were necessary to

ensure algorithmic stability.




0.1 0.1


0.08 n 0.08
E .0 -- analytic solution
Sc computed solution
0.06 0.06


0.04 0.04


0.02 0.02

0 0
(a) 0 0.2 0.4 0.6 0.8 1 (b) 0 0.1 0.2 0.3 0.4
x/c n


Figure 4.5. Computed equilibrium configuration, (a), and midpoint coordinate, (b), for
an initially flat, uniformly loaded elastic membrane. Computed results used 100 grid
points along the membrane chord.

Figure 4.5 (b) shows a comparison of the computed midpoint coordinate of the

membrane with an analytic solution for an elastic membrane given by Seide (1977). Here,

uniform pressure has been substituted for the stagnation pressure in the definition of I. The

computed solutions and the analytic results are essentially indistinguishable.















CHAPTER 5
MEMBRANE WINGS IN STEADY FLOW

5.1 Rigid Wing in a Viscous Fluid

In addition to the elementary test cases discussed in the previous chapter several other

aspects of the computational procedure need to investigated before the algorithm is applied

to a membrane wing problem. Specifically, the effect of locating the outer computational

boundary and the effect of grid refinement on the steady state solution is investigated in this

section.

5.1.1 Effect of Outer Boundary Location

Since the physical problem under consideration is the flow over a wing in an

unbounded domain the placement and boundary conditions imposed on the outer boundary

of the computational domain require investigation. Consequently, a sequence of

computations was performed to assess the sensitivity of the computed results to the location

and the boundary conditions imposed at the outer flow boundary.

Figure 5.1 (a) shows a typical grid used in the membrane wing computations. The

wing is situated in the center of a square computational domain typically 10c X 10c in size.

The freestream velocity is imposed on all outer flow boundaries except the downstream

boundary where a zero gradient condition is imposed on the velocity components. An

enlarged view of the grid in the vicinity of the wing is shown in Fig 5.1 (b). The refinement

of the grid near the membrane was based primarily on experience and the degree to which

flow separation was expected in each situation. Attention to grid refinement is particularly

important near the trailing edge and in the wake where large regions of recirculating flow

often occurred.








Dirichlet







ill IlI IiI Ill I I Ii

















Dirichlet









(b ) ..... .. .

Figure 5.1. Typical grid used for membrane wing computations. The complete
computational domain and boundary conditions are shown in (a), and an enlarged view
of the grid near the membrane is shown in (b).

A series of computations was performed using different conditions at the outer

domain boundary prior to selecting the boundary condition set shown in Fig. 5.1 (a). These

included computations using Dirichlet conditions at all outer domain boundaries as well as

computations where a zero gradient condition was imposed on all boundaries except the







upstream boundary where freestream conditions were imposed. In the latter case with exit

conditions specified at the north, south, and east domain boundaries, the algorithm was

unable to admit a uniform flow field as a solution to the fluid dynamic conservation laws.

Consequently, the exit boundary condition was retained only at the downstream domain

boundary as shown in Fig. 5.1 (a).

Table 5.1. Effect of outer boundary location on computed aerodynamic properties of a
rigid, 2% camber, zero thickness, circular arc airfoil at a = 50, Rn = 4x103.
grid domain size CL CD CM
185x 77 5cx 5c .637 .0761 .0442
207x 91 10cX 10c .626 .0752 .0425
221x 101 15cx 15c .623 .0750 .0422
247x 121 30cx 30c .621 .0749 .0420


The effect of moving the outer boundary away from a rigid, zero thickness flat plate

at 40 angle of attack and a Reynolds number of 4x 103 may be seen in Table 5.1. The boundary

was moved outward by adding additional grid points to the basic 5c X 5c grid so that the outer

boundary location is isolated from other grid dependency issues. For all grids 100 grid points

were used along the wing chord. As may be seen from the table there is very little change

- less than 1% in the computed lift, drag and aerodynamic moment beyond a domain size

of 10c X 10c. Consequently, a computational domain size of 10c X 10c with exit conditions

imposed at the downstream boundary is adopted for all further computations.


5.1.2 Effect of Grid Refinement

In order to asses the effect of grid resolution on the computed flow field near a

membrane wing two test cases were chosen for investigation. The first test case was a rigid,

zero thickness flat plate at 40 angle of attack and a Reynolds number of 4x103. Figure 5.2

shows the computed streamlines for the flat plate using 100, 200, and 300 grid points along

the wing chord. As may be seen in the figure a small separation bubble emerges at the leading

edge and grows as the grid is refined. The computed lift, drag, and moment taken about the

wing quarter chord are presented in Table 5.2 for the grid refinement sequence of Fig 5.2.







Table 5.2. Effect of grid refinement on computed aerodynamic properties of a rigid flat
plate at a = 40, Rn = 4x103
grid points on c CL CD CM
161 x 81 100 .401 .0673 .0012
281x 101 200 .404 .0604 .0009
401x 121 300 .413 .0572 .0022


The table shows that drag, as expected, is the quantity most strongly affected by grid

resolution. The computed lift varies by less than 3% over the range of grid resolution

investigated. The computed aerodynamic moment is effectively zero for all grids. This is

consistent with the analytic result from thin wing theory (Katz and Plotkin 1991).


Figure 5.2. Effect of grid refinement on streamfunction contours for a rigid flat plate at
a=4", Rn=4x103.








Figure 5.3 shows the computed aerodynamic surface pressures at Rn = 104 for a rigid

flat plate at 30 angle of attack. The computed surface pressures used a grid with 100 points

along the wing chord. Also shown in the figure for comparison is the potential flow solution

obtained from conformal mapping for the flat plate. The analytical solutions for the potential

flow is taken from Katz and Plotkin (1991). The computed surface pressures for the flat plate

are in good agreement with the potential flow solution since the flow is largely attached and

the boundary layer is fairly thin at this Reynolds number.


0 0.2 0.4 0.6 0.8 1
x,/c


Figure 5.3. Comparison of computed surface pressures at Rn=104 and a=30 with
potential flow solution for a rigid, zero thickness flat plate.


Figure 5.4 shows the computed lift curves at Reynolds numbers of 104 and 4 x 103

for a rigid flat plate using 100 grid points along the wing chord. Also shown in the figure

for comparison is the potential flow lift curve taken from Katz and Plotkin (1991). Again,

as in Fig. 5.3, the viscous flow solution shows close agreement with the potential flow

solution at these Reynolds numbers and angles of attack.













Rn=10
--- Rn=--4x 10
0.8 potential flow


0.6


0.4


0.2


0
0 2 4 6 8 10
a



Figure 5.4. Comparison of computed lift coefficient with potential flow solution for a
rigid, zero thickness flat plate.




The second test case for investigating the effect of grid refinement on the computed

flow field is a rigid, 2% camber, circular arc airfoil at 50 angle of attack and a Reynolds

number of 4x 103. Figure 5.5 shows the computed streamlines for the circular arc airfoil using

100, 200, 300 and 400 grid points along the wing chord. As was seen with the flat plate grid

refinement sequence a leading edge separation bubble appears as the grid is refined.

However, the most significant flow feature to emerge as the grid is refined is the recirculation

region near the trailing edge and the attendant departure of the streamlines from the upper

surface of the airfoil. The wake behind the airfoil may also be seen to thicken significantly

with the appearance of the trailing edge separation region. Based on the grid refinement

sequence shown in Figure 5.5, it may be concluded that for circular arc configurations -

associated with either flexible or rigid wings potential flow will not be capable of

accurately describing the flow field.




































400 points on c








Figure 5.5. Effect of grid refinement on streamfunction contours for a rigid, 2% circular
arc airfoil at a= 5, Rn=4xl03.


The computed lift, drag, and aerodynamic moment for the circular arc airfoil are

presented in Table 5.3 for the grid refinement sequence shown in Fig 5.5. In contrast to the

flat plate results the computed values for the lift, drag, and aerodynamic moment are quite

sensitive to grid refinement. Even at the highest grid resolution there is a 5% 10% change

in computed values when compared to the previous grid. However, it is expected that a grid

independent solution will be difficult to obtain for this problem since a stringent requirement

is placed on grid resolution at this Reynolds number.


I I




44


Table 5.3. Effect of grid refinement on computed aeroelastic properties for a rigid, 2%
circular arc airfoil at a = 50, Rn = 4x103
grid points on c CL CD CM
207x 91 100 .626 .0752 .0425
311x 93 200 .595 .0654 .0394
411x 95 300 .567 .0603 .0378
531 x105 400 .531 .0558 .0352


Based on the results for a circular arc airfoil it may be expected that computations

with a flexible membrane will exhibit a similar sensitivity to grid refinement since the

equilibrium configuration for a pressure loaded, end restrained membrane is quite similar

to a circular arc and will likely exhibit the same type of trailing edge separation.


0.2 0.4 0.6 0.8 1
xl/C


Figure 5.6. Comparison of computed surface pressures at Rn=104 and a=3 with
potential flow solution for a rigid, zero thickness, 2% camber, circular arc airfoil.

Figure 5.6 shows the computed aerodynamic surface pressures at Rn = 104 for a rigid,

2% camber circular airfoil at 30 angle of attack. The computed surface pressures used a grid





45


with 100 points along the wing chord. Again, the potential flow solution obtained from

conformal mapping for a circular arc airfoil taken from Katz and Plotkin (1991) is shown

in the figure for comparison. The computed surface pressures for the arc airfoil are in

reasonable agreement with the potential flow solutions. Figure 5.7 shows the computed lift

curves using 100 grid points along the wing chord at Reynolds numbers of 104 and 4x103

for the rigid, circular arc airfoil. Inspection of Figs. 5.4 and 5.7 indicate there is a general

loss of circulation, compared to potential flow theory, around the airfoil due to viscous

effects near the trailing edge. The influence of viscosity is much more pronounced for the

circular arc than for the flat plate at this angle of attack and Reynolds number.


1.2

1

0.8

S0.6

0.4

0.2


2 4 6 8 10


Figure 5.7. Comparison of computed lift coefficient with potential flow solution for a
rigid, zero thickness, 2% camber circular arc airfoil.


5.2 Flexible Membrane Wing in a Viscous Fluid

Membrane wings may be broadly classified by considering the following three

limiting cases of the parameters 171, 172, and E. The first limiting case may be referred to

as the elastic wing case with








11 finite



e =0
E=O

In this case the wing is initially flat and taut but with no pretension. Consequently, the

membrane behavior is determined exclusively by 11.

The second limiting case for membrane wings may be referred to as the constant

tension wing case with

H1 =0

172 finite

E =0

As with the first case the wing is initially flat and taut but with now with sufficient pretension

so that the material stiffness does not participate in determining the equilibrium membrane

configuration.

The third limiting case may be referred to as the inextensible wing case with

7; infinite

172 =0

S- finite

In this limiting case the wing is not initially flat and taut but has an excess length of

material defined by the parameter E. Of course, the Reynolds number and the angle of attack

are additional parameters that characterize a steady flow about an airfoil, but they are to be

distinguished from the three parameters listed above that characterize flexible membrane

wing behavior.

5.2.1 Elastic Membrane Case

The computed streamlines and constant pressure contours for an elastic membrane

wing typical of the first limiting case with HI1=10, 12-=0, E=0, are shown in Fig. 5.8. The

solution shown in the figure was computed using a grid with 100 points along the wing chord.

The recirculation region near the trailing edge is reminiscent of the flow over a rigid circular








arc airfoil as discussed above. Consequently, it may be anticipated that further refinement

of the grid would lead to the emergence of additional flow features, such as a leading edge

separation bubble, as was demonstrated with the rigid arc airfoil.


Figure 5.8. Streamlines and constant pressure contours for an initially flat elastic
membrane wing, Rn=4x103, a=60 171 =10, 12=0, E=0.


Figure 5.9 illustrates the convergence path of the four coupled governing equations

for the the elastic wing case of Fig. 5.8. It may be seen that the residuals, as defined by Eq.

(3.2) and Eqs. (3.37) through (3.39) are reduced to single precision machine accuracy after

a few thousand iterations. Also, the terminal level of the residuals shown in the figure is

consistent with the single precision floating point accuracy of the arithmetic used in the

calculation. It may be noted that the stability of the aeroelastic algorithm is largely dictated

by the choice of relaxation parameter used in the solution of the membrane equilibrium

equation. For the solution shown in the Fig. 5.8 and Fig. 5.9, a relaxation value of 10-2 was

used during the solution of Eq. (3.1).


I












le+00
Sv, residual
Sv residual
le-01 --- mass residual
--- x, residual

le-02


l, le-03


le-04



0 1000 2000 3000
iteration


Figure 5.9. Convergence path of the momentum, continuity, and equilibrium equations
for an initially flat elastic membrane wing, Rn=4x103, a=60, 1I =10, 12 -0, E=0.





Figure 5.10 compares the computed membrane surface pressures for the case of Fig.

5.8 with the surface pressures predicted by a potential flow calculation for the same

membrane configuration. Again, as with Fig. 5.7, the discrepancy between the potential

solution and viscous flow solutions appears to be attributable to a general loss of circulation

about the wing due to viscous effects primarily near the trailing edge. However, for this set

of aeroelastic parameters, potential flow does reasonably well in predicting the membrane

surface pressures. For other choices of the dimensionless parameters, potential flow is not

as successful in accurately describing the fluid dynamics of the problem. In the following

section describing the aeroelastic behavior of an inextensible membrane wing of moderate

excess length, the limitations of a potential based description of the flow field will be

demonstrated.












-3
---Rn=4.xlO
potential flow
-2



U -1



0



0 0.2 0.4 0.6 0.8 1
X,/C


Figure 5.10. Equilibrium surface pressures for an initially flat elastic membrane wing,
Rn=4x103, a=60, H7=10, 72=0, E-=. Also shown is the potential flow surface pressure
computed for the same membrane configuration.



Figure 5.11 (a) shows the equilibrium lift coefficient for the initially flat elastic wing

of Fig. 5.8. The figure shows the dependency on Reynolds number to be more pronounced

for the elastic wing when compared to the rigid circular arc airfoil of Figure 5.7. This may

be explained by observing that the maximum wing camber, shown in Fig. 5.11 (b), follows

the same trend as the lift coefficient. It may be concluded that this enhanced Reynolds

number dependency is largely a result of the coupling between the membrane deformation

and the flow field. Also shown in the figure is the equilibrium lift coefficient predicted by

a potential based membrane wing model for the same set of aeroelastic parameters. It may

be seen that the potential based solution considerably overpredicts the lift on the wing even

at small camber ratios and angles of attack when compared to the viscous flow based solution

at these Reynolds numbers. As stated previously this discrepancy can be attributed to a

general loss of circulation about the wing due to viscous effects near the trailing edge.











i0.04
Rn=4.xl0 -- Rn= 1i
po.l -Ra=4xlO
0.8 I -- potential
0.035

0.6


0.4

0.025
0.2


0 0.02
0 2 4 6 8 0 2 4 6 8
(a) a (b) a


Figure 5.11. Equilibrium lift coefficient, (a), and x2 coordinate at midchord, (b), for an
initially flat elastic membrane wing, 7i=10,712=0, =-0.



5.2.2 Constant Tension Membrane Case

Figure 5.12 shows the equilibrium lift coefficient and maximum wing camber for the

second limiting case of a constant tension membrane with 17 =0, 172=2, and =-0. Again, 100

grid points were used along the membrane chord. This figure may be contrasted with Fig.

5.11. The contrast between constant tension and elastic wing characteristics may be

reconciled by recalling that the response of an initially flat, constant tension membrane to

a specified uniform pressure load is qualitatively different from the response of a membrane

with finite material stiffness. For small deflections the constant tension membrane will

exhibit linear elastic behavior whereas an elastic membrane with finite material stiffness will

exhibit a 1/3 power law response to the applied pressure load (Seide 1977) as shown in Fig.

4.5. An inspection of Figs. 5.12 and 5.11 illustrates the different behavior of the two

membrane wings near zero degrees incidence. The initially flat, unpretensioned elastic

membrane wing of Fig. 5.11 has very little initial resistance to deformation, and

consequently, a substantial amount of camber exists in the equilibrium configuration near

zero degrees angle of attack. Conversely, the pretensioned membrane wing of Fig. 5.12 has









adequate stiffness in its initially flat unstrained configuration to postpone the development

of camber until higher angles of attack are reached. For the constant tension wing of Fig. 5.12

the lift curve closely mimics the the development of camber with angle attack as was the case

with the elastic wing of Fig. 5.11. Streamlines and pressure contours are not presented for

this limiting case since they are similar to the first limiting case shown in Fig. 5.12.




1 ,0.08
Rn=10 Rn=10'
Rn=4.xl03 -- Rn=4.xl0'
0.8
0.06

0.6
1 /0.04
0.4

0.02
0.2


0 2 4 6 8 0 2 4 6 8
(a) a (b) a


Figure 5.12. Equilibrium lift coefficient, (a), and x2 coordinate at mid chord, (b), for an
initially flat, constant tension membrane wing, 171=0,172=2, E=O.


5.2.3 Inextensible Membrane Case

Finally, the limiting case of a flexible but inextensible membrane wing with excess

length E=.017,171 =46 and 72=0 is illustrated in Fig. 5.13. The excess length and aeroelastic

parameters were chosen to reproduce a portion of the experimental data reported by Newman

and Low (1984). Experience has shown that for this excess length any further increase in H7

will not substantially effect the aeroelastic solution. Consequently, the results shown in Fig.

5.13 are for a membrane wing that may be considered essentially inextensible. Due to the

assumption of steady laminar flow the numerical computations were limited to Reynolds

numbers below 4x103 whereas the Reynolds number reported by Newman & Low (1984)

was 1.2x105 which is in the turbulent flow regime.




























300 points on c







400 points on c









Figure 5.13. Effect of grid refinement on streamfunction contours for an inextensible
membrane wing at a=0, 17 =46,172=0, E=.017, Rn-4xl03.


The effect of grid refinement on the computed streamlines for the inextensible

membrane wing at Rn-4xl03 and 0 angle of attack is shown in Fig. 5.13. As may be seen

in the figure a large recirculation region appears on the lower surface of the membrane as

the grid is refined. The integrated aerodynamic properties and the nominal membrane

tension are given in Table 5.4 for the grid refinement sequence of Fig. 5.13. Even for the







highest grid resolution shown the tabulated values of lift and tension are not yet grid

independent. Of course, the membrane lift and tension will follow the same trend as the grid

is refined as required by equilibrium. The highest resolution grid shown in the figure

required 24 hours of CPU time on DEC 3000 AXP workstation.


Table 5.4. Effect of grid refinement on computed aeroelastic properties for an
inextensible membrane wing at a-=00, 117=46,172--0, =.017, Rn=4xl03
grid points on c CL CD CM CT max x2
261x 101 100 .378 .0716 .152 .893 .0836
361x 121 200 .342 .0689 .144 .823 .0835
461x 141 300 .317 .0685 .139 .773 .0835
561x 141 400 .300 .0702 .136 .741 .0834


Figure 5.14 contrasts the viscous flow based aeroelastic solution using 400 points

along the wing chord with a potential flow based aeroelastic solution for the inextensible

membrane wing of Fig. 5.13. Although the equilibrium membrane profiles are similar, as

shown in Fig. 5.14 (a), the surface pressures, shown in Fig. 5.14 (b), are dramatically

different. In particular near the leading edge the pressure jump across the membrane, Ap,

predicted by the two flow models are of different sign. The negative Ap computed by the

viscous flow model implies an inflection point in the equilibrium membrane profile which,

upon close inspection, may be seen in the figure. In contrast the potential based model

predicts a flow pattern and a membrane profile that is completely symmetric about the

midchord. It may be concluded that for this set of dimensionless parameters, a potential flow

based membrane wing theory essentially fails.

Figure 5.15 (a) compares the computed, laminar flow equilibrium lift coefficient at

Reynolds numbers of 2x103 and 4x103 with the experimental data of Newman and Low

(1984) who reported a Reynolds number of 1.2x105 for their experiments. The solution

shown in the figure was computed using a grid with 100 points along the wing chord. The









disagreement between the computational and experimental results is largely due to the

assumption of laminar flow in the present computation.


0.04


0 0.2 0.4 0.6 0.8 I 0 0.2 0.4
,/c (b) /c


Figure 5.14. Comparison of potential and viscous flow based membrane configurations,
(a), and surface pressures, (b), for an inextensible membrane wing with excess length.
The dimensionless parameters are 7lj=46, 12=0, E=.017, a=0.


2



1.5







0.5



0
0 2 4 6
0 2 4 6


8 10 12 14 0 1 2 3
a ac


Figure 5.15. Comparison of lift coefficient with experimental data, (a), and separation
point, (b), for an inextensible membrane wing with excess length, 171=46, 172=0, E=.017








The decrease in lift with increasing Reynolds number shown in Fig. 5.15 (a) is

contrary to the trend observed in the previous computations on membrane wings. This trend

may be explained by observing that the separation point, xlSP (defined as the point of

vanishing shear stress), moves forward along the membrane as the Reynolds number is

increased as shown in Fig. 5.15 (b). The presence of turbulent mixing in the boundary layer

would undoubtedly postpone flow separation and lead to higher lift coefficients consistent

with the experimental data. Also shown in Fig. 5.15 (a) is the computed equilibrium lift using

a potential based membrane wing theory. The potential flow based aeroelastic solution is

seen to considerably overestimate the lift on the membrane wing when compared to the

experimental results. Conversely, the viscous, laminar flow based model considerably

underestimates the lift reported by Newman and Low (1984).














CHAPTER 6
MEMBRANE WINGS IN UNSTEADY FLOW

In this chapter the role of viscosity in the unsteady flow over a membrane wing is

investigated. In many situations viscosity may play a pronounced role in the unsteady

aerodynamics of membrane wings since the membrane will adjust to changes in the

freestream velocity as required by equilibrium. This adjustment in the membrane shape may

lead to the periodic appearance or disappearance of regions of separated flow. In addition

to investigating the role of viscosity in the aeroelastic boundary value problem, the

simulation of a membrane wing subjected to a harmonically varying freestream velocity also

serves to demonstrate the differences between the three limiting aeroelastic cases discussed

in the previous chapter. On a more practical note, the simulation of a membrane wing

subjected to a harmonically varying freestream is a rudimentary yet meaningful model

problem for investigating the behavior of a marine sail in a wind gust. For each of the three

limiting cases presented in the following sections the freestream velocity varies by 20%

about the mean value.

The time step required to resolve adequately the major features of the unsteady flow

about a membrane wing may be determined by comparing the aeroelastic solutions

computed using several different time steps. The constant tension wing was chosen for this

study since experience has shown this case to be quite responsive to a harmonically varying

freestream velocity.

Figure 6.1 shows the time series of the aeroelastic solution using dimensionless time

steps of .075, .15, and .30 with a 281 x 101 spatial grid. The freestream velocity,

aerodynamic lift coefficient, and membrane x2 coordinate at the midchord point for each of

the three time steps are shown in the figure. The dimensionless parameters that govern the








solution are given in the figure caption. It may be seen from Fig. 6.1 that even for the smallest

time step the aeroelastic solution is not yet time step independent. It is interesting to note that

the solution is essentially independent of the time step during the portion of the cycle when

the freestream is accelerating, and conversely, strongly dependent on the time step when the

freestream is decelerating. The reason for this strong timestep dependency during periods

of freestream deceleration will be discussed in the following section. However, on a more

practical note, it is known from the previous chapter that the steady state solution is also not

completely grid independent at this level of spatial resolution. Therefore, the spatial and

temporal accuracy afforded by a 281 x 101 grid and a nondimensional time step of .075 is

considered adequate for the present investigation and is adopted for all unsteady

computations. This choice is simply an accommodation to limitations in the available

computing resources.






1.4 -.-e v
CL, At=.075
1.2 x2* 10 @ c12, At=.075
CL, At=-15
,*- -- 10 @ c2, At=.15
1 p C At=.30
S--- x 10 @ c/2, At=.30
0.8

0.6

0.4

0.2

0

-0.2
0 2 4 6 8 10 12
time

Figure 6.1. Variation in the computed aerodynamic lift and midchord x2 coordinate
using three different timesteps for a constant tension membrane wing subjected to
harmonic freestream oscillation. The nondimensional parameters are E=0, Rn=4x 103,
St= 1.5, 11 =0, H2 =2, and a=4.








6.1 Constant Tension Membrane Case

The first case to be investigated is that of an initially flat, constant tension wing.

Figure 6.2 (a) shows the time series of the free stream velocity, the aerodynamic lift, drag

and moment, and the x2 coordinate at the midchord point for the membrane wing at 40 angle

of attack. The initial condition for this simulation is taken to be the steady state solution to

the aeroelastic problem with the same dimensionless parameters. Also, the time shown in

Fig. 6.2 (a) is the dimensionless time defined by Eq. (2.26). Figure 6.2 (b) and 6.2 (c) show

the membrane surface pressures and membrane profile, respectively, at several instants

during the second complete cycle of the freestream velocity.

Several observations may be made concerning the solution shown in Fig. 6.2 (a).

First, it may be seen that the peak in the membrane deflection lags the peak in the freestream

velocity by approximately 600. This phase lag, as well as the large amplitude of the motion

of the membrane the membrane camber varies from negative 1% to positive 8% is

characteristic of a system that is being driven near, but below the system natural frequency.

This observation is supported by the value of the frequency ratio 922 (Eq. 2.49), which has

a value that is approximately 1.0 for this case. Secondly, the membrane profiles shown in

Fig. 6.2 (c) demonstrate that the algorithm is capable of successfully simulating the

dynamics of membrane wings when substantial changes in the membrane profile occur.

The variation in the membrane profile may also be seen in Fig. 6.3 where the

streamfunction contours are shown at equally spaced intervals during the second complete

cycle of the freestream velocity. In this figure it may be seen that as the freestream velocity

decelerates the flow separates along the upper surface of the membrane and two regions of

recirculating flow develop near the trailing edge. The periodic appearance and collapse of

these recirculation zones, as seen in Fig. 6.3, suggests the role of viscosity is enhanced in the

unsteady flow scenario since the acceleration and deceleration of the freestream velocity

strongly influences the separation and reattachment of the flow. The periodic appearance and

collapse of these flow features, along with an attendant adjustment in the membrane








configuration, results in an aeroelastic response which may not be characterized as simple

harmonic response at the freestream forcing frequency. A leading edge separation bubble

may also be seen to appear near the end of the freestream cycle and then collapse as the

freestream velocity accelerates. The associated pressure contours are shown in Fig. 6.4.

6.2 Elastic Membrane Case

The second unsteady case to be investigated is that of an initially flat elastic wing.

In this case the membrane develops tension as a result of material strain rather than

pretension. The controlling aeroelastic parameter 17j was chosen so that the steady state

aeroelastic solution is the same as the constant tension case previously discussed. Here, as

with the previous case, the steady state solution was used to initialize the solution at time t=O.

The time series of the aerodynamic and membrane properties are shown in Fig. 6.5

(a). Again, several observations may be made concerning this time series. First, the lift and

membrane deflection show very little deviation from the steady state value when compared

to the deviations observed in the constant tension case. The difference between the two cases

may be explained by recalling that the deflection of an initially flat, elastic membrane

without pretension is proportional to the cube root of the applied pressure. Consequently, an

elastic wing will become progressively stiffer as the freestream velocity is increased. This

effect is usually referred to as geometric or stress stiffening (Leonard 1990). In this case the

effect of geometric stiffening serves to restrict the development of wing camber, and

consequently, the appearance of regions of separated flow is less pronounced than in the

constant tension case of the previous section. The membrane deflection, as well as the

aerodynamic properties, may also be seen to be essentially simple harmonic response at the

freestream forcing frequency.

An interesting feature of the elastic wing case is the apparent lack of a phase shift

between the freestream velocity and the membrane deflection. This situation is consistent

with a system that is being driven well below its natural frequency. This observation is

supported by the frequency ratio Q1 (Eq. 2.48), which for this problem is approximately 0.1.





60












0---- *Co
1.4 V,


\ x,*10 @ c/2




0.6



0.2



-0.2
(a) 0 5 10 15
time



-3
F----=- o 0.09 o* ]
I 2 -- = 14e144
-2 ---=26
0.07 288

i
\ 0.05
0 0.03 \
1 0.013


2 0.01

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8
(b) xI/c (c) x /c


Figure 6.2. Time series of system response for a constant tension membrane wing
subjected to harmonic freestream oscillation is shown in (a). Surface pressures and
membrane profiles at several times during one complete oscillation of the freestream
velocity are shown in (b) and (c), respectively. The nondimensional parameters are
E-O, Rn=4x103, St=1.5, 17 =0, 12 =2, and a=4.































= 72.










S= 144.0


L. j


0 = 216.


= 288. 0


Figure 6.3. Streamfunction contours at equally spaced time intervals during one
complete oscillation of the freestream velocity for the case of the previous figure.


111! It


T-70-011


s
























= 72. --/ -///-\\\\ \ \








c=144.








S= 216.








= 288.0

Figure 6.4. Constant pressure contours at equally spaced time intervals during one
complete oscillation of the freestream velocity for the case of the previous figure.





63











2


CL
CM
CD
)=0 ,*1 c0@/2









(a) O
0 5 10 15
time



-3 0.1

----* = 216' -- = 72
-2 0.08 *---= 144
*=-- 216'
-- =288'

-1 0.06


0 0.04


10.02


2 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8
(b) x,/C (c) x,/c


Figure 6.5. Time series of system response for an elastic membrane wing subjected to
harmonic freestream oscillation is shown in (a). Surface pressures and membrane
profiles at several times during one complete oscillation of the freestream velocity are
shown in (b) and (c), respectively. The nondimensional parameters are E=0, Rn=4x103,
St=1.5, 1i =7.9, 72 =0, and a=4.










o = 72.0 1


S=144.0 1


(q = 216.


S= 288.0


Figure 6.6. Streamfunction contours at equally spaced time intervals during one
complete oscillation of the freestream velocity for the case of the previous figure.


a~f;--~


~-~

























S= 72.








p = 144.0








p = 216.








p = 288.


Figure 6.7. Constant pressure contours at equally spaced time intervals during one
complete oscillation of the freestream velocity for the case of the previous figure.








The stream function contours for this case are shown in Fig. 6.6 and the associated

pressure contours are shown Fig. 6.7. It may be seen by comparing Fig. 6.3 and Fig. 6.6 that

although the two cases have nearly identical steady state behavior (recall the solution at t=O

is the steady state solution in both figures) the unsteady response of the two cases is

dramatically different. The differences between the two cases may be broadly attributed to

the geometric nonlinearity inherit in the elastic problem.

6.3 Inextensible Membrane Case

The final class of problems to be investigated is the case of an elastic wing with the

material stiffness tending to infinity. The controlling dimensionless parameter E, was chosen

to coincide with a portion of the experimental data reported by Newman and Low (1984).

The value of the aeroelastic parameter I17, was chosen large enough so that the membrane

may be considered essentially inextensible. Figure 6.8 (a) shows the time series of the

aerodynamic and membrane properties for the inextensible wing at 0 angle of attack. It may

be seen from the figure that the x2 coordinate of the membrane at the midchord point is

largely unaffected by variation in the freestream velocity. However the membrane profile

does change shape in response to the unsteady flow field. This shape change is shown is Fig.

6.8 (c) at two time instants during the second complete cycle of the freestream velocity. For

this class of problems the membrane profile is responding primarily to the periodic

appearance and collapse of regions of recirculating flow on both the upper and lower

surfaces of the membrane. These regions of recirculating flow may be seen in Fig. 6.9. It may

be seen that as the freestream begins to decelerate a region of separated flow appears on the

lower surface of the membrane and a pair of counterrotating recirculating zones appear on

the upper surface of the membrane near the trailing edge. As the freestream velocity begins

to accelerate these recirculating zones quickly collapse and the flow once again becomes

fully attached. As was the case with the constant tension wing the aeroelastic response shown

in Fig. 6.8 may not be characterized as simple harmonic response at the freestream forcing

frequency. The associated pressure contours are shown in Fig. 6.10.





67








1.5


--


-- x,*10 @ c/2







0.5






(a) o
0 5 10 15
time


-2 .0.1
---- =0" ---- >=0"
= 216" -O= 216'
0.08
-1

0.06
0
0.04


0.02


2 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
(b) x1/C (c) x1/c


Figure 6.8. Time series of system response for an inextensible membrane wing
subjected to harmonic freestream oscillation is shown in (a). Surface pressures and
membrane profiles at several times during one complete oscillation of the freestream
velocity are shown in (b) and (c), respectively. The nondimensional parameters are
e=.017, Rn=4x103, St=1.5, 17 =46, 172 =0, and a=0.
















0 =o0.


S= 144.


p = 216.


, = 288.


Figure 6.9. Streamfunction contours at equally spaced time intervals during one
complete oscillation of the freestream velocity for the case of the previous figure.
















q=0.0


= 144.








S= 216.








= 288.0


Figure 6.10. Constant pressure contours at equally spaced time intervals during one
complete oscillation of the freestream velocity for the case of the previous figure.













CHAPTER 7
A THREE-DIMENSIONAL MEMBRANE WING MODEL
7.1 The Elastic Membrane Problem in Three-Dimensions

In this section the principle of virtual work is used to formulate the finite element

matrices for a bilinear quadrilateral, plane stress membrane element. The total Lagrangian

formulation is used in deriving the element which may undergo large displacements and

rotations but is assumed to undergo only small strains. The formulation and nomenclature

used follows the work of Bathe (1984). Index notation with the summation convention is

used throughout and a comma implies partial differentiation with respect to the indicated

independent variable.



SYt+

'V

X2 J configuration at time 0

configuration at time t + At

X1 configuration at time t
0X3

Figure 7.1. Deformation trajectory of material volume OV

Figure 7.1 shows a material volume undergoing large deformation from a natural

unstrained state at time t=O through an intermediate state at time t to a final deformed

configuration at time t+A t. The time-like coordinate t is used here only as a convenient

parameter for describing the loading and deformation state and does not imply the existence

of a dynamic process. The statement of equilibrium for the continuum at configuration t+A t

is given by the principle of virtual work as








St+AtSi 6 ( t+Atj) odV = t+AtR (7.1)
o oV

where Si is the 2nd Piola-Kirkoff stress tensor, ei is the Green's strain tensor, 6 is an
admissible variation and t+AtR is the virtual work performed by external forces on the
material volume at time t+A t.
Consider next the incremental decomposition of the stresses, strains, displacements,
and material coordinates given by
t+AtE = tij + E~i (7.2)

t+AtSij = t'i + Sij (7.3)

t+Atui= 'ui + ui (7.4)

t+txi = txi + Ui (7.5)

where ei Sy and ui are increments in strain, stress, and displacement, respectively. Using
these decompositions in Eq. (7.1) along with the following approximations

Sij Cirs ers (7.6)

eij 6 e i (7.7)

gives the following linearized incremental equation of equilibrium

Cirs ers. ei, dV+ f tSi, 6 j dV = t+AtR t Si 6 eij dV (7.8)
J oV J oV J oV

where C,rs is the material property tensor, and ei and r/i are the linear and nonlinear
components, respectively, of the incremental strain tensor given by

e = .I( + U + t k,i Ukj + Uki kj ) (7.9)

S= ( Uki Ukj ) (7.10)
2 (U k








Having stated the general incremental equilibrium condition in Eq. (7.8) consider

now the development of the finite element matrix expressions for a plane stress,

isoparametric, bilinear quadrilateral element shown in Fig. 7.2.


local system


global system


Figure 7.2. Plane stress quadrilateral element with local and global coordinate systems

The element coordinates, displacements, and displacement increments are

interpolated from their nodal values in the element coordinate system as follows


x = [H][X XI X2 X .X4 T (7.11)
31


U1l

t U3



u2 = [ H ] 1[
u3


t U tUl .U2 t4 T



U2 3 U2 3
ulu u ..u4i


h l0 Oh2 Oh3 0 Oh4 0 01
[H] = 0 h2 0 0h20 0h30 Oh40
S00 h 00 0h2 0h0 Oh4


with the isoparametric interpolation functions hi, h2, h3, and h4 given


(7.12)


where


(7.13)


(7.14)








hi = ( 1 r ) ( 1 s ) (7.15)


h = ( 1 + r ) ( 1 s) (7.16)


3 = 1 ( + r ) ( 1 + s ) (7.17)


h4 = ( r) ( 1 + s) (7.18)

where r and s are natural coordinates defined on the bi-unit square. Substituting these
coordinates and displacements interpolations into Eq. (7.8) and invoking the principle of
virtual work for each nodal displacement in turn produces the following set of discrete
equilibrium equations in index free form

( [ 'KL ] + [ KN ] ) [ U ] = [+tR ] (i-) [ +AtF ](i-1) (7.19)

where the superscript indicates the iteration level of the solution procedure. The element
stiffness, force, and load matrices appearing in Eq. (7.19) are defined as follows

[ 'KL ] = f [ 'BL, ][ C ] [ 'BL ] dV (7.20)
J oV


[ tKN ] = J[ 'BN]T[ 'S ][ 'tBN] 0dV (7.21)
J oV


[ t+AtF ] = [ t+AtBL]T[ +S] dV (7.22)



[t+'tR] = I [H ] [t+Atf]odS (7.23)


where the explicit expressions for both the linear and nonlinear strain displacement matrices,
[ BL ] and [ BNL i, respectively, may be found in Bathe (1984).
At this point a departure is made from the general three-dimensional continuum
formulation by assuming a state of plane stress to exist in the element and employing a plane








stress consitutive law. These assumptions lead to the following choice for the material

property matrix


(7.24)


1 0
[c] E ) v 1 0
V 2 0 0 0- v
2


,where v is Poisson's ratio, and for the stress state matrices


[tS] =


S11
'S12
0
0
0
0


tS12
tS22
0
0
0
0


0
0
tSil
t"11
tS12
0
0


0
0
tS12
IS22
0
0


0
0
0
0
tSil
S112


0
0
0
0
tS12
t S22


(7.25)


and


[t+AtS] =


t+AtS12

S+At S22
t+AtS12


(7.26)


With these choices for the constitutive and stress state matrices the linear strain

transformation matrix becomes


['B ] = ['BO ] + ['BL ] (7.27)

where


[ 'BLO =


Oh 0 0
Ox1

0 h,0

Ox21
ax2 ax1


ah2
Ox


ah2
Ox2


(7.28)


and








ah,


Z 12

S h+1 Oh
114 2 x124


121
xl

122 9h,
Oh, ah,
121 +2 O
x21 + 122 1


Oah, h4
131 131x
Oh O h4


ahx ah, h4 ah4
hI + 1 h14
131 X2 32 X* 31 + a 32 1

(7.29


S h t ,) i,j = 1,2,3 and k = 1,2,3,4



Similarly, the nonlinear strain transformation matrix becomes


[BNL ] =


ah1
ax,
ah.


37 0
3x2
a0h,
0 O
ax,


0

0 -

0


1 0
0 0
Ox2



0 0
dx2


Since the elastic membrane model described above will be used in conjunction with

an inviscid flow model the surface traction vector is taken to be normal to the element and

equal to the pressure difference, dp, across the membrane. This pressure difference is

assumed to be uniform over each element. Consequently, the surface traction vector is given

by


[ 'BL] =


with


(7.30)


h,2

ai 0
3x1

ah2
0 a
Oh2
o0
Ox1
ah2
Ox2

0 0

0 0


ah3
0
Xh3
03x1

Oh3
3x2

0 h3
9xl
ah3,


0 0

0 0
0 0


0

0

0


0

Oh2
ax,
ah2
Ox2


0

0 -

0


0

ah3
Ox1
3Oh
Ox2


Oh4
-0



o ah4
-0

Ox1


aX2
Ox2
Ox1


0 0

0 0


(7.31)


0

ah4
dax
ah4
Ox2








rO' 0 *
[t+Atf] = 0 = 0 (7.32)
Ap -p+

where p + and p- are the aerodynamic surface pressures on the upper and lower membrane
surfaces respectively.
The transformation relating the derivatives in physical coordinates to the derivatives

in isoparametric coordinates is given by the chain rule. This transformation is

a E ax2 Ox2] -
ax1 (1 ~as x1 r r (7.33)

aX2s r j L

where


J x a x2 ax1 ax2 (7.34)
dr as as ar

The volume integrals appearing in Eqs. (7.20) through (7.23) are evaluated using
3 x 3 Gauss quadrature on the bi-unit square according to
+1 +1
()odV = t I J ( ) dr ds (7.35)
OV -1 -1

where t is the thickness of the element.
Since the element matrices were evaluated in a local element coordinate system the
element matrices must be transformed to the global coordinate system prior to assembly.
These transformations are accomplished using the conventional 1st and 2nd order tensor
transformations given by

[K] = [T]T[K][T] (7.36)

[ F] = [ T]T[ F] (7.37)


(7.38) 8)


[R] = [TIT[R]








where [ T] is the matrix of direction cosines relating the global system to the local system

(Malvern 1969). After transformation the system of global equilibrium equations is

assembled using the direct stiffness method Bathe (1984). The nodal displacement

increment residual may be defined as follows

St+AtU
Ru = c (7.39)
all dof

where c is the membrane chord and the sum is over all unconstrained degrees of freedom.

It should be noted that the linearizing approximations given by Eqs. (7.6) and (7.7)

require that an iterative procedure be used to satisfy the exact equilibrium condition given

in Eq. (7.1). The procedure described above for incrementally increasing the aerodynamic

pressure on the membrane and iterating until a converged solution is achieved provides

considerable flexibility in the choice of a solution strategy for aeroelastic problems. The

particular strategy adopted for the membrane wing model will be discussed in the section

on computational results.

A very appealing feature of the finite element formulation described above is that the

assembled stiffness matrix is symmetric and narrowly banded. This feature is largely

responsible for the usefulness of the method in solving problems where the number of

unknowns is large. This is in direct contrast to the situation encountered in the aerodynamic

problem where the assembled influence coefficient matrix is neither symmetric nor banded.

Consequently as the number of aeroelastic elements increases the aerodynamic problem

requires an increasingly larger portion of the total solution time.

7.2 The Aerodynamic Problem in Three-Dimensions

In this section a method is presented for determining the forces and moments as well

as the pressure distribution resulting from the inviscid flow over thin wings of finite span.

The fundamental assumption concerning the flow field around the wing is that it is

irrotational, and consequently the velocity field may be derived from a scalar potential.

Implicit in this assumption is that the flow is inviscid, incompressible, and free of vorticity





78



far upstream. The use of vortex filament singularities to model the lifting potential flow

around thin wings has its historical origin and justification in work by James (1972). A

description of the method in its modem form for wings of finite span may be found in Katz

and Plotkin (1991).


vortex filament FI


control point


typical panel


I I
I

* *


*00





0 i
II


trailing vortex filament




I'




----- ---------------------
I I
II

*I


* I



Il


Figure 7.3. Discretization of a elliptical planform thin wing into a number of
quadrilateral vortex filament panels.




Figure 7.3 shows a thin wing that has been discretized into a number of quadilateral

panels. A typical panel is composed of a horseshoe shaped vortex filament and a control point

where the flow tangency condition is enforced. The figure also shows the vortex filament

to extend downstream to infinity as required by the Kelvin-Helmholtz vorticity theorem

(Katz and Plotkin 1991). The existence of this trailing vortex wake is the unique feature of

a finite span wing that distinguishes it from a two-dimensional airfoil.








By modelling the wing as an assemblage of vortex singularity segments of unknown

strength and enforcing the zero normal relative velocity condition at each control point the

following set of linear algebraic equations may be formed and solved for the vortex strengths

F

t+' t [ A] [r = [ -v, n ] (7.40)

where [A ] is the matrix of vortex influence coefficients, v is the freestream velocity vector,

and n is the panel unit normal vector. If the vortex filaments and control points are located

at the local quarter and three-quarter chord points of the panel, respectively, the Kutta

condition is implicitly satisfied and no additional boundary condition is needed at the trailing

edge (James 1972, Katz and Plotkin 1991). For the sake of brevity the superscript in Eq.

(7.40) indicating the deformation state will be omitted in the following development of the

aerodynamic problem since it is clear that the aerodynamic pressures are to be evaluated in

the final equilibrium state at time t+A t.

Once the circulation strength for each panel has been determined from Eq. (7.40) the

force vector, F acting on the bound segment of each vortex filament is given by the

generalized Kutta-Joukowski theorem (Katz and Plotkin 1991) as


F = el (v x r ) (7.41)

where Q is the fluid density, I is the length of the bound segment of the vortex filament, v is

the velocity vector at the control point, and F is the circulation vector per unit length of the

bound vortex filament. Having determined the net force acting on the bound vortex filament

using Eq. (7.41), the net pressure difference acting on the aeroelastic element is then simply

given by


Ap = (7.42)


where F is the magnitude of F and a is the panel area.











F* x3


L\
a v

Xl

Figure 7.4. Resolution of far field panel forces into lift and drag components

The force acting on each panel is again determined by using the generalized

Kutta-Joukowski theorem except now, in order to account correctly for the downwash and

induced drag on the wing, the velocity that is used v*, is the sum of the freestream velocity

and the wake induced velocity (Katz and Plotkin 1991). Figure 7.4 illustrates the far field

panel force resolved into lift and drag components and are seen to be given by


L = F*cosa F sina (7.43)


D = Fsina + F cosa (7.44)

where a is the angle of attack and the far field panel force is given by


F* = (QIv* F ) (7.45)

The total lift and drag on the wing is obtained by summing up the contributions from each

bound vortex filament over the entire wing.

The aerodynamic analysis procedure described above offers a simple method for

approximating the aerodynamic characteristics of thin wings of finite span. Unfortunately

the potential flow model is a very restrictive fluid dynamic model. Many important features

of real fluid behavior such as boundary layer growth and separation are unapproachable from

a potential flow viewpoint. However the potential flow model does serve as a reliable tool

for developing a computational procedure for the aeroelastic analysis of membrane wings








of finite span. The synthesis of the plane stress quadrilateral element and the quadrilateral

vortex panel to form a quadrilateral membrane wing element is described in the next section.

7.3 The Aeroelastic Problem

The development of a computational procedure for the aeroelastic analysis of

membrane wings may be divided into two parts the development of a viable membrane

wing element and the development of a solution strategy for the coupled aeroelastic

boundary value problem. In developing such a procedure there are many options available

in terms of both the aeroelastic element technology and the solution strategy adopted for the

solution of the coupled boundary value problem. Of course, the problem formulation and

discretizations that were described in the previous sections were chosen such that a hybrid

aeroelastic element could be easily formed by combining the plane stress quadrilateral

element of section 7.1. with the quadrilateral vortex filament panel of section 7.2.

Recalling the load stepping procedure described in section 7.1 for solving the elastic

problem, it is natural to extend this methodology to a velocity stepping procedure for the

solution of the aeroelastic problem. Furthermore, since a number of iterations are required

at each load step in the elastic problem in order to satisfy equilibrium, it is again natural to

evaluate the aerodynamic pressures based on the updated wing configuration during each

iteration. This procedure guarantees that the wing configuration and the aerodynamic

loading are in equilibrium (within some prespecified tolerance) at the end of each velocity

step. However, iterating at each load step in order to achieve aeroelastic equilibrium is quite

expensive since the computation of the aerodynamic surface pressures involves the solution

of a large nonbanded system of equations. Consequently, the CPU time required to solve the

aerodynamic problem increases much more rapidly with problem size than the CPU time

required to solve the elastic problem. This situation is a direct result of the bandedness of

the stiffness matrix and the fullness of the aerodynamic coefficient matrix.

The iterative aeroelastic procedure imbedded within a velocity stepping procedure

as described above is a very general strategy that, while guaranteeing aeroelastic equilibrium








at each velocity step, may not be the most efficient or reliable strategy for solving certain

membrane wing problems. Indeed the strategy may fail to converge in some situations.

Consequently some tailoring of the strategy may be necessary in some cases. Some of these

situations will be described in the next section on computational results.

7.4 Three-Dimensional Test Cases

Before applying the procedure described in the previous sections to the analysis of

membrane wings two limiting cases of the aeroelastic boundary value problem are

investigated and the results of the computational procedure are compared with well known

analytic solutions. The two limiting cases presented in this section are the deformation of

a square, edge constrained, elastic membrane subjected to uniform pressure load and the

aerodynamic characteristics of a rigid, zero thickness elliptical planform wing of moderate

aspect ratio.


X2

(a)


computed
analytic

0.2



0.1



0
0 --0.2 0.4 1 0 -6 0.8 -
0 0.2 0.4 1 0.6 0.8
(b) 1,


F X3


Figure 7.5. Comparison of computed and analytic midpoint displacement, (b), of a
square, initially flat membrane subjected to uniform pressure load. Computed solution
used 256 square elements as shown in (a). Poisson's ratio is 0.3.

Figure 7.5 (a) shows an initially flat, square, edge constrained membrane that has

been discretized into 256 square plane stress elements. Figure 7.5 (b) shows the computed








transverse midpoint deflection of the membrane subjected to a uniform pressure load.

Results are shown in the figure for several different values of the dimensionless stiffness

parameter II1, where uniform pressure has been substituted for the stagnation pressure in the

definition of li. The analytic result for this problem is shown by the solid line in the figure

and is taken from Seide (1977). The agreement between the computed solution and the

analytic solution is seen to be quite good even with such a a small number of elements. For

completeness, all three components of the displacment field, are shown in Fig. 7.6 where the

displacement profiles are taken along lines where the x2 coordinate is constant.

Since the membrane is initially flat and free of any pretension, the tangent stiffness

matrix is singular for the problem of an initially flat membrane without pretension, and the

load stepping procedure fails for the first load step. Consequently, the solution strategy

adopted for the test case of Figs. 7.5 and 7.6 was to pretension the membrane for the first

load step in order to remove the singularity in the tangent stiffness matrix, and subsequently

remove the pretension before assembling the tangent stiffness at the next load step.

The solution shown in the Fig. 7.6 was obtained by taking 5 load steps to develop the

final pressure load and allowing iteration only at the last load step when the pressure had

reached its final value. The convergence path for this problem may be seen in Fig. 7.6 (d)

It may be seen that, once the iterative procedure begins after load step 5, the covergence is

very rapid. However, if fewer load steps were taken before iteration was allowed, the

convergence would have been much slower since the tangent stiffness matrix, at load step

2 or 3, would have been a much less accurate estimate of the stiffness of the membrane than

at load step 5. Typically, between 5 and 10 explicit load steps are adequate to insure a stable

and rapid convergence path for uniform pressure loading. However for an aeroelastic

problem, more explicit load steps may be required prior to iteration since the net pressure

acting on the membrane will then depend on the membrane configuration.





84






0.01 0.01



0.005 0.005



0-0._x 0 -- 0--



-0.005 -0 005



-0.01 -0.01
(a) 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 I
(a) x,/c (b) x,/c




0.2 le+00
X- = s/3.
X-,= S2. le-01
0.15
le-02

S0.1 Ic-03

lc-04
0.05
le-OS

le-06 ---
0 0.2 0.4 0.6 0.8 0 2 4 6 8 10
(c) x,/c (d) loadstep/iteration


Figure 7.6. Computed displacements of a square, edge restrained, initially flat
membrane subjected to uniform pressure load with 171 = 2.5. Computed solution used
256 square elements. Poisson's ratio is 0.3.



Having investigated the accuracy of the elastic formulation with a uniform pressure

load in the previous section consider now the aerodynamic properties of a thin, rigid,

cambered, elliptic planform wing in uniform flow. Figure 7.7 shows the computed

chordwise loading at several spanwise locations for an elliptical wing of aspect ratio 10 with

the xj coordinate defined by a NACA 6-series meanline (Abbott and von Doenhoff 1959).

The angle of attack was adjusted until the difference between the geometric angle and the









induced angle was equal to the ideal angle of attack for the chosen meanline. In Fig. 7.7 the

computed chordwise loading at several spanwise locations along the wing is shown along

with the two-dimensional analytic result from thin wing theory for the same NACA

meanline at the ideal angle of attack (Abbott and von Doenhoff 1959). The computed loading

on the midspan cross section may be seen to be in good agreement with the analytic result

except very near the leading edge where the potential flow model is singular.


0.8 k


0.6


0.4
0.4


0 0.2 0.4 0.6 0.8 1
X /C


Figure 7.7. Computed pressure coefficients at several spanwise locations for a rigid,
ellipitical planform wing of aspect ratio 10 operating at the ideal angle of attack
(a=2.580). The profile of the wing is defined by a NACA 6-series meanline (a=.50,
CLi=.50). The computed solution used 15 spanwise, and 30 chordwise panels.


The accuracy of the lifting surface model may be further tested by computing the lift

curve slope, CLa, for an elliptic planform wing of varying aspect ratio. A comparison is

shown in Fig. 7.8 for three different aspect ratios between 3 and 10. The analytic result is

taken from Van Dyke (1975). The agreement with the analytic result may be seen to be quite

good.


o ,=s/16
S,=s*3/16
x,=s*5/16
a x,=s/2
-- thin wing theory



- o
0


0O


o o

















4


,2




computed
S analytic
0
0 5 10 15
aspect ratio


Figure 7.8. Comparison of computed and analytic lift curve slopes for an elliptical
planform wing of varying aspect ratio. Computed solution used 24 spanwise, and 12
chordwise panels for a total of 288 quadrilateral elements.



7.5 Application to an Elliptical Planform Elastic Wing

Having established the accuracy of the elastic and aerodynamic computational

procedures independently in the previous section consider next the application of the

procedure to a flexible membrane wing. Figure 7.9 shows the computed equilibrium

chordwise loading and membrane geometry for an initially flat, pretensioned, edge

constrained, elliptic planform wing of aspect ratio 4. The controlling aeroelastic parameter

isJH2 and is taken to be equal to 2 and the angle of attack is 5. For comparison, the chordwise

loading for rigid wing of the same geometry and operating under the same conditions is also

shown in Fig 7.9 (b). As may be seen in the figure the effect of stretch in the membrane is

to move the center of pressure toward the trailing edge of the wing and to increase the lift.

The velocity stepping strategy used for this problem was to take 10 steps to develop the final

velocity and to allow iteration only at the last load step. The convergence path shown in Fig.

7.9 (d) is at the 10th load step.





87







3 3
-- x2= -s/16 = s/16
-- = s*3/16 L = s*3/16
x* = s*5/16 x= sS5/16
x1 = s/2 -- x, = s/2

2 2




11




0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8
x /c x1/c
(a) (b)


0.05 10'
I--*. x 2 s110.6
0.04

10"
0.03


0.02
S10'

0.01 10

0 10 .

0 0.2 0.4 0.6 0.8 1 0 5 10 15 20
X, iteration
(c) (d)


Figure 7.9. Equilibrium chordwise aerodynamic loading for (a), an initially flat,
pretensioned, ellipitical planform, constant tension membrane wing of aspect ratio 4, 1
= 2, HI = 0, and a = 5. and (b), a similar rigid wing at the same angle of attack, (c),
equilibrium membrane profiles, (d), residual of incremental nodal displacements.
Poisson's ratio is 0.3. Computed solutions used 15 spanwise, and 30 chordwise panels.


The effect of membrane deformation on wing lift can be seen in Fig. 7.10 where the

lift curves for both a constant tension wing and an elastic wing are shown along with the lift

curve for a rigid wing of the the same aspect ratio. The effect of membrane stretch on the

constant tension wing is simply to increase the lift curve slope by an amount inversely








proportional to 172, otherwise the variation of lift with angle of attack is identical to a rigid

wing. The behavior of the lift curve for the elastic wing is somewhat different due the

intrinsic nonlinearity in the deformation of the membrane. At very small incidence angles

the elastic wing develops lift very rapidly since the membrane cannot initially resist

deformation. However once the wing inflates and geometrically stiffens there is very little

further deformation and the lift curve slope is nearly identical to a rigid wing of the same

aspect ratio.






1.2 Hn,=o1, nr20
S n,=o, n1=2
rigid


0.8

S0.6

0.4

0.2

0-
0 2 4 6 8 10 12
a


Figure 7.10. Variation of lift coefficient with angle of attack for initially flat, ellipitical
planform, elastic and constant tension membrane wings of aspect ratio 4. Poisson's ratio
is 0.3. Computed solutions used 16 panels spanwise and chordwise.


Finally, the sensitity of the aeroelastic solution to the discretization scheme is shown

in Table 7.1. As may be seen in the table the computed lift is essentially the same for all three

discretizations tabulated. The computed drag is somewhat more sensitive to the refinement

of the grid since it is proportional to square of the lift coefficient. However the relative

insensitivity of results to the discretization is largely an artifact of the edge constrained

boundary condition used in the computations. For other edge boundary conditions, e.g., the





89


free trailing edge of a marine sail, the computed solution would undoubtedly demonstrate

a greater sensitivity to grid refinement.


Table 7.1. Effect of grid refinement on the computed lift and drag for an initially flat,
ellipitical planform, elastic membrane wing of aspect ratio 4. 11 = 10, 172 = 0, and a =
5. Poisson's ratio is 0.3.
layout of elements number of elements CL CD
8x8 64 .504 .0125
16x16 256 .503 .0119
24x24 576 .502 .0116















CHAPTER 8
SUMMARY AND CONCLUSIONS

In the present work a numerical model simulating the aeroelastic characteristics of

a flexible two-dimensional membrane wing has been presented. The use of the

Navier-Stokes equations as the fluid dynamic model in the present model is a substantial

departure from previous work on sail aerodynamics which has, almost universally, adopted

a potential based description of the flow field.

The two-dimensional aeroelastic boundary value problem was nondimensionalized

and a set of six basic dimensionless parameters was derived which govern the solution of the

problem. An additional parameter, the frequency ratio, was proposed as a meaningful

parameter for characterizing the harmonically driven unsteady aeroelastic problem.

A numerical procedure was then developed for the solution of the coupled aeroelastic

problem and was shown to yield results that are in agreement with available analytic

solutions for several appropriate limiting cases. The numerical procedure was also shown

to satisfy certain identities as dictated by the fundamental fluid dynamic conservation laws.

The role of viscosity in membrane wing aerodynamics was investigated using the

numerical model for both steady and unsteady flows. These investigations were facilitated

by distinguishing three classes of problems which are associated with corresponding limiting

cases of the dimensionless parameter set. The aerodynamic characteristics of membrane

wings at Reynolds numbers between 103 and 104 were shown to be substantially different

from those predicted by a potential based membrane wing theory.

The role of viscosity was shown to be preeminent in the harmonically forced

unsteady flow about a membrane wing. In this case the influence of viscosity is enhanced

since the acceleration and deceleration of the freestream velocity strongly influences the








separation and reattachment of the flow. The periodic appearance and collapse of

recirculation zones, along with an attendant adjustment in the membrane configuration,

results in an aeroelastic response which may not be characterized as simple harmonic

response at the freestream forcing frequency.

Finally, a three-dimensional potential based membrane wing element was derived

and shown to yield results that, for several limiting cases, are in agreement with available

analytic solutions. This finite span aeroelastic model, although it makes use of a simplified

fluid dynamic model, nonetheless serves as a useful first step in developing a viscous flow

based model for finite span membrane wings.

In terms of future work, the extension of the present two-dimensional model to

include the effects of turbulence in the flow is very appealing. This extension is especially

attractive since a direct, meaningful comparison of the numerical predictions could then be

made with available experimental data. Reasonable agreement between the measured and

predicted aeroelastic characteristics of membrane wings over the full range of aeroelastic

parameters would mark a turning point in the state of affairs described by Marchaj (1979)

in the introduction.













REFERENCES


Abbott, I., and A. von Doenhoff 1959, Theory of Wing Sections, Dover, Mineola, NY.

Bathe, K. J. 1984, Finite Element Procedures in Engineering Analysis, Prentice-Hall,
Englewood Cliffs, NJ.

Boschitsch, A. H. and T. R. Quackenbush 1993, High Accuracy Computation of
Fluid-Structure Interaction in Transonic Cascades, AIAA paper 93-0485 presented at the
31st AIAA Aerospace Sciences Meeting, Reno, NV.

Braaten, M., and W. Shyy 1986, A Study of Recirculating Flow Computation using
Body-Fitted Coordinates: Consistency Aspects and Mesh Skewness, Numerical Heat
Transfer, 9, pp. 559-574.

Chambers, L. I. 1966, A Variational Formulation of the Thwaites Sail Equation, Quarterly
Journal on Mechanical and Applied Mathematics, 19, pp. 221-231.

de Matteis, G., and L. de Socio 1986, Nonlinear Aerodynamics of a Two-Dimensional
Membrane Airfoil with Separation, AIAA Journal, 23, pp. 831-836.

Greenhalgh S., H. C. Curtiss and B. Smith 1984, Aerodynamic Properties of Two
Dimensional Inextensible Flexible Airfoils, AIAA Journal, 22, pp. 865-870.

Holla, V. S., K. P. Rao, C. B. Asthana and A. Arokkiaswamy 1984, Aerodynamic
Characteristics of Pretensioned Elastic Membrane Rectangular Sailwings, Computer
Methods in Applied Mechanics and Engineering, 44, pp. 1-16.

Jackson, P. S. 1983, A Simple Model for Elastic Two-Dimensional Elastic Sails, AIAA
Journal, 21, pp. 153-155.

Jackson, P. S. 1985, The Analysis of Three-Dimensional Sails, Proceedings of CANCAM,
University of Western Ontario.

Jackson, P. S., and G. W. Christie 1987, Numerical Analysis of Three-Dimensional Elastic
Membrane Wings, AIAA Journal, 25, pp. 676-682.

Jackson, P. S., and S. Fiddes 1993, Two-Dimensional Viscous Flow Past Flexible Sail
Sections Close to Ideal Incidence, unpublished paper.

James, R. M. 1972, On the Remarkable Accuracy of the Vortex Lattice Method, Computer
Methods in Applied Mechanics and Engineering, 1, pp. 55-79.

Katz, J., and A. Plotkin 1991, Low Speed Aerodynamics, McGraw-Hill, New York, NY.

Kroo, I. 1986, Aeroelasticty in Very Light Aircraft, Recent Trends in Aeroelasticity, P.
Hajela, Editor, University Presses of Florida, Gainesville, FL, pp. 315-321.








Leonard, J. 1988, Tension Structures: Behavior and Analysis, McGraw-Hill, New York, NY.

Malvern, L. 1969,An Introduction to the Mechanics ofa Continuous Medium, Prentice-Hall,
Englewood Cliffs, NJ.

Marchaj, C. A. 1979, Aero-Hydrodynamics of Sailing, Dodd, Mead and Co., New York, NY.

Murai, H., and S. Maruyama 1980, Theoretical Investigation of the Aerodynamics of Double
Membrane Sailwing Airfoil Sections, Journal ofAircraft, 17, pp. 294-299.

Mural, H., and S. Maruyama 1982, Theoretical Investigation of Sailwing Airfoils Taking
Account of Elasticities, Journal ofAircraft, 19, pp 385-389.

Murata, S., and S. Tanaka 1989, Aerodynamic Characteristics of a Two-Dimensional Porous
Sail, Journal of Fluid Mechanics, 206, pp. 463-475.

Newman, B. G. 1987, Aerodynamic Theory for Membranes and Sails, Progress in
Aerospace Sciences, 24, pp. 1-27.

Newman, B. G., and H. T. Low 1984, Two-Dimensional Impervious Sails: Experimental
Results Compared with Theory, Journal of Fluid Mechanics, 144, pp. 445-462.

Nielsen, J. N. 1963, Theory of Flexible Aerodynamic Surfaces, Journal of Applied
Mechanics, 30, pp. 435-442.

Ormiston, R. A. 1971, Theoretical and Experimental Aerodynamics of the Sailwing,
Journal of Aircraft, 8, pp. 77-81.

Patankar, S. V. 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington,
DC.

Patankar, S. V., and D. B. Spalding 1972, A Calculation Procedure for Heat, Mass and
Momentum Transfer in Three-Dimensional Parabolic Flows, Int. J. Heat and Mass Transfer,
15,pp. 1787-1806.

Schlichting, H. 1979, Boundary Layer Theory, McGraw-Hill, New York, NY.

Seide, P. 1977, Large Deflections of Rectangular Membranes Under Uniform Pressure,
International Journal of Nonlinear Mechanics, 12, pp. 397-406.

Shyy, W. 1994, Computational Modelingfor Fluid Flow and Interfacial Transport, Elsevier,
Amsterdam, The Netherlands.

Shyy, W, S. Thakur and J. Wright 1992, Second-Order Upwind and Central Difference
Schemes for Recirculating Flow Computation, AIAA Journal, 30, pp. 923-932.

Shyy, W., and T. C. Vu 1991, On the Adoption of Velocity Variable and Grid System for Fluid
Flow Computation in Curvilinear Coordinates, Journal of Computational Physics, 92, pp.
82-105.

Smith, R. W., and W. Shyy 1993, A Viscous Flow Based Membrane Wing Model, AIAA
paper 93-2955 presented at the 24th AIAA Fluid Dynamics Conference, Orlando, FL.








Sneyd, A. D. 1984, Aerodynamic Coefficients and Longitudinal Stability of Sail Airfoils,
Journal of Fluid Mechanics, 149, pp. 127-146.

Sugimoto, T., and J. Sato 1988, Aerodynamic Characteristics of Two-Dimensional
Membrane Airfoils, Japan Society for Aeronautical and Space Sciences Journal, 36, pp.
36-43.

Thomas, P. D., and C. K. Lombard 1979, Geometric Conservation Law and its Application
to Flow Computations on Moving Grids, AIAA Journal, 17, pp. 1030-1037.

Thompson, J., Z. Warsi, and C. Mastin 1985, Numerical Grid Generation, Elsevier,
Amsterdam, The Netherlands.

Thwaites, B. 1961, Aerodynamic Theory of Sails, Proceedings Royal Society of London,
261, pp. 402-42.

Van Dyke, M. 1975, Perturbation Methods in Fluid Mechanics, Parabolic, Stanford, CA.

Vanden-Broeck, J. M. and J. B. Keller 1981, Shape of a Sail in a Flow, The Physics ofFluids,
24, pp. 552-553.

Vanden-Broeck, J. M. 1982, Nonlinear Two-Dimensional Sail Theory, The Physics of
Fluids, 25, pp. 420-423.

Voelz, K. 1950, Profil und Luftriebeines Segels, ZAMM, 30, pp. 301-317.