UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

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 TwoDimensional Potential Flow Based Models ....... 2 1.2.2 TwoDimensional Viscous Flow Based Models ........ 4 1.2.3 ThreeDimensional 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 THREEDIMENSIONAL MEMBRANE WING MODEL ..... 70 7.1 The Elastic Membrane Problem in ThreeDimensions ...... 70 7.2 The Aerodynamic Problem in ThreeDimensions .......... 77 7.3 The Aeroelastic Problem ............................. 81 7.4 ThreeDimensional 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 NavierStokes 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 twodimensional, 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 threedimensional 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 hitormiss 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 TwoDimensional 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), VandenBroeck 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. VandenBroeck (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 TwoDimensional 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 NavierStokes 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 bodyfitted, finite volume formulation of the NavierStokes 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 ThreeDimensional Potential Flow Based Models The development of a threedimensional aeroelastic membrane wing model introduces several complicating factors in the analysis which do not have to be considered in the twodimensional analysis. In contrast to the twodimensional case where the inextensible approximation is often valid, the effects of elasticity are preeminent in threedimensional 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 threedimensional 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 twodimensional 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 threedimensional membrane wing problem. The twodimensional 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 threedimensional 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 twodimensional 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 NavierStokes 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 NavierStokes 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 threedimensional 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 twodimensional 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 bodyfitted 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 NavierStokes equations which may be written in twodimensional 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) SX1 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 NavierStokes 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 quarterchord 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 quarterchord 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 twodimensional 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 bodyfitted 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 lDO 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 ^4t 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 finitedifference 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 Ras 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.152.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 + X2,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 bodyfitted as the wing deforms. The maintenance of the bodyfitted 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 bodyfitted 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, twodimensional 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 threequarter 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 NavierStokes equations and an exact solution to the conservation laws may be obtained (Schlichting 1979). le+02 le+01 le+00 le01 le02 le03 le04 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 104 and the continuity equation converges to a terminal level of 5 x 105. These levels of residuals are consistent with the single precision floating point accuracy of the arithmetic used in the calculation. le+02 le+00 le01 I R R l e02 r R le03 Sle04 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, bodyfitted 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 107 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, bodyfitted 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(107) (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 104 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 102 was used during the solution of Eq. (3.1). I le+00 Sv, residual Sv residual le01  mass residual  x, residual le02 l, le03 le04 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, Rn4xl03. The effect of grid refinement on the computed streamlines for the inextensible membrane wing at Rn4xl03 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,1720, =.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 EO, 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 T70011 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 THREEDIMENSIONAL MEMBRANE WING MODEL 7.1 The Elastic Membrane Problem in ThreeDimensions 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 timelike 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 PiolaKirkoff 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 biunit 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 ](i1) (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 threedimensional 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 biunit 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 ThreeDimensions 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 KelvinHelmholtz 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 twodimensional 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 threequarter 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 KuttaJoukowski 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 KuttaJoukowski 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 ThreeDimensional 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 00._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. le01 0.15 le02 S0.1 Ic03 lc04 0.05 leOS le06  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 6series 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 twodimensional 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 6series 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 twodimensional membrane wing has been presented. The use of the NavierStokes 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 twodimensional 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 threedimensional 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 twodimensional 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, PrenticeHall, Englewood Cliffs, NJ. Boschitsch, A. H. and T. R. Quackenbush 1993, High Accuracy Computation of FluidStructure Interaction in Transonic Cascades, AIAA paper 930485 presented at the 31st AIAA Aerospace Sciences Meeting, Reno, NV. Braaten, M., and W. Shyy 1986, A Study of Recirculating Flow Computation using BodyFitted Coordinates: Consistency Aspects and Mesh Skewness, Numerical Heat Transfer, 9, pp. 559574. Chambers, L. I. 1966, A Variational Formulation of the Thwaites Sail Equation, Quarterly Journal on Mechanical and Applied Mathematics, 19, pp. 221231. de Matteis, G., and L. de Socio 1986, Nonlinear Aerodynamics of a TwoDimensional Membrane Airfoil with Separation, AIAA Journal, 23, pp. 831836. Greenhalgh S., H. C. Curtiss and B. Smith 1984, Aerodynamic Properties of Two Dimensional Inextensible Flexible Airfoils, AIAA Journal, 22, pp. 865870. 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. 116. Jackson, P. S. 1983, A Simple Model for Elastic TwoDimensional Elastic Sails, AIAA Journal, 21, pp. 153155. Jackson, P. S. 1985, The Analysis of ThreeDimensional Sails, Proceedings of CANCAM, University of Western Ontario. Jackson, P. S., and G. W. Christie 1987, Numerical Analysis of ThreeDimensional Elastic Membrane Wings, AIAA Journal, 25, pp. 676682. Jackson, P. S., and S. Fiddes 1993, TwoDimensional 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. 5579. Katz, J., and A. Plotkin 1991, Low Speed Aerodynamics, McGrawHill, 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. 315321. Leonard, J. 1988, Tension Structures: Behavior and Analysis, McGrawHill, New York, NY. Malvern, L. 1969,An Introduction to the Mechanics ofa Continuous Medium, PrenticeHall, Englewood Cliffs, NJ. Marchaj, C. A. 1979, AeroHydrodynamics 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. 294299. Mural, H., and S. Maruyama 1982, Theoretical Investigation of Sailwing Airfoils Taking Account of Elasticities, Journal ofAircraft, 19, pp 385389. Murata, S., and S. Tanaka 1989, Aerodynamic Characteristics of a TwoDimensional Porous Sail, Journal of Fluid Mechanics, 206, pp. 463475. Newman, B. G. 1987, Aerodynamic Theory for Membranes and Sails, Progress in Aerospace Sciences, 24, pp. 127. Newman, B. G., and H. T. Low 1984, TwoDimensional Impervious Sails: Experimental Results Compared with Theory, Journal of Fluid Mechanics, 144, pp. 445462. Nielsen, J. N. 1963, Theory of Flexible Aerodynamic Surfaces, Journal of Applied Mechanics, 30, pp. 435442. Ormiston, R. A. 1971, Theoretical and Experimental Aerodynamics of the Sailwing, Journal of Aircraft, 8, pp. 7781. 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 ThreeDimensional Parabolic Flows, Int. J. Heat and Mass Transfer, 15,pp. 17871806. Schlichting, H. 1979, Boundary Layer Theory, McGrawHill, New York, NY. Seide, P. 1977, Large Deflections of Rectangular Membranes Under Uniform Pressure, International Journal of Nonlinear Mechanics, 12, pp. 397406. Shyy, W. 1994, Computational Modelingfor Fluid Flow and Interfacial Transport, Elsevier, Amsterdam, The Netherlands. Shyy, W, S. Thakur and J. Wright 1992, SecondOrder Upwind and Central Difference Schemes for Recirculating Flow Computation, AIAA Journal, 30, pp. 923932. 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. 82105. Smith, R. W., and W. Shyy 1993, A Viscous Flow Based Membrane Wing Model, AIAA paper 932955 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. 127146. Sugimoto, T., and J. Sato 1988, Aerodynamic Characteristics of TwoDimensional Membrane Airfoils, Japan Society for Aeronautical and Space Sciences Journal, 36, pp. 3643. Thomas, P. D., and C. K. Lombard 1979, Geometric Conservation Law and its Application to Flow Computations on Moving Grids, AIAA Journal, 17, pp. 10301037. 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. 40242. Van Dyke, M. 1975, Perturbation Methods in Fluid Mechanics, Parabolic, Stanford, CA. VandenBroeck, J. M. and J. B. Keller 1981, Shape of a Sail in a Flow, The Physics ofFluids, 24, pp. 552553. VandenBroeck, J. M. 1982, Nonlinear TwoDimensional Sail Theory, The Physics of Fluids, 25, pp. 420423. Voelz, K. 1950, Profil und Luftriebeines Segels, ZAMM, 30, pp. 301317. 