Interactive Nonlinear Analysis of Wind-Loaded Pneumatic Membrane Structures by Peng Siew Han B.Eng., The University of Singapore, 1978 M.Eng., Asian Institute of Technology, 1980 A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies Department of Civil Engineering We accept this thesis as conforming to the required standard The University of British Columbia May 1986 © Peng Siew Han, 1986 CP In presenting this thesis in partial fulfilment of the requirements for an advanced degree at The University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Civil Engineering The University of British Columbia 2324 Main Mall Vancouver, British Columbia Canada V6T 1W5 May 1986 ii Abstract A nonlinear analysis of pneumatic membrane structures loaded by wind, taking into account the fluid-structure interactions, is presented here. The analysis comprises two parts: the flow part which is solved using a boundary element procedure and the structures part which is analysed with a finite element approach. These two numerical techniques are then coupled to capture the effects of thefluid-structureinteraction in the analysis. The wind is represented by potential flow and thus Laplace's equation governs. The boundary element method employed in the analysis, features moving singularities of the fundamental solutions. The singularities are placed on an auxiliary boundary located outside the domain of interest. Results from solved examples show the technique to be highly accurate for the computation of wind loads. The method is however, nonlinear with the nonlinearity occuring in the determination of the auxiliary boundary. The membrane, assumed to be perfectly flexible and to exhibit only extensional rigidity, is modeled by a small strain-finite deformation theory. A Galerkin finite element approach is formulated for the analysis of two-dimensional pneumatic membrane structures. The analysis is based on an infinitely long cylindrical membrane with an initially circular cross-section. For three-dimensional pneumatic membrane structures, the analysis is accomplished with a continuum mechanics-based incremental finite element formulation using a convected coordinate description. The metric tensor for any deformed configuration is derived and from this, a consistent measure of deformation is obtained. In both analyses, the governing finite element equations are solved via Newton-Raphson iteration with an 'automatic' line-search scheme. A pertubation technique is developed to check the accuracy of the proposed two-dimensional finite element analysis for low wind speeds. At higher wind speeds, comparison with available solutions are made. To assess the performance of the proposed three-dimensional finite element analysis, a square membrane and a cylindrical membrane are analysed. Good agreement with published results is observed. Having established the accuracies of the boundary element and the finite element methods on independent flow and structures problems respectively, interactive analysis is then carried out. Two pneumatic membrane structures are analysed: an infinitely long cylindrical membrane as a two-dimensional example, and the air-supported roof of the BC Place Stadium as a three-dimensional example. Contents iii Contents Abstract Tables Figures Acknowledgement 1 Introduction 1.1 1.2 1.3 1.4 2 Prologue Literature Review Scope of Investigation and Thesis Layout Nomenclature 1 1 3 5 8 Two and Three Dimensional Boundary Element Analyses 2.1 2.2 2.3 2.4 2.5 2.6 3 Introduction Concept of Moving Singularities Theoretical Formulation Computer Implementation Numerical Results Summary 10 10 11 12 15 16 30 Two Dimensional Finite Element Analysis 3.1 3.2 3.3 3.4 3.5 3.6 3.7 Introduction Assumptions and Theoretical Formulation Galerkin Finite Element Formulation Evaluation of Finite Element Matrices Computer Implementation Numerical Results Summary . 33 33 36 42 45 48 51 67 Contents 4 iv Three Dimensional Finite Element Analysis 4.1 4.2 4.3 4.4 4.5 4.6 4.7 5 Introduction Mathematical Preliminaries Theoretical Formulation Isoparametric Finite Element Solutions Computer Implementation Numerical Results Summary . . . . 68 69 74 84 94 100 128 Summary and Conclusions 5.1 5.2 5.3 5.4 The Flow Analysis The Structural Analysis The Interactive Analysis Recommendations for Further Study 68 130 131 132 133 134 References 135 Appendix A 140 Appendix B 144 Appendix C 146 Appendix D 155 Appendix E 164 V Tables 2.1 Location of singularities and the associated coefficients for rectangular membrane 19 2.2 Comparison of membrane deflection: Exact vs. B.E.M 20 2.3 Location of singularities and the associated coefficients forflowpast rigid circular cylinder 23 Location of singularities and the associated coefficients for heat conduction in a rectangular prism 25 Comparison of temperature and its gradient for heat conduction in a rectangular prism 25 2.4 2.5 2.6 Location of singularities and the associated coefficients forflowpast rigid sphere 30 3.1 Accuracy of hoop stress computations 55 3.2 Comparison of hoop stress predictions 58 3.3 Discretization schemes 62 4.1 Fortran variable names used in the main data input 95 4.2 Comparison of deflection response for the square membrane . . . . 102 4.3 Comparison of tension increments for the square membrane . . . . 104 4.4 Comparison of hoop stress resultants for Example 1 107 4.5 Comparison of hoop stress resultants for Example 2 108 4.6 Discretization schemes .112 4.7 Stress resultants at the apex of the roof, for the inflation problem .114 4.8 Stress resultants at the apex of the rigid roof (V = 25 m/s) .119 4.9 Maximum roof elevations (initial height=28.81m) 120 4.10 Stress resultants at the apex of theflexibleroof 125 4.11 Maximum shear stress resultants and their locations 128 vi Figures 1.1 Wind loads on a pneumatic membrane 2.1 Problem definition and auxiliary boundary 13 2.2 Rectangular membrane 17 2.3 Discretization scheme and initial singularity locations for rectangular membrane, • Boundary Node, x Singularity 18 Loci of singularity movement in the direction of arrows, for rectangular membrane (domain of problem shown shaded) 19 Discretization scheme, initial singularity locations, flow domain and boundary conditions for cylinder, • Boundary Node, x Singularity 21 Loci of singularity movement in the direction of arrows, for flow past circular cylinder (domain of problem shown shaded) 22 2.7 Pressure coefficient, C on a rigid circular cylinder 23 2.8 Discretization scheme and initial singularity locations for rectangular prism, • Boundary Node, x Singularity 24 Loci of singularity movement in the direction of arrows, for heat conduction (domain of problem shown shaded) 26 2.10 Problem definition and coordinate system forflowpast rigid sphere. 27 2.11 Boundary elements for sphere and theflowdomain 29 2.12 Pressure coefficient on a rigid sphere 31 3.1 Wind loads on a cylindrical membrane 34 3.2 Positive stress resultants 37 3.3 Deformation of an initially circular string 41 3.4 Typical membrane element . 43 3.5 Structure of INFEAM2 . 49 3.6 Membrane roof subjected to low wind pressure 55 2.4 2.5 2.6 2.9 p 1 Figures vii 3.7 Comparison of membrane deflection 56 3.8 Membrane deformations for Example 1 (£„ = 1.91, C = \ x 10 ) . . . 59 3.9 Membrane deformations for Example 2 (<f = 3.40, C = ± x IO ) . . . 59 3.10 Flow domain and boundary conditions for cylindrical membrane . . 61 3.11 Membrane deformations at each global iteration ( £ = 0.276, C = -4 -4 a . a £xl0" ) 63 4 3.12 Membrane deformations at each global iteration ( £ = 1.875, C = 0 \ x IO" ) 63 4 3.13 Pressure loads acting on membrane in equilibrium (C = \ x 10 ) 3.14 Influence of the loading ratio, <f on membrane deflections (C = -4 . . . 0 J x 10~ ) 65 4 3.15 64 Variation of membrane stress resultants N , with internal pressure, p (C = i x IO" ) . . a 66 4 a 4.1 Coordinate system for a surface in three dimensional space 4.2 Configurations in the path of deformation 74 4.3 Cauchy's triangle 76 4.4 Typical finite element 85 4.5 Structure of NFEAM3 96 4.6 Storage scheme used for global stiffness matrix 98 4.7 Discretization scheme for the square membrane 101 4.8 Deflected profile for the square membrane 103 4.9 Deflection response of the centerline with increasing load increments: square membrane 103 Variation of central deflection with internal pressures for the square membrane 105 Variation of maximum tension increments with internal pressures for the square membrane 105 4.12 Discretization schemes for the cylindrical membrane 106 4.13 Membrane deformations for Example 1 (<f = 1.91, C = ^ x l O ) 4.10 4.11 . -4 a . . . . 70 107 Figures viii 4.14 Membrane deformations for Example 2 (£„ = 3.40, C - | x l O ) 4.15 Stadium and flow domain 4.16 BC Place Stadium, Vancouver, Canada 4.17 Undeformed and initial roof profiles due to internal pressure 4.18 Stress resultants along the minor axis for the inflation analysis 4.19 Model test of BC Place Stadium (Parkinson (1980)) 116 4.20 Pressure and roof profile for the minor axis 117 4.21 Pressure and roof profile for the major axis 118 4.22 Stress resultants along the minor axis for theflowmodel investigations 119 Pressure and roof profile: (A) and (B) refer to the minor axis, (C) and (D) refer to the major axis (V=15 m/s) 121 Pressure and roof profile: (A) and (B) refer to the minor axis, (C) and (D) refer to the major axis (V=25 m/s) 122 Pressure and roof profile: (A) and (B) refer to the minor axis, (C) and (D) refer to the major axis (V=35 m/s) 123 Membrane deformations along the Minor Axis for varying loading ratios 124 4.27 Deformed profiles (internal pressure = 345 Pa) 126 4.28 Stress resultants along the minor axis: (A) for V=15m/s, (B) for V=25m/s, (C) for V=35m/s 127 D.l The two line search conditions 158 D.2 Quadratic backtrack 160 D. 3 Cubic backtrack E. l Finite element grid of BC Place Stadium Roof 4.23 4.24 4.25 4.26 - 4 . . . 108 . 110 I l l . . . . . . 113 114 . 161 165 ix Acknowledgement The author wishes to express his immense gratitude and thanks to his advisor, Dr. M.D. Olson, for his invaluable advice and assistance in the preparation of this thesis. He would also like to thank Dr. R.L. Johnston for his collaboration in the boundary element research and to Dr. H. Vaughan for reviewing the section on the continuum mechanics of large deformations. The assistance of Mr. A.G. Martin with the color graphics and Mr. Fon Chow for providing free drafting offigures,too difficult to produce on a computer, are also much appreciated. Thefinancialsupport of the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Finally, the author would like to express his gratitude to his wife, Pat and children, Flore and Merv for their support and the frustrations they have endured with great fortitude. CHAPTER 1 Introduction 1.1 P R O L O G U E Membrane structures, commonly use air pressure differential to stiffen or stabilize the membrane, resulting in a structural shape as shown in Figure 1.1. The structural envelope, being supported by air, permits large areas to be spanned with the minimum amount of F I G U R E 1.1 Wind Loads on a Pneumatic Membrane. l Introduction / l.l 2 material. It combines the inherent strength of the material acting as a membrane in tension with the structural efficiency of a shell. Hence, pneumatic membrane structures 1 are probably the most efficient structural forms available to date. In recent years, this has resulted in an increase in the construction of very large pneumatic membrane structures such as air supported roofs for stadiums. A major design problem for these structures is their response to moving loads like those due to wind. It is not uncommon for pneumatic membrane structures to deform visibly under wind loads and in fact, many of them have been severely damaged by wind. The main difficulty here is to obtain an accurate description of the action of wind on these structures. Designers generally rely on solutions based on rigid versions of these structures. In this approach, the wind load is described by pressure coefficients obtained empirically in wind tunnel experiments and used together with the known velocity field of air in free space. The structural deformation, however large, i 3 assumed to have no effects on the aerodynamic pressure variation. In other words, the wind loads are computed as though the structure is perfectly rigid. This may result in the wind loads being grossly underestimated. Sophisticated computer programs based on such rigid models and incorporating both static and dynamic analyses have been developed (Morris 1975). A flexible membrane as explained, will deform quite substantially and the large structural motion must distort the wind velocity field. Thus an accurate determination of wind loads is only possible by resorting to the interaction between the structure and the wind. In the state-of-art report on pneumatic membrane structures by the ASCE Task Committee on Air Supported Structures (1979), it concluded that the behavior of pneumatic membrane structures subjected to wind loads can be significantly different from that of rigid structures. From this point of view, one can say that the literature presently available is inadequate. It is the objective of this thesis to make a contribution in overcoming this deficiency by formulating an interactive structural analysis of pneumatic membrane structures loaded by wind. The material from which pneumatic membranes are made usually possesses little or no bending stiffness. Introduction 1.2 / 1.2 L I T E R A T U R E 3 R E V I E W The British engineer, Frederick Lanchester (1917) appears to have been the first to introduce the concept of air supported structures. The basic principles of pneumatic membrane structures were proposed in his patent for a field hospital. In a lecture on long span structures, Lanchester (1938) mentioned the potential of pneumatic membrane structures for enclosing large areas like those in aircraft hangers and sports arenas. Unfortunately, Lanchester was not able to realize the achievements of his designs. This was attributed to the lack of suitable materials and construction techniques necessary to accomplish these proposals (Dent 1971). During the Second World War, attempts to advance the developments of pneumatic membrane structures were made by Herbert Stevens (1942), an American engineer, with a design proposal to house the manufacture of aircraft. The fulfilment of Lanchester's ideas came in the late forties when Walter Bird (1967) built many pneumatic membrane structures, mainly of spherical or cylindrical shapes. Frei Otto (1962) presented in his book several very imaginative pneumatic forms and applications. By the late sixties with an increasing interest generated by these developments, the first International Symposium on pneumatic membrane structures was organised in Stuttgart, Germany in 1967. The construction of very large pneumatic membrane structures originated with the air-supported roof by David Geiger (1970) at the 1970 Exposition in Osaka, Japan. Since then, several large stadiums with air-supported roofs have been built in North America. Take for example, the BC Place Stadium in Vancouver, British Columbia, was completed in 1984. Such large structures fully exploit all the advantages associated with pneumatic membrane structures. The first significant analytical treatment of the finite deformation analysis of membranes was presented by Adkins and Rivlin (1952) using the neo-Hookean and Mooney forms of strain energy function. Further work on this classical approach based on membrane shell theory has been given in Klingbeil and Sheild (1964), Hart-Smith and Crisp (1967) and Green and Adkins (1970). Introduction / 1.2 4 With the advent of high-speed computers, numerical methods based on finite element modelling became the main tools of analysis. Becker (1966) solved the problem of finite in-plane deformation of elastic sheets subjected to prescribed conditions. Oden (1966, 1967), Oden and Kubitza (1967), Oden and Sato (1967) presented several solutions on the finite deformation analysis of elastic membranes using constant strain triangular elements. An excellent summary is given in Oden's book (1972). Employing a convected coordinate description approach, Haug and Powell (1972a) analysed problems of large displacements but small strains, using warped isoparametric quadrilateral membrane elements with straight sides. Dynamic response investigation of membranes have been reported by Oden, Key and Frost (1974) using constant strain triangular elements, by Benzley and Key (1976) using cubic quadrilateral elements and by Morris (1981) using a finite difference approach. Morris (1981) also attempted a random vibration analysis of flexible roofs in wind and approximated the aerodynamic pressure by a stationary Gaussian process. Stability analyses of pneumatic membrane structures have been presented by Malcolm and Glockner (1978, 1981), Ahmadi and Glockner (1983, 1984) and Lukasiewicz and Glockner (1983). Their investigations are however, restricted to ponding instability of cylindrical and spherical air supported membranes. Ross (1970) and Uemura (1972) reported some results on the large membrane deformations of cylindrical pneumatic structures under wind loads. Their analyses are however, based on rigid structure modelling and hence there is no wind-structure interaction. In the field of interactive analysis of wind-loaded pneumatic membrane structures, very few published results are available. In an attempt to overcome this deficiency, Kunieda, Yokoyama and Arakawa (1981) presented some solutions obtained from an interactive analysis of cylindrical membrane structures subjected to wind. They introduced simple flow models capable of reproducing flow separation and wake formation to predict the wind loads. Their results indicated a maximum difference of 4% in the deflection response between a real fluid model and an inviscid fluid model of the wind. Introduction 1.3 5 / 1.3.1 S C O P E O F INVESTIGATION A N D THESIS L A Y O U T As mentioned earlier, the objective of this thesis is to develop a method of analysis for wind-loaded pneumatic membrane structures, taking into account thefluid-structureinteractions. This proposed interactive analysis is carried out in two parts: solution to the flow part and solution to the structures part. A boundary element scheme is used to solve the flow problem, while a finite element procedure is employed in the solution of the structures problem. The performances of both these techniques are assessed on independent flow and structures problems with known solutions. Having established satisfactory accuracies, the boundary element program and the finite element program are then coupled to produce the interactive results. This is done as follows: The wind loads are computed based on an assumed membrane profile in the flow analysis. With the loading known, a new membrane profile is then obtained from the structural analysis. We then go back to the flow analysis to compute new wind loads based on the new profile. This iterative process is continued until an equilibrium profile consistent with the wind loading is attained. Two and three-dimensional pneumatic structures are analysed in this manner. The thesis layout is arranged in this manner: a) Chapter I — Introduction b) Chapter 2 — Two and Three-Dimensional Boundary Element Analyses c) Chapter 3 — Two-Dimensional Finite Element Analysis d) Chapter 4 — e) Chapter 5 — Summary and Conclusions 1.3.1 Three-Dimensional Finite Element Analysis T H E F L O W ANALYSIS Here we are interested in computing the wind loads acting on the membrane for any deformed configuration. To model the wind, potential flow is assumed and hence, Laplace's equation governs. The boundary element scheme used, features moving singularities of the Introduction / 1.3.2 6 fundamental solution. This method was first introduced by Mathon and Johnston (1977) and a subsequent variation was reported by Han, Olson and Johnston (1984). Unlike conventional boundary element methods where the singularities are located on the boundary of the problem, the singularities here, are placed on an auxiliary boundary located outside and encompasses the domain of interest. This is done to avoid singular integrals and to improve on the accuracy of the computed derivatives. Also, by allowing the singularities to move, we realized yet another advantage. Unlike conventional boundary element methods which require the same number of singularities as the number of boundary elements used, here the number of singularities required can be significantly reduced. The determination of the auxiliary boundary is done via the Galerkin minimization criterion. The discrete least square criterion has also been used to generate the auxiliary boundary. See for example, Johnston, Fairweather and Han (1984). Thus by allowing the singularities to move freely until they attain optimum positions, a highly adaptive and accurate method for computing wind loads is realized. However, the price one pays for this feature is that the resulting boundary element equations are nonlinear, with the nonlinearity occuring in the determination of the coordinates of the singularities. This nonlinearity is not seen as a major handicap as good algorithms for solving systems of nonlinear equations are readily available (see Powell (1970), Hopper (1973), Blue (1976), More (1977) and Blue (1978)). The details of the foregoing formulation are reported in Chapter 2. Some of the examples used in the numerical experimentation are: deflection of a stretched membrane, flow past a rigid circular cylinder, heat conduction in a rectangular prism and flow past a rigid sphere. Their analytical solutions served as the basis of comparison for the boundary element results. The agreement obtained is excellent. 1.3.2 T H E S T R U C T U R A L ANALYSIS In this section, we are interested in computing stresses and deformations in the membrane. The analysis is based on a small strain-finite deformation theory. The membrane, whose Introduction / 1.3.2 7 self-weight is considered small compared to the internal pressure, is assumed to be perfectly flexible and exhibits only extensional rigidity. The differential equations of the structural analysis are then solved via a finite element technique using Newton-Raphson iteration with line search. The line search parameter encountered ranges from 0.2 to 0.6. The structural analysis can be conveniently grouped into two categories, namely, two and threedimensional pneumatic membrane structures. In Chapter 3, the solution of two-dimensional pneumatic membrane structures are reported. An infinitely long pneumatic cylindrical membrane whose cross-section is initially circular and subjected to a uniform wind blowing from the broadside is considered. This results in a second-order integrodifferential equation which is solved via a Galerkin finite element procedure as shown in Olson and Han (1984). To check the accuracy of the finite element results at low wind speeds, a perturbation technique is developed. The perturbation parameter, c used here, is defined as the ratio of a characteristic membrane displacement to the diameter of the cylinder. For moderate to high wind speeds, the finite element results are compared to the solutions of Uemura (1972). Good agreement is observed in both cases. To avoid unnecessary complications, the fluid-structure interactions have not been included in these checks. In Chapter 4, we describe the analysis of three-dimensional pneumatic membrane structures. Solution to the structural deformations is obtained through a contiuum mechanics-based incremental finite element formulation using a convected coordinate description. In this description, the coordinate system displaces and deforms with the body such that the coordinates of a point on the body remain unchanged. The metric tensor for any deformed configuration is derived and from this, a consistent measure of deformation is obtained. Since the deformation, however large, is assumed to be of the small-strain type, linear constitutive law becomes applicable. Expressions for the elastic and geometric stiffness matrices and the unbalanced loads are presented. An unsymmetric stiffness matrix arising from the follower-load nature of the external loading has also been incorporated into the analysis. Introduction / 1.4 8 The membrane element used is a warped isoparametric quadrilateral element given by Haug and Powell (1972a, 1972b). This element formed by straight boundaries, is warped in space. It has 12 degrees of freedom and is defined with respect to a nonorthogonal convected coordinate system. To assess the accuracy of the porposed formulation, a square membrane and a long cylindrical membrane are analysed. Here again, for simplicity, the fluid-structure interactions are ignored. The solutions are checked against published results and good agreement is observed. 1.3.3 T H E I N T E R A C T I V E ANALYSIS By coupling the two numerical methods in the manner discussed earlier, we generated some interactive results. All results pertaining to the two-dimensional pneumatic membrane structures are reported in Chapter 3, while those belonging to the three-dimensional pneumatic membrane structures are given in Chapter 4. An infinitely long cylindrical membrane is considered in the former while the air-supported roof of the B C Place Stadium, in the latter. They are also discussed in Olson and Han (1984), and Han and Olson (1986) respectively. Since the accuracies of the boundary element method and the finite element method are assessed on independent flow and structures problems respectively and good agreement is attained, we believe that the coupled results are just as accurate. 1.4 N O M E N C L A T U R E Commonly used abbreviations like i.e.(that is), l.h.s.(left hand side) and r.h.s.(right hand side) are adopted in this thesis. All other notations will be defined where they first appear. In the two dimensional formulation where a rectangular cartesian coordinate system is adopted, the notations will closely follow those given in Flugge (1973). Thus N , N^, N f, x denote membrane stress resultants; c xi x< s^, 7 ^ for membrane strains and so on. In the three X( dimensional formulation where a convected coordinate description is used, the notations will follow as close as possible those given in Green and Zerna (1968). Hence Greek indices Introduction / 1.4 9 Hence Greek indices are reserved for surface quantities and only take the values 1, 2 while all lower case Latin indices take the values 1, 2, 3 and all upper case Latin indices, the values 1, 2, 3,•••, 12. In addition, the upper case calligraphic indices will take values 1, 2, 3, 4. The summation convention for repeated indices applies unless indicated otherwise. Also it should be pointed out that since different types of indices sum to different range of values, it does not make sense to sum over repeated indices of different types. To illustrate the use of different types of indices, consider the following examples: dr a d6 a a & r = a dO + a d9 x 2 2 l = * Si + J Si + ! S3 M £ / =Af e + M 6 + M 6 + - + Af 7 1 2 3 1 N N 12 £12 y\ = N y\ + N y^ + N yj + N y\ 1 2 3 A One problem that was encountered in presenting the continuum mechanics relations for the finite deformation analysis was the need of an effective notation. Since we are dealing with so many quantities simultaneously, it is essential that the symbols used, should display all the necessary information in a compact and easy to read manner. Thus the following nomenclature has been adopted. A left superscript, symbolized exclusively by i/, indicates the configuration in which a quantity is referred to, while a right superscript and subscript would denote a tensorial quantity unless mentioned otherwise. A zero left superscript is reserved for the undeformed configuration. Incremental quantities are shown without left superscripts. For example, o »• j J This quantity occurs in the undeformed configuration 0. The indices i,j 1 as shown denote the quantity is a mixed second order tensor. This quantity occurs in the deformed configuration 1. The indices i,j as depicted show the quantity is a covariant second order tensor. This is an incremental quantity. The indices indicate the quantity is a contravariant second order tensor. CHAPTER 2 Two and Three Dimensional Boundary Element Analyses 2.1 I N T R O D U C T I O N The boundary element method is now a well-established technique for the approximate solution of differential equations. The method employs a fundamental solution of the differential equation to reformulate the problem as an integral equation on the boundary. The generation of this boundary integral equation is usually accomplished by locating the singularities of the fundamental solution on the boundary. The equation is then discretized by boundary elements and the resulting finite system of singular integrals is evaluated to yield a system of algebraic equations. This procedure of locating the singularities on the boundary presents some serious drawbacks. For example, special attention is required to evaluate accurately the singular integrals in the neighborhood of the singularities. Another is that, the accuracy of the computed derivatives on the boundary is generally very poor. A simple solution would be to ensure that the boundary nodes and the singularities are not located on the same boundary. Kupradze (1967) proposed locating the boundary nodes on an auxiliary boundary enclosing the original boundary. This approach however, has some difficulties. Christiansen (1967) has shown that with this method, the solution may not be unique and that the location of the auxiliary boundary must be chosen carefully to yield accurate results. Oliveira (1968) suggested moving the singularities rather than the boundary 10 Two and Three Dimensional Boundary Element Analyses / 2.2 11 nodes onto the auxiliary boundary. Patterson and Sheikh (1981) also presented similar proposals. There is however, a difference between these two suggestions. In the former, the integration of the kernel function is carried out with respect to the auxiliary boundary while in the latter, the integration is with respect to the original boundary. Heise (1978) investigated the numerical behavior of Kupradze's and Oliveira's methods. He concluded that care must be exercised in the choice of the auxiliary boundary, just as Christiansen had found with Kupradze's method. In a later paper, Patterson and Sheikh (1982) also investigated the best location of the auxiliary boundary. In all these approaches, the singularities are fixed in location. This results in a system of linear algebraic equations. Also the number of singularities coincide with the number of boundary nodes. In higher dimensional problems, this can lead to numerical difficulties due to the very large system of equations to be solved. 2.2 CONCEPT OF MOVING SINGULARITIES Having seen that there exist a best location for the auxiliary boundary that yields accurate results, two questions naturally arise: how can this best location which varies with different problems be established without making several guesses and more crucially, for problems where analytical solutions do not exist, how can this best location be measured? The obvious solution is to replace this idea of fixed location of the singularities with one that is automatically determined as part of the process of solving the problem. In short, a moving auxiliary boundary containing the singularities. By allowing the singularities to move freely until they attain optimum positions, we end up with a highly adaptive and accurate method. Also as shown later, the number of singularities can be significantly less than the number of boundary nodes with no apparent deterioration in accuracy. This is an important advantage in higher dimensional problems. However, the price one pays for this feature is that the system of algebraic equations will be nonlinear, with the nonlinearity occurring in the determination of the coordinates of the singularities. Two and Three Dimensional Boundary Element Analyses 12 / 2.3 Mathon and Johnston (1977) were the first to introduce this concept of moving singularities. They used a discrete least squares technique to minimise the boundary residue and obtained very good results. Further examples are found in Ho-Tai, John- ston and Mathon (1979), Johnston and Fairweather (1983), and Johnston, Fairweather and Han (1984). In the next section, the general boundary element formulation for the /-dimensional space are presented. Unless mentioned otherwise, the boundary residue here is minimized on the basis of the Galerkin criterion. The innovation introduced in the Galerkin criterion is that, we include as weighting functions, the derivatives of the approximate solution with respect to the singularity coordinates as well as the usual ones with respect to the coefficients. This derivation is taken from Han, Olson and Johnston (1984). 2.3 T H E O R E T I C A L F O R M U L A T I O N Consider a problem defined in Figure 2.1 with its differential equation represented by a linear operator i , such that, L{u{x)} = 0 xeD (2.1) with boundary conditions B{u{x)} = 0 xeT (2.2) where B is an operator of one of the following forms, { a(x) + u(x) (Dirichlet) a(x) + du(x)/dn (Neumann) a(x) + /?(x)u(x) + 7(x)c5u(x) fdn (Mixed) (2.3) where a, /3, and ^ are prescribed functions. Observe in Figure 2.1 that the domain of the problem is enclosed by a fictitious boundary which we will call, the auxiliary boundary. It contains the singularities of the fundamental solution. Its shape is allowed to vary as the solution process marches on, thereby yielding optimum locations for the singularities. The criterion that governs the 'motion' of the auxiliary boundary will now be described. The proposed formulation involves selecting a function that identically satisfies the differential Two and Three Dimensional Boundary Element Analyses ^Auxiliary F I G U R E 2.1 / 13 2.3 Boundory Problem Definition and Auxiliary Boundary. equation but not the boundary conditions. The free-space Green's function exhibits this type of behavior and will thus be employed here. Let £(tj,x) denote the free-space Green's function. Then a solution to the differential equation in Eq. (2.1) can be expressed as, N u(jO = £ C y f ( « , x ) j'=l y t^D,T (2.4) where Cj and tj are unknown coefficients and locations of the singularities respectively. The boundary conditions are then enforced in a weighted residual sense so as to yield the best solution for a given problem. To achieve this, we compute the unknown coefficients Cj and the unknown locations of the singularities tj such as to minimize the error on the Two and Three Dimensional Boundary Element Analyses 14 / 2.3 boundary which is denned by, SR = Bu ± 0 (2.5) B For convenience, let Z k = (C , t ) k be the vector of unkowns. Then discretizing the T k boundary with M boundary elements and minimizing the boundary residue fk yields, B A/ ]T J[w x ]idr B = o f[w Bu}idT = 0 k or M JT ,=i k ff where T,- is the boundary for the I l<k<N (2.6) element, and the weighting functions w are given th k by, (d(Bu)/dZ for least square minimization \ du/dZ for Galerkin minimization k k k where d dZ k \dC ' k dt / k Note that the usual Galerkin weighting functions have been augmented with functions defined as derivatives with respect to the singularity locations. Substituting the Galerkin weighting functions into Eq. (2.6) we have, l/IKSH.^ 0 - (2 8) Eq. (2.8) represents a system of (I + 1)N nonlinear boundary element equations for the unknown coefficients and singularity locations, where / represents the spatial dimension of the domain D. If integration is omitted in Eq. (2.8) and the least squares weighting function is used, we will obtain the boundary element equations associated with the discrete least square minimization technique proposed by Mathon and Johnston (1977). Two and Three Dimensional Boundary Element Analyses 2.3.1 SPECIALIZATION T O / 2.4 LAPLACE'S 15 E Q U A T I O N For simplicity, we consider potential problems. This is consistent with the assumption used in the flow analysis in which the wind is modeled as an inviscid fluid. The differential equation for the problem is thus the Laplace's equation and its normalized free-space Green's function is given by, c(t-,x) = —In (\=^—r— in 2-dimensions: U+W U j in /-dimensions: ${tj,x) (2.9) = |t - x| ' - (R + |t |) "', 2- y 2 y />3 where R denotes the radius of a fictitious circle or sphere circumscribing the domain D. It has been found by Ho-Tai, Johnston and Mathon (1979) that the addition of the term (R+ \ tj\) to normalize the free-space Green's function helps to accelerate the convergence rate in the solution of the nonlinear boundary element equations. These equations are given in the Appendix A. 2.4 C O M P U T E R I M P L E M E N T A T I O N Like the finite element method, the boundary element mesh should be such that the boundary T, is sufficiently well discretized. The number of singularities used can be very much smaller than the number of boundary nodes employed in the discretization without any apparent loss of accuracy. This important advantage can be attributed to the fact that the singularities are allowed to move. In contrast, the conventional boundary element methods which have fixed singularities, require equal number of singularities and boundary nodes. It has been found through numerical experimentation that the following relation governing the relative number of singularities and boundary nodes to hold: M Net— (for 2-D) and M N w —, /> 3 (2.10) If the number of singularities employed greatly exceed the empirical relation of Eq. (2.10), it has been observed that the singularities actually migrate into the domain D of the Two and Three Dimensional Boundary Element Analyses / 2.5.1 16 problem. While no deterioration of accuracy of the results is likely, this practice is not 1 encouraged. Firstly, the cost of computation increases quite significantly and secondly, singular solutions exist in the domain D. The following alternating algorithm is used to solve the nonlinear boundary element equations. The initial locations of the singularities are prescribed in some fashion. A good guess would be to locate them close to the boundary and distributed out evenly. With the positions of the singularities known, the associated coefficients are computed via a linear solver. These coefficients are then fed into a nonlinear solver to compute the new locations of the singularities. This back and forth iteration between the linear and nonlinear solvers is continued until convergence is achieved. 2.5 N U M E R I C A L RESULTS To assess the numerical response of the proposed formulation, these examples are used: a) Deflection of a stretched membrane b) Flow past a rigid circular cylinder c) Heat conduction in a rectangular prism d) Flow past a rigid sphere Dirichlet boundary conditions are prescribed in Example a while Examples b and d have only Neumann boundary conditions. Mixed boundary conditions are used in Example c. Exact solutions are available in these examples and will be used as the basis of comparison. 2.5.1 D E F L E C T I O N O F A S T R E T C H E D M E M B R A N E Consider as an example involving the Dirichlet boundary conditions, a homogeneous membrane stretched on a rectangular base by a uniform tension and subjected to a uniform may even result in a slight improvement of accuracy Two and Three Dimensional Boundary Element Analyses / 17 2.5.1 y F I G U R E RecfcanguJar Membrane. 2.2 lateral pressure, as depicted in Figure 2.2. The differential equation and boundary conditions for the deflection z are, dz dx 2 2 z= 0 + <Pz dy ~ _ = 2 q ( x = ±a; on 2 n ) y = ±6 where q is the ratio of the uniform pressure to the uniform tension in the membrane. Eq. (2.11) is transformed into Laplace's equation by defining, « = * + Uix + y ) 2 (2.12) 2 The differential equation and boundary conditions in terms of u are then given by, 3 u dx 2 « = \q(x 2 + y) 2 du 2 2 + on dy 2 = 0 x = ±a; (2.13) y = ±b Two and Three Dimensional Boundary Element Analyses F I G U R E 2.3 / 2.5.1 18 Discretization Scheme and Initial Singularity Locations for Rectangular Membrane, • Boundary Node, x Singularity. The solution to Eq. (2.13) is, 16a g 2 ^ (_i) ! cosh(n7r«/2a)l 1 n=l,3,5., /nirx\ , . r?—T7TTT cos ( — - J + \q(x cosh(nrr6/2a)J V 2a / ^ 4 , ? 2 . x + y ) (2.14) ' ' 2 y y The following values are used in the problem: o = 1, b = 5 , q = 2. Symmetries about the x and y axes are exploited. The problem is analysed using only three singularities located outside the domain as shown in Figure 2.3. The initial and final locations of the singularities and their coefficients are tabulated in Tabie 2.1. Since the proposed formulation exhibits moving singularities, it is interesting to trace the loci of their movement as the solution approaches convergence. This is shown in Figure 2.4. It is seen that two of the singularities did not move very much whereas the third one moved dramatically. A l l three showed little motion as the solution becomes nearly converged. Two and Three Dimensional Boundary Element Analyses / 2.5.1 19 o X-Axis Loci of Singularity Movement in the Direction of Arrows, for Rectangular Membrane (Domain of Problem shown Shaded). F I G U R E 2.4 Initial Values Table 2.1 F i n a l Values X y Coefficient X y Coefficient 1.10 1.15 0.40 0.10 0.65 0.60 -0.0045 -0.1145 0.0473 0.08 1.21 0.27 3.26 0.69 0.89 0.1552 -0.1760 0.0380 Location of Singularities and the Associated Coefficients for Rectangular Membrane Two and Three Dimensional Boundary Element Analyses Location / 2.5.2 Exact Solution 20 2-D BEM Solution Boundary Points Bl(0.0,0.50) 0.00 0.0006 B2(0.67,0.5) 0.00 -0.0008 B3(1.0,0.25) 0.00 0.0006 B4( 1.0,0.0 ) 0.00 -0.0008 Interior Points Table 2.2 11(0.25,0.20) 0.18613 0.18604 12(0.75,0.20) 0.11387 0.11391 13(0.50,0.30) 0.12707 0.12722 14(0.75,0.40) 0.05258 0.05258 Comparison of Membrane DeSections: Exact Vs. BEM The deflection computed at the boundary and interior points are given in Tabie 2.2. Very good agreement with the exact solutions is observed. Using this membrane deflection analogy, one could similarly analyse the Saint Venant torsion problems. 2.5.2 FLOW PAST A RIGID CIRCULAR CYLINDER As an example using the Neumann boundary conditions, the potential flow past a rigid circular cylinder is chosen. The flow over such a cylinder in an infinite domain is approximated by afinitedomain as depicted in Figure 2.5. Symmetry about the x-axis is invoked. The discretization scheme, initial locations of the singularities and the boundary conditions of flow are also shown in Figure 2.5. Since only the Neumann boundary conditions are prescribed, the computed potential function u, is accurate to an arbitrary constant. Hence to get meaningful results, we will obtain quantities that are given by the derivatives of the potential function. One such quantity is the pressure coefficient on the cylinder. The Initial and final locations of the singularities and their coefficients are depicted in Tabie 2.3. The movement of singularities with solution iteration is given in Figure 2.6. In Two and Three Dimensional Boundary Element Analyses / 2.5.2 21 Jk y X X X dn \ du_ ^ dn du 2 V / x dn ^ - R = I.O ' X A \ 0 x X -« F I G U R E 2.5 IC .0 » f » , V u=0 X X X 10.0 * X Discretization Scheme, Initial Singularity Locations, Flow Domain and Boundary Conditions for Cylinder, • Boundary Node, x Singularity. Figure 2.6(A), the motion of the singularities inside the cylinder is depicted, while Figure 2.6(B) shows the movement of the singularities outside the flow domain. Again, we see minimum movement as the solution is very near convergence. A comparison of the pressure coefficient on the cylinder with the exact solution is shown in Figure 2.7. The exact pressure coefficient is given by, 2 C {9) = l - 4 s i n 0 2 p (2.15) The agreement with the exact solution is excellent. Apparently, the boundary element representation of the finite domain gives an excellent approximation of the exact infinite domain problem. Also, since the pressure is given by the derivative of the potential function, we see that the derivatives on the boundary are very accurately computed. The exact pressure coefficient is computed assuming an infiniteflowdomain. Since the diameter of our cylinder is very small compared to the dimension of ourfinitedomain, we can justifiably make the comparison. Two and Three Dimensional Boundary Element Analyses F I G U R E 2.6 / 2.5.2 22 Loci of Singularity Movement in the Direction of Arrows, for Flow Past Circular Cylinder (Domain of Problem shown Shaded). Two and Three Dimensional Boundary Element Analyses / 2.5.2 Initial Values 23 Final Values X y Coefficient X y Coefficient 11.00 11.00 8.00 -8.00 -11.00 -11.00 -0.70 -0.40 0.40 0.70 3.00 11.00 11.00 11.00 11.00 3.00 0.20 0.40 0.40 0.20 -0.7080 -1.3351 -0.0624 0.0624 1.3351 0.7080 0.0840 0.1677 -0.1677 -0.0840 14.43 14.28 8.73 -8.73 -14.28 14.43 -0.27 -0.32 0.32 0.27 3.06 9.89 14.40 14.40 9.89 3.06 0.49 0.14 0.14 0.49 -0.7507 -1.1937 -0.6406 0.6406 1.1937 0.7507 0.0298 0.3598 -0.3598 -0.0298 Table 2.3 Location of Singularities and the Associated Coefficients for Flow Past Rigid Circular Cylinder F I G U R E 2.7 Pressure Coefficient, C on a Rigid Circular Cylinder. p Two and Three Dimensional Boundary Element Analyses / 2.5.3 24 Ay B3 B4 B2 -T = 3 0 0 V T =0 dT =0 dn 2 T = 0- Bl 3.0 B6-12 13 -i— -3.0 F I G U R E 2.8 2.5.3 Discretization Scheme and Initial Singularity Locations for Rectangular Prism, • Boundary Node, x Singularity. H E A T C O N D U C T I O N IN A R E C T A N G U L A R PRISM Heat conduction in a rectangular prism is used as an example of a prescribed mixed boundary condition problem. This example has been used by Brebbia and Dominguez (1977) and Patterson and Sheik (1982) to illustrate the application of their proposed boundary element formulations. The discretization scheme, the initial locations of the singularities and the boundary conditions for the problem, are shown in Figure 2.8. The initial and final locations of the singularities and their coefficients are given in TaWe 2.4. The movement of the singularities is traced in Figure 2.9. From the plot, we see that the solution has not reached convergence as the singularities are still moving. In fact, one of them has taken off to what appears to be infinity. Inspite of the still continuing motion of the singularities, the solution process was halted. This is because there was very little change in the solutions themselves at that stage of iteration. Also, the RMS value of the boundary residue has reached the prescribed value of 10 . From all three plots of singularity motion, we see -3 the singularities exhibit a general trend of moving away from the domain of the problem Two and Three Dimensional Boundary Element Analyses / 2.5.3 25 as the solution approaches convergence. Initial Values X y Final Values Coefficient 3.20 1.00 0.0198 3.20 1.00 -1.00 -3.20 -3.20 3.20 3.20 3.20 3.20 1.00 0.0093 -0.0120 -0.0236 -0.1030 -0.0599 Table 2.4 Coefficient y 10.85 2.92 x 10 1.44 -2.99 -7.77 -8.90 8 0.0560 3.34 5.32 x 10 10.30 8.71 6.93 2.07 1.99 x 10 -0.0702 -0.0676 -0.2107 -0.1296 7 6 Location of Singularities and the Associated Coefficients for Heat Con duction in a Rectangular Prism Location Boundary Points Bl(-3.0,0.75) B2(-3.0,2.25) B3(-1.0,3.00) B4( 1.0,3.00) B5( 3.0,2.25) B6( 3.0,0.75) Interior Points Il(-2.0,0.0) I2( 0.0,0.0) I3( 2.0,0.0) Table 2.5 X Temperature Temperature Gradient Exact 2-D BEM Exact 2-D BEM 300.0 300.0 200.0 100.0 0.0 0.0 300.0001 299.9995 200.0002 100.0003 0.0003 0.0001 -50.0 -50.0 -50.0 -50.0 -50.0 -50.0 -50.0003 -49.9986 -50.0004 -49.9995 -49.9999 -50.0001 250.0 150.0 50.0 250.0000 150.0000 50.0001 -50.0 -50.0 -50.0 -50.0003 -49.9998 -50.0002 Comparison of Temperature and its Gradient for Heat Conduction in a Rectangular Prism T h e r e s u l t s f o r t e m p e r a t u r e a n d t e m p e r a t u r e g r a d i e n t c o m p u t e d a t v a r i o u s boundary a n d i n t e r n a ] p o i n t s a r e s u m m a r i z e d i n T a b i e 2.5 a l o n g w i t h t h e e x a c t s o l u t i o n s . C l e a r l y , t h e agreement between the c o m p u t e d a n d the exact results is very good. F I G U R E 2.9 Loci of Singularity Movement in the Direction of Arrows, for Heat Conduction (Domain of Problem shown Shaded). Two and Three Dimensional Boundary Element Analyses 2.5.4 F L O W PAST A RIGID / 27 2.5.4 S P H E R E As an example in 3-dimensional space, the flow of an inviscid fluid past a rigid sphere is considered. It is a direct extension of the 2-D flow past a rigid circular cylinder that was presented earlier. As before, the infinite domain of the flow is approximated by a finite domain as shown in Figure 2.10. Hence a sphere of unit radius is embedded in the center of a cube of size 16. Symmetries about the x-y and x-z planes are invoked to reduce the y F I G U R E 2.10 Problem Definition and Coordinate System for Flow Past Rigid Sphere. Two and Three Dimensional Boundary Element Analyses / 2.5.4 28 computational cost. This implies that it is sufficient to consider a quarter of the original problem. Only Neumann boundary conditions are prescribed on all the boundaries. Unlike the other examples considered previously, this problem will be solved via the discrete least square minimization technique developed in Ho-Tai, Johnston and Mathon (1979). The discretization scheme used in the problem is depicted in Figure 2.11. The number of boundary nodes used is 279, with 108 on the sphere and 171 on the cube. The singularities were distributed as follows: 10 inside the sphere and 12 outside the cube. Observe that the number of singularities used here i 3 not given by Eq. (2.10). The reason is that, this problem is solved using the discrete least square minimization technique and for this procedure, Ho-Tai, Johnston and Mathon have suggested the following relation, ""wh <- ' 2 16 For the same reason as mentioned in the circular cylinder example, we are interested only in the pressure coefficient on the sphere. The exact solution derived for 3 an infinite domain, is given by, C {9) = l-\sm 0 2 p (2.17) where 6, defined in Figure 2.10, is the direction angle with the x-axis or the colatitude. The initial and final locations of the singularities are displayed in Table 2.6. The computed pressure coefficient on the sphere, along with the exact solution are plotted in Figure 2.12. The direction angle with the z-axis or the meridional angle, is indicated by a (see Figure 2.10). Results at four a-planes are shown. Since not all the nodes are located at the a = 30° and a = 45° planes, the results at these two planes are obtained by interpolation between two neighboring nodes. The exact pressure coefficient on any a-plane is identical, as Eq. (2.17) does not depend on a. Good agreement is seen, especially at the a = 0° and a — 90° planes. 5 See p. 149 of Bird, Stewart and Lightfoot (1960) F I G U R E 2.11 Boundary Elements for Sphere and the Flow Domain. Two and Three Dimensional Boundary Element Analyses / 2.6 Initial Values X 0.50 0.42 0.36 0.28 0.19 -0.19 -0.28 -0.36 -0.42 -0.50 8.50 8.50 4.80 3.20 -3.20 -4.80 -8.50 -8.50 -4.80 -3.20 3.20 4.80 Table 2.6 2.6 y z 0.00 0.00 0.27 0.07 0.30 0.18 0.33 0.26 0.20 0.42 0.20 0.42 0.32 0.26 0.30 0.18 0.27 0.07 0.00 0.00 3.20 1.60 3.20 6.40 4.80 8.50 1.60 8.50 1.60 8.50 4.80 8.50 3.20 6.40 3.20 1.60 8.50 1.60 8.50 4.80 8.50 4.80 8.50 1.60 30 Final Values Coefficient -0.0019 0.0030 -0.0098 0.0046 -0.0031 0.0031 -0.0046 0.0098 -0.0030 0.0019 -0.3339 -0.3163 0.0284 0.0224 -0.0224 -0.0284 0.3163 0.3339 -0.0380 -0.0196 0.0196 0.0380 X 0.38 0.26 0.12 0.26 0.13 -0.13 -0.26 -0.12 -0.26 -0.38 24.35 19.12 10.83 7.55 -7.55 -10.83 -19.12 -24.35 -16.85 -6.00 6.00 16.85 y z Coefficient 0.00 0.20 0.07 0.24 0.06 0.06 0.24 0.07 0.20 0.00 6.03 6.19 16.13 4.36 4.36 16.13 6.19 6.03 16.32 17.19 17.19 16.32 0.00 0.09 0.06 0.25 0.29 0.29 0.25 0.06 0.09 0.00 4.43 15.37 15.23 18.64 18.64 15.23 15.37 4.43 5.45 3.52 3.52 5.45 - 0.0003 0.0056 - 0.5420 0.0004 - 0.0037 0.1760 -0.0004 0.5420 - 0.0056 0.0003 -29.8991 -31.7628 -17.2196 - 4.2419 4.2419 17.2196 31.7628 29.8991 21.4340 1.7581 - 1.7581 -21.4340 Location of Singularities and the Associated Coefficients for Flow Past Rigid Sphere S U M M A R Y The distinctive feature of the boundary element method presented here is that the singularities of the free-space Green's function are positioned on an auxiliary boundary which is located outside the problem domain. This boundary is allowed to move via a Galerkin minimization criterion. This results in a highly accurate, adaptive but nonlinear method. Since good algorithms for solving systems of nonlinear equations are available, the latter Two and Three Dimensional Boundary Element Analyses F I G U R E 2.12 / 2.6 Pressure Coefficient on a Rigid Sphere. 31 Two and Three Dimensional Boundary Element Analyses / is n o t seen a s a m a j o r d r a w b a c k . 2.6 32 U n l i k e conventional b o u n d a r y element m e t h o d s w h i c h require the same n u m b e r o f singularities a n d b o u n d a r y elements, t h e n u m b e r o f singularities i n t h e proposed formulation c a n be significantly reduced w i t h o u t a n y apparent loss o f a c c u r a c y . T h i s i s a n i m p o r t a n t a d v a n t a g e i n h i g h e r d i m e n s i o n a l p r o b l e m s . the From p l o t s o f t h e l o c i o f s i n g u l a r i t y m o t i o n , we see t h a t a s t h e s o l u t i o n p r o g r e s s e s t o w a r d s convergence, the singularities e x h i b i t a general t r e n d of m o v i n g away f r o m t h e d o m a i n of the p r o b l e m . W h e n t h e s o l u t i o n is v e r y n e a r c o n v e r g e n c e , w e g e n e r a l l y o b s e r v e m i n i m u m m o v e m e n t of t h e s i n g u l a r i t i e s . C o m p u t e d v a l u e s o n t h e b o u n d a r y , a s w e l l a s t h o s e i n t h e i n t e r i o r , a r e a c c u r a t e l y o b t a i n e d . I n t h e c o m p u t a t i o n o f d e r i v a t i v e s o n t h e b o u n d a r y , conv e n t i o n a l b o u n d a r y e l e m e n t t e c h n i q u e s c a n y i e l d v e r y p o o r r e s u l t s w h i c h i s n o t t h e case for the present m e t h o d . It s h o u l d b e p o i n t e d out that the p r o p o s e d m e t h o d w o u l d b e m o r e e x p e n s i v e t o r u n t h a n c o n v e n t i o n a l b o u n d a r y e l e m e n t a p p r o a c h b e c a u s e t h e m e t h o d is n o t only nonlinear b u t also requires e x t r a iterations i n c o m p u t i n g t h e o p t i m u m s i n g u l a r i t y locations. O n e w o u l d therefore have t o balance this e x t r a c o m p u t i n g cost w i t h the need to o b t a i n h i g h l y a c c u r a t e d e r i v a t i v e s o n t h e b o u n d a r y . T h e p r o p o s e d m e t h o d w i l l b e u s e d i n t h e flow a n a l y s i s f o r c o m p u t i n g t h e fluid l o a d i n g o n p n e u m a t i c m e m b r a n e s t r u c t u r e s . CHAPTER 3 Two 3.1 The Dimensional Finite Element Analysis INTRODUCTION i n t e r a c t i v e a n a l y s i s of a c y l i n d r i c a l p n e u m a t i c m e m b r a n e s u b j e c t e d t o w i n d l o a d s i s p r e s e n t e d here. S i n c e p n e u m a t i c m e m b r a n e s e x h i b i t s u c h l a r g e d e f o r m a t i o n s i n w i n d , i t is e s s e n t i a l t h a t one i n c o r p o r a t e i n t o t h e a n a l y s i s , t h e fluid-structure interactions. In this t h e s i s , t h i s w i l l b e d o n e b y c o u p l i n g t h e finite e l e m e n t p r o g r a m w h i c h p r e d i c t s t h e s t r u c t u r a l d e f o r m a t i o n s w i t h the b o u n d a r y element p r o g r a m t h a t c o m p u t e s the w i n d loads. To derive the pertinent differential equations for the s t r u c t u r a l analysis, a c y l i n d r i c a l pneum a t i c m e m b r a n e s t r u c t u r e whose cross-section is c i r c u l a r i n the u n d e f o r m e d c o n f i g u r a t i o n is c o n s i d e r e d . B y t a k i n g t h e c y l i n d e r t o be s u f f i c i e n t l y l o n g c o m p a r e d t o i t s c r o s s - s e c t i o n a l d i m e n s i o n s , t h e p r o b l e m b e c o m e s two dimensional, thereby p e r m i t t i n g some simplifica- tions. T h u s an i n i t i a l l y c i r c u l a r c y l i n d e r is d e f o r m e d i n t o a non-circular c y l i n d e r by means of a p l a n e m o t i o n . The pressure p a c y l i n d r i c a l m e m b r a n e h e l d d o w n a t t h e edges, i s s u p p o r t e d b y i n t e r n a l and l o a d e d by e x t e r n a l pressure p , w t i o n as s h o w n i n Figure 3.1. due to w i n d b l o w i n g i n the broadside direc- G r a v i t y f o r c e s are n e g l e c t e d i n c o m p a r i s o n w i t h t h e pressure loadings. A s m a l l strain-finite d e f o r m a t i o n m e m b r a n e theory is derived f o r the s t r u c t u r a l analysis. The finite r e s u l t i n g second-order integrodifferential equation is solved v i a a G a l e r k i n element scheme using N e w t o n - R a p h s o n i t e r a t i o n w i t h line-search. C o m p a r i s o n w i t h p u b l i s h e d results is m a d e whenever possible. 33 Two Dimensional Finite Element Analysis / 3.1 34 Pw{<i>) (Aerodynamic Pressure) F I G U R E 3.1 Wind Loads on Cylindrical Membrane. Some results of the non-interactive analysis of cylindrical pneumatic membranes have been presented by Ross (1969) and Uemura (1972). In Ross's formulation, the second-order integrodifterential equation is expressed in a radial coordinate defining the deformed configuration. Solutions were then obtained by numerical integration using a trial and error approach. Among some of the problems he encountered, not suprisingly, was the problem of convergence. Using a Galerkin scheme with the Fourier sine series as trial functions, Uemura converted the integrodifferential equation into algebraic equations which are then solved numerically. The integrodifferential equations here is expressed in terms of the hoop stress and the radial displacement component. This enables the membrane tension and deformation to be obtained directly from the solution of the algebraic equations. His approach is, however, approximate since the effects of a variable radius of curvature due to a deforming membrane are not captured. The results of the Two Dimensional Finite Element Analysis / 3.1 35 non-interactive study presented here, are compared with those obtained in a perturbation analysis for the case of low wind speed and by Uemura's solutions for the case of moderate to high wind speed. Also, as an independent check of this two-dimensional formulation, the cylindrical membrane will be resolved in Chapter 4 where the analysis used, is based on a different formulation. The results will be discussed in that chapter. Analysis that accounts for thefluid-structureinteractions of pneumatic membranes in wind is rare. Interaction effects are usually obtained by iterating between the wind pressure computation from an assumed wind flow model, and the structural deformation computation from the solution of the integrodifferential equation. The iterative process is stopped when the structure attains an equilibrium configuration that is consistent with the fluid loading. One such interactive analysis of wind loading on cylindrical pneumatic membranes has been performed by Kunieda, Yokoyama and Arakawa (1981). To facilitate analytical treatment of real flow around bluff bodies, they introduced simple flow models that predicted accurately, pressure loads arising from the viscous effects of flow separation and wake formation. With the pressure loads known, the integrodifferential equation for the structural analysis, expressed in radial coordinates as in Ross, is solved via a numerical technique. The iteration scheme mentioned above is then invoked until equilibrium is attained. The wind flow model used in this thesis is governed by incompressible potential flow theory. Real flow effects are not considered here as they are outside the scope of the present study. It should be pointed out that almost all air-supported stadium roofs that have been built so far, have very low rise-to-span ratios. This creates minimum disturbance to the fluid flow field. Also, in the study involving cylindrical pneumatic 1 membranes, Kunieda et. al. presented results that showed at most, a 4% difference between the deflection response predicted by the separated flow model and the potential flow model. The potential flow model, however, predicts an inward deflection at the leeward side of the membrane which is sometimes opposite to that obtained by a realflowmodel. But as discussed in Chapter 4,flowseparation still occurs in these structures, with re-attachment then taking place at the leeward side Two Dimensional Finite Element Analysis 3.2 ASSUMPTIONS AND / 3.2.1 T H E O R E T I C A L 36 F O R M U L A T I O N Figure 3.1 shows a cylindrical membrane inflated by internal pressure, p and acted upon a by external pressure, p . Despite the large structural deformations, the internal pressure w is taken to be constant. This is justified by assuming that the inside air is motionless since the structure envelops an area so large that the change of volume of air due to the finite deformation can be considered negligible. The cylinder wall is postulated to be perfectly flexible in bending, but possesses finite extensional rigidity. Small strains are assumed in the analysis and this permits the use of a linear orthotropic law in the constitutive relations. The rotations and displacements may however, be large and thus a finite deformation membrane theory will be derived. To simplify the resulting integrodifferential equations, a long cylindrical membrane is assumed. Finally, the self weight of the membrane is assumed to be small compared to the inflation pressure. 3.2.1 EQUILIBRIUM EQUATIONS Figure 3.2 depicts the shell element where the positive stress resultants are defined. The equilibrium equations, taking large rotations and displacements into account are given in Donnell (1976) as, dN x 1 dN x<l> dx R d<f> dN + 1 dN+ dx R d<j> x 0 (3.1) 0 (3.2) (3.3) where N , x N ^ are the usual membrane stress resultants. x< Two Dimensional Finite Element Analysis / 3.2.2 NI F I G U R E 3.2 3.2.2 37 R<M> jc,U Positive Stress Resultants. STRAIN-DISPLACEMENT RELATIONS For a thin c y l i n d r i c a l shell undergoing large deformations, t h e strain components e xi are g i v e n by, H[6) © 0] 2+ dx * - 5 (Ij ") [(£) ^ + = Ua? + di) + + RTxa* 2+ + RTX\M + w ( 3 4 ) - ») S ")']» <> ' 2+ ) + Rdi w h e r e u, w, tw a r e d i s p l a c e m e n t s i n t h e x ,y, z d i r e c t i o n s a n d s 0 determined later. I( - 2+ dx £ e^, 7 ^ + +£ {at - ) v -> (3 6 is a prestrain that will be 53 Two Dimensional finite Element Analysis 3.2.3 C O N S T I T U T I V E / 38 3.2.4 LAWS A s s u m i n g s m a l l strains despite the large deformations, t h e constitutive relations for a m e m b r a n e w i t h o r t h o t r o p i c properties are, N. N. NA. N, ^ = f e ( 3 - 9 ) w h e r e E, G, / i a r e t h e e l a s t i c a n d s h e a r m o d u l i , a n d P o i s s o n ' s r a t i o s r e s p e c t i v e l y a n d t i s the m e m b r a n e thickness. 3.2.4 SIMPLIFICATIONS If t h e c y l i n d e r i s s u f f i c i e n t l y l o n g c o m p a r e d t o i t s c r o s s - s e c t i o n a l d i m e n s i o n , t h e d e f o r m a tions c a n be considered t o b e two dimensional. T h e c y l i n d e r t h u s degenerates into a s t r i n g i n t h e y — z p l a n e . H e n c e e a n d a l l d e r i v a t i v e s w i t h r e s p e c t t o x c a n b e set t o z e r o . A l s o , x we c a n assume t h a t t h e r a d i a l displacement component, w is larger than the circumferen- t i a l d i s p l a c e m e n t , v. T h u s w e w i l l i g n o r e t h e v a n d i t s d e r i v a t i v e s i n a l l q u a d r a t i c t e r m s i n c o m p a r i s o n w i t h t h e to a n d i t s d e r i v a t i v e s . T h e v a l i d i t y o f t h i s a p p r o x i m a t i o n will be checked w h e n t h e m e m b r a n e is resolved using t h e three-dimensional f o r m u l a t i o n where such approximations are not introduced. T h e results o f this checking w i l l b e discussed i n C h a p t e r 4. T h e e q u i l i b r i u m e q u a t i o n s t h e r e f o r e s i m p l i f y t o , 1 dN R _ = 0 d<f> x<t> (3.10) ldN _ 4> = 0 R d<f> M R « " w F n a) 1+ + (3.11) L =RiPa ~ Pw) - (3 12) Two Dimensional Finite Element Analysis / 3.2.4 39 T h e strain-displacements relations reduce to, e = 0 x = 5S) +w + + O'] ° = 0 +e (3.13) (3ll4) (3.15) F r o m E q s . (3.10) a n d (3.11) w e have, N<t» Nx4> = c o n s t a n t s (3.16) a n d f r o m E q s . (3.9) a n d (3.15) w e c o n c l u d e t h a t N+ = 0 (3.17) x I n v i e w o f E q . (3.13) we h a v e f r o m E q . (3.7), ' N = "+{E£) N+ - (3 18) S u b s t i t u t i n g E q . (3.18) i n t o E q . (3.8) y i e l d s , e^^CN^ where _ /1 - M X M A V Eliminating (3.19) EA ) f r o m E q s . (3.14) a n d (3.19) w e have, T h e p r e s t r a i n , e i s e v a l u a t e d a t n o w i n d c o n d i t i o n (i.e. v = tw = 0). T h u s , w e h a v e f r o m 0 E q s . (3.12) a n d (3.20), e = CRp 0 E q . (3.20) t h e n b e c o m e s , a (3.21) Two Dimensional Finite Element Analysis Eliminating / 3.2.5 40 (dv/d<P) f r o m E q s . (3.12) a n d (3.22), y i e l d s t h e d i f f e r e n t i a l e q u a t i o n f o r t h e problem, ,d?w 2i2—rr+ d<f> 2 (dw — V _- 22 —— ' \ \d<f> ) - — — + w,.2 + 2Rw d<f> d<j> = f(<j>) (3.23) where /(*) = 2R> [ l + C{N+ - R} Pa - T h e c y l i n d r i c a l m e m b r a n e i s fixed a t t h e edges, r e s u l t i n g i n t h e f o l l o w i n g b o u n d a r y c o n ditions, namely = v{9 ) »(*!) 2 = «;(*!) = u»(* ) = 0 (3.24) 2 T h e s t r u c t u r a l a n a l y s i s o f a p r e s s u r e - l o a d e d s t r i n g i s h e n c e d e f i n e d b y t h e n o n l i n e a r equat i o n d e p i c t e d i n E q . (3.23), u s e d i n c o n j u n c t i o n w i t h t h e b o u n d a r y c o n d i t i o n s g i v e n i n E q . (3.24). S o l u t i o n t o t h e n o n l i n e a r e q u a t i o n w i l l b e o b t a i n e d v i a a finite e l e m e n t s c h e m e . T h e r e are however two u n k n o w n quantities i n E q . (3.23) a n d these a r e t h e h o o p s t r e s s , a n d t h e r a d i u s o f c u r v a t u r e R((f>). T h e y w i l l b e d e t e r m i n e d n e x t . 3.2.5 H O O P T h e h o o p stress, STRESS, is obtained conditions listed i n Eq. by integrating Eq. (3.22) a n d i n v o k i n g t h e b o u n d a r y (3.24) t o e l i m i n a t e t h e v. T h i s y i e l d s , J + w + {27? [ ( $ ) d<f> (3.25) C*?R{4>)d<f> *1 Observe f r o m Eq. (3.11) t h a t N$ is a c o n s t a n t a n d c a n t h u s b e t a k e n o u t o f t h e i n t e g r a t i o n s i g n as s h o w n i n E q . (3.25). Two Dimensional Finite Element Analysis i / 3.2.6 41 i Deformed Configuration F I G U R E 3.3 3.2.6 Deformation of an Initially Circular String. R A D I U S O F C U R V A T U R E , R(<f>) Consider an initially circular string of radius RQ , undergoing large deformation as shown in Figure 3.3. A point, initially on the undeformed configuration, is located on the deformed configuration by, y = (R + rv) cos <j> + v sin <j> (3.26) z — (RQ + to) sin <j> — v cos <f> (3.27) 0 The radius of curvature at any point, R(<f>) is defined as, (3.28) Two Dimensional Finite Element Analysis / 3.3 42 Observing that dz dz d(f> dy d<j> dy %=($)/ (10 - {%) (0) / (10 and Differentiating Eqs. (3.26) and (3.27) and substituting the results into Eq. (3.28) we have, \(w' + vf + {i\ + w-v'f R(<j>) = -(RQ 3/2 L + W-V') W" + (tw' + v) (2tw' + v - v") + (RO + W-V'){RQ + W- 2V') However, to be consistent with the approximations introduced in the simplified straindisplacement relations, the v motions in the above quadratic terms are neglected in comparison with the tw motions. Thus we have, 13/2 W = R , „ + tw) tw" + (2tw' ,.,„. ,.,aw) {RQ ) + (RQ.+ (3-29) 2 As before, the validity of this approximations wil be checked in Chapter 4. 3.3 GALERKIN FINITE ELEMENT FORMULATION The integrodifferential equations defined by Eqs. (3.22), (3.24) and (3.29) are solved using the finite element method. Since the profile of the membrane is unknown apriori, an iteration scheme will be used to obtain the solutions. The formulation will now be described. The string is modelled by piecewise linear finite elements. Over a typical finite element as depicted in Figure 3.4, the local displacement tw(^), is interpolated as, w(<f>) = NjWj i.e. j = 1, 2 = ^jKh ~ t) ( * - * ! ) ] { ! £ } (3-30) Substituting Eq. (3.30) into Eq. (3.23) yields the residue 3fc given by, dt = 2{RN'j - R'N'j + RNj)wj + (JV?JVj + NiN^wj - f (3.31) Two Dimensional Finite Element Analysis / 3.3 43 Typical Membrane Element. F I G U R E 3.4 Minimizing the residue over the domain of a typical finite element results in the equation of equilibrium namely, (3.32) where the linear and nonlinear stiffness matrices, and load vector are given respectively by, <h KW = 2 j [RWNj - NlN'j) - R'iNiN})] d<f> <t>i \NL) K (3.33) (3.34) = Pi = j fNid* <t>i (3.35) Two Dimensional Finite Element Analysis / 3.3 44 To solve the nonlinear equation shown in Eq. (3.32), a Newton-Raphson iteration scheme is derived by defining as follows, Ti = Pi - SijWj (3.36) where Given solutions at previous configuration ° w , the solutions at current configuration w l are obtained by, 'w = ° w + A w (3.37) Expand Eq. (3.36) by Taylor's series about ° w and neglecting quadratic and higher order terms, we get dTi dwj Awj = -Ti (°w) (3.38) 'w Differentiating Eq. (3.36) yields, i.e. dTi d awj uWj —± = - '* w - 3.39 k Defining the following quantities we have, Gii = ^ ^ h J = 1, 2 (3.40) i, j, k = 1, 2 (3.41) d K ^ Hij = d * w k Eq. (3.39) becomes = Gij ~ Hij - Sij (3.42) In view of Eqs. (3.36) and (3.42), Eq. (3.38) becomes, (Gij - Hij - Sij)^ Awj = (SijWj - P , ) k i,j = 1,2 (3.43) The hoop stress, JV^ is obtained by substituting Eq. (3.30) into Eq. (3.25). Discretizing the string by N finite elements, we have, e / 3.4 Two Dimensional Finite Element Analysis {it! E (2* 45 [WW + W)»i"j\ N + °R Pa) #} 2 + Similarly, substituting Eq. (3.30) into Eq. (3.29) yields the radius of curvature, 3.4 EVALUATION OF FINITE ELEMENT MATRICES In order to be able to integrate some of the terms derived previously, it is necessary to v assume some form of wind pressure distribution over the domain of a finite element. We will assume here that the wind pressure remains constant over the domain of a finite element, but will be different for different elements. Similarly, the radius of curvature will also be assumed constant over the domain of a finite element; this is achieved by computing the average radius of curvature. Substituting Eq(3.30) into Eq. (3.44) yields the hoop stress given by, £{*[ ^W fi n=l ^ + ^) + % a « i ^ ' + A^(tw + w) + 1 2 2CR p A<j>\ 2 a > n L 2C(* 2 n=l (3.46) {R} n Similarly, substituting Eq. (3.30) into Eq. (3.45) and computing the average radius of curvature we have, 3 { + ( * o + " 1 ) . 2 } The various stiffness matrices, I 2 3 / 2 " {(^) +(*>+^ ) } 2 2 - + (*o + "i) x 2 — — 2 ( ^ ) (, /2 2 j - • 2(=^) —-= (3.47) + {Ro+ 102)1 G,y, Hij are then computed by substituting Eq. (3.30) into Eqs. (3.33) and (3.34) and the results are integrated to yield, Two Dimensional Finite Element Analysis / 3.4 46 = 2R i K(NL) ] = T2 (3.48) (3u, 1 + to^) u; 2+ («, 1+ u , -5I^i) 2 (3.49) Since 5.-v = we have, [5] = ^12 (3.50) .^21 ^22. where {3w + tu ) l -u>2) 2 2 <S R + 12 ( u» + u>) 1 -<o ) 2 2 2 - W2) ^21 2 (l -«>*) w + 2R + 2R From Eq. (3.40), the stiffness matrix G,y is given by, and this yields 1 (3.51) <7wf TJtcJJ where 7\ = R A# <C + 2 • L dwi dw-t -R (P« ~ P J } 7 = RA<f> | + C (2JV, - 3 * ) - * f a ^ - P « ) j 3 2 2 Pa and (dN^/dwi), (dR/dwi) are defined in Appendix B. The stiffness matrix, 1 see Appendix B. is given Two Dimensional Finite Element Analysis / 3.4 47 from Eq. (3.41) as, (NL) Hi >' ~ -dw- W K But as shown in Eq. (3.49), K\^ ^ is a linear function of w . Thus we have, L k (NL) dK) ^ - W k = K^ NL) Hence [H\ = [jr<™>] = (3.52) 12 The external load vector P,, is obtained by from Eqs. (3.30) and (3.35) as, P P = R &<f> (3.53) 2 l = 2 The unbalanced load vector can thus be computed from Eqs. (3.50) and (3.53) as, (3.54) Rewriting the incremental equilibrium equations for a typical finite element in Eq. (3.43), we have, (Gij - Hij - Au, J = (P,) k i,j = 1,2 (3.55) The elemental stiffness matrices and load vectors are then assembled in the usual finite element manner, to yield their corresponding global quantities. The incremental displacements are then solved and added to the previous displacements to get the current displacements. This is indicated by Eq. (3.37). Further details of the computer implementation including the convergence criteria will be discussed next. Two Dimensional Finite Element Analysis 3.5 C O M P U T E R / 3.5.1 48 I M P L E M E N T A T I O N The governingfiniteelement equations derived for the structural analysis were implemented in program N F E A M 2 . To capture thefluid-structureinteraction effects, the finite element 3 program is coupled to the boundary element program presented in Chapter 2. This coupled program is named I N F E A M 2 . Figure 3.5 depicts the structure of I N F E A M 2 . 4 Starting with an assumed profile of the cylindrical membrane, the flow analysis is performed using the boundary element program. This enables the wind loads on the membrane to be calculated. With the loading known, the structural analysis is carried out next, using the finite element program to compute the new membrane profile. This new profile is then compared with the assumed profile. If the two profiles agree within a specified tolerance, the iterative process is stopped. Otherwise, the iteration is continued by calling upon the boundary element program to compute new wind loads corresponding to the new membrane profile. The process is repeated until convergence is achieved. The converged profile is now in equilibrium with the final wind loads. 3.5.1 C O N V E R G E N C E CRITERION By specifying that the two profiles agree within a set tolerance, we are using a displacementtype convergence criterion. Convergence could also be defined in terms of the relationship between the internal loads and the externally applied loads, the so-called unbalanced load concept. This is indicted by the r.h.s. of Eq. (3.43). When an equilibrium configuration, consistent with the fluid loading is attained, the unbalanced load term should vanish. However, Bergan and Clough (1972) have argued that a criterion based on displacements is preferable. Their arguments are as follows: A convergence criterion based on unbalanced loads does not always make sense, because the force quantities to be compared can be of completely different orders or even different dimensions. For instance, how does one decide if the moments and the out-of-plane forces in a buckled plate have converged, if Nonlinear finite Element Analysis of Membranes: 2-D * Interactive Nonlinear Finite Element Analysis of Membranes: 2-D s Two Dimensional Finite Element Analysis / 3.5.1 F I G U R E 3.5 49 Structure of INFEAM2. one compares them with the external loads comprising only the in-plane forces? A more appropriate measure of convergence is to compare the unbalanced forces with the stiffness properties of the structure. But this corresponds to working with displacements, and thus a direct use of a displacement criterion would be preferred. Two Dimensional Finite Element Analysis 3.5.2 ITERATION 50 / 3.5.2 T E C H N I Q U E It should be mentioned that a separate iterative process is required in both the boundary element program and the finite element program, for them to reach their respective solutions. The local 4 iterative process for the boundary element program has already been presented in Chapter 2. Hence we will describe only the local iteration technique for the finite element program. The Newton-Raphson iteration is used to solve the non-linear structure in the finite element program. The user has a choice of prescribing either the standard or the modified Newton-Raphson iteration technique. To prevent solutions from diverging, a line search parameter, fa is used in the solution algorithm. This idea of a line search 6 parameter is well known and has been used by several people, notably, Mattheis and Strang (1979), Crisfield (1979, 1983), and Simon and Powell (1982), to solve difficult nonlinear problems. A good description of this solution technique can be found in a recent book by Dennis and Schnabel (1983). Basically, the concept of a line search algorithm 6 involves seeking a new solution twj. +1 at the fc iteration: th in a given descent direction d as follows, k ™k+i = v>k + where fa > 0 is the line search parameter that makes tuj. The term line search refers to a technique for selecting the fad-k +1 (3.56) an acceptable next iterate. . When fa = 1, we have the standard Newton-Raphson iteration. A line search algorithm based on the suggestion of Dennis and Schnabel has been implemented in NFEAM2 and INFEAM2. 7 When the solution estimate is far from the true solution (this usually happens at the beginning stages of the iteration), the value of fa is found to range typically from 0.15 to 0.3. In this case of p < 1, a retardation effect is produced in the iterative process. On the other hand, k * To distinguish between the iterative process for the overall program to reach solution and the iterative process within either the boundary element program or the finite element program, the former will be called a global iteration while the latter will be known as a local iteration. Also known as an under-relaxation parameter when 0 < < 1, see Cook (1981) cit. loc, p. 359 cit. loc, p. 116 See Appendix D for more details 5 6 7 Two Dimensional Finite Element Analysis / 3.6.1 51 when the solution estimate is near the true solution (as usually happens towards the end of the iteration stage), it is possible to have fa > 1. 3.6 N U M E R I C A L RESULTS This section is presented in two parts: non-interactive analysis from NFEAM2 and interactive analysis from INFEAM2. Under non-interactive analysis, we will report results pertaining to the checking of the finite element program, NFEAM2. Having established the validity of NFEAM2, it is then coupled to the boundary element program to generate some numerical solutions depictingfluid-structureinteractions. The results are given in the part on interactive analysis section. 3.6.1 NON-INTERACTIVE ANALYSIS As the title implies, thefluid-structureinteraction effects are not considered. Hence the fluid loading will remain constant as the structure undergoes large deformations. The rationale for adopting this simplification is simply that we are only interested in verifying the accuracy of the proposed finite element formulation programmed in NFEAM2. The checking will be carried out in the following manner: in the case of low wind speeds, the finite element solutions are compared with results from a perturbation analysis, and for the case of moderate to high wind speeds, comparison with solutions from Uemura (1972) are used. However as mentioned earlier, in Uemura's formulation, the effects of a variable radius of curvature were ignored. This resulted in a simpler integrodifferential equation. One should therefore bear this difference in mind when making the comparison. Case a: Law Wind Speed — Perturbation Analysis The details of formulation for the perturbation analysis are given in Appendix C. Only the pertinent points will be presented here. At low wind speeds, the cylindrical membrane is not expected to exhibit large deformations. This is especially true when the internal pressurization is much larger Two Dimensional Finite Element Analysis / 3.6.1 52 than the externally applied wind loads. Under this condition, the radius of curvature does not change appreciably and can thus be assumed constant. With this approximation, the integrodifferential equation becomes analyzable by the proposed perturbation technique. Observing that the deformation of the cylindrical membrane is small in relation to its diameter, a convenient definition of the perturbation parameter e can be made in terms of the relative magnitudes of these two quantities namely, 8 where w is an unknown characteristic displacement introduced to non-dimensionalize the t integrodifferential equation in terms of a dimensionless quantity g(<f>) given by, «W) = «»* 9(<f>) (3-58) Also the wind loads acting on the membrane is assumed to be given by an arbitrary cosine law, that is, cos(o^) Pw(4>)=P. (3.59) where p» is a characteristic wind pressure and a is an arbitrary parameter. This quantity p t has to be prescribed for a given wind. At low wind speeds, p. must be necessarily small while the internal pressure in the membrane p , can be set arbitrarily large. This a would then permit the introduction of another small parameter r/, defined as the ratio of the wind pressure to the internal pressure. Thus we have, r, = ^ (3.60) Pa It is then shown in Appendix C that for terms of the same order of magnitude to balance, we must have rj = s. Hence, equating Eqs. (3.57) and (3.60) yields, w =2R[^) t 8 (3.61) Actually this definition of e is not only convenient bnt is necessary and obvious as it multiplies the nonlinear terms in the integrodifferential equation depicted in Eq. (C.9) of Appendix G Two Dimensional Finite Element Analysis / 3.6.1 53 In view of the discussion just presented, the integrodifferential equation for the structural analysis can be expressed as follows, 1 + C(JV^ - R ) 9" + e[9 + 9 ]+9 = Y ,2 2 Pa £ - Rp (l a — ecosoxp) (3.62) 4> where $2 + N c{* - = » / ^ 2 ^ l} + ^ + 2eg }* d (3.63) < + Rp It was found that both the non-dimensionless displacement g((f>) and the hoop stress must be perturbed and this is done as follows, g{4>) N = gM) + eg {4>) + e &{<!>) + 0{e ) 2 (3.64) = N + eN + e N + O{e ) 2 <j> z x 0 1 3 (3.65) 2 Substituting Eqs. (3.64) and (3.65) into Eqs. (3.62) and (3.63) result in the following sets of linear recussive relations, 9Z + g = -[cN + $] + l 0 2 1 l 1 g'l + g = I [CN + g - (&)] 2 2 x g ^ + g2 = i[cN 3 + +2 ^ + - 1 ^ cos(a^) - g' - g 2 0 a 2 0 (3.66) ( ^ y - 2 ^ - 2^1 [ ( f t ) cos(o^) ~ ft cos( <fi) - 2g gx Q and AT = RPa 0 1 C(* 2 s -^7)S $2 %d4> (3.67) Two Dimensional Finite Element Analysis / 3.6.1 54 The homogeneous boundary conditions for Eq. (3.66) are given by, <7o($i) = 0o($2) = 0 <7i($i) = 9i($2) = 0 (3.68) 02($i) = <fe($2) = 0 Solutions of Eq. (3.66) with the hoop stress and the boundary conditions defined by Eqs. (3.67) and (3.68) respectively are not too difficult to obtain. They are given in Appendix C. However, the amount of effort required to solve this set of linear integrodifferential equations increases rapidly, as higher and higher order solutions of the perturbation analysis are desired. For example, in 0 ( e ° ) solutions, only 4 terms need to be integrated. In 0(e ) solutions, the number of terms requiring integration jumps to 20. For 1 this reason, the solutions of the perturbation analysis will be obtained to at most 0(e ) 2 accuracy. The proposed perturbation technique is then applied on a numerical example. The finite element solutions for this example are also generated. In this way the results of the two analyses can be compared. Consider the problem depicted in Figure 3.6. The membrane structure is subjected to a cosine wind load. Because the wind load distribution is asymmetric, a more realistic wind pressure is simulated. From the given data, the following values are obtained for the computation of hoop stress and membrane deflections using the perturbation equations': a = 1.2, p« = 3, £ = 0.10 and u>, = 2. Some comments on the units adopted in the example. Since all the examples given here are derived in one way or another from Uemura (1972), it is convenient to retain the metric units used by him. This problem is next discretized with 12 elements and solved via NFEAM2 for membrane stresses and deflections. Five iterations were required to achieve convergence. Tabie 3.1 shows the hoop stress predictions by both solution methods. Excellent agreement Two Dimensional Finite Element Analysis / 3.6.1 55 = 30° * = 150° 2 C = | x IO" 4 m/kg Membrane Roof Subjected to Low Wind Pressure. F I G U R E 3.6 is observed. Since the nonlinearities are very small for this problem, the hoop stress due to wind is slightly above the initial stress due to the internal pressure only. Perturbation Analysis o{* ) 0( >) Analysis 307.74 303.62 304.34 1 300.00 Table 3.1 Finite Element S Accuracy of Hoop Stress Computations The membrane deflections computed by these two methods are plotted in Figure 3.7. Unlike the hoop stress results which are necessarily computed to a higher Two Dimensional Finite Element Analysis F I G U R E 3.7 / 3.6.1 56 Comparison of Membrane Deflections. order of accuracy in the perturbation analysis, the membrane deflections here are obtained to 0 ( e ° ) accuracy only. The finite element solutions agree closely with the perturbation results. A note on the boundary conditions used in the perturbation analysis is worth mentioning. All homogeneous solutions in the proposed perturbation method have in their denominators, the term sin($2 (eg. $ x = 0°, $ = 2 — $i)- 9 This term vanishes at some combinations of $'s 180°). This in turn causes the homogeneous solutions to blow up. Of course one could simply prevent this by avoiding those combinations of $'s that make the denominator vanish. On the other hand, if one has to work with those combinations see Eqs. (G.27) and (G.34) of Appendix G Two Dimensional Finite Element Analysis / 3.6.1 57 of $'s then one should realize that one cannot prescribe as boundary conditions only the kinematic boundary conditions. One has to include the natural boundary conditions in lieu of some kinematic boundary conditions. Take for example, at <f> = one can prescribe the kinematic boundary condition of zero displacement. But at <f> = $ > 2 o n e n a s *° prescribe the natural boundary condition of zero slope. This situation is analogous to the problem encountered in the membrane theory for cylindrical shells, where one cannot satisfy all the boundary conditions. This problem does not arise in the finite element analysis, as we are solving the nonlinear differential equations. Having seen that the solutions from the finite element analysis are accurate for the case of low wind speeds, the program NFEAM2 will now be checked for the case of moderate to high wind speeds. Case b: Moderate to High Wind Speed — Uemura's Solutions Uemura (1972) presented some results for the non-interactive analysis of cylindrical membranes in wind. The wind speeds used in his examples ranged from moderate to high (namely from 20m/s to 40m/s). The nonlinear terms shown in Eq. (3.22) are absent in Uemura's formulation. These terms are due to the effects of the variable radius of curvature, neglected in Uemura's theory. Initially, in the undeformed state, the membrane is assumed to have a constant radius of curvature. As the membrane begins to deform, the radius of curvature starts to vary. This variation in the radius of curvature increases as the membrane deformation gets larger. One would thus expect the nonlinear terms to become increasingly important as the membrane experiences bigger and bigger loads. To obtain a quantitative measure of the loading on the membrane, one has to consider the internal pressure p , in addition to the wind pressure p . a w Since the wind pressure can be computed from a given wind speed U, the following dimensionless quantity which represents the ratio of dynamic pressure to the internal pressure, is introduced to quantify the strength of loading on the membrane, ±o U 2 Two Dimensional Finite Element Analysis / 58 3.6.1 where p is the density of air. a Another quantity that will determine the deformation response of the membrane is the elastic constant, C. Since C is inversely proportional to the Young's modulus of the membrane material, a larger C would therefore result in a larger membrane deformation at a fixed f . 0 The following two cylindrical membranes, designated as Examples 1 and 2 and shown in Figures 3.8 and 3.9 are taken from Uemura. The stresses and deformations computed by Uemura for these two membranes will form the basis of comparison for similar results predicted by NFEAM2. Both membranes, discretized with 15 elements, have initial radius of R = 10m. Comparison of hoop stress predictions is tabulated in Table 3.2. Hoop Stress, in Kg/m Membrane Uemura (1972) 2-D F E M Example 1 818.7 864.1 Example 2 1227.3 1271.9 Table 3.2 Comparison of Hoop Stress Predictions The loading ratio in Example 1 is, £ a = 1.91. This figure represents a moderate wind. For Example 2, we have £, = 3.40 and this is in the region of a strong wind. Also Example 1 is more 'rigid' than Example 2, as indicated by the smaller value of C. Hence we expect Example 2 to exhibit a larger deflection than Example 1. Their deformed profiles, together with Uemura's solutions, are plotted in Figures 3.8 and 3.9. The agreement with Uemura's result for Example I i 3 close but for Example 2, there is a visible difference. This difference is due to the nonlinear terms which are retained in the present more exact analysis. In Example 1, the nonlinearities are small and hence the Two Dimensional Finite Element Analysis / 3.6.1 <~~~~> 2 - D F E M U E M U R A (1972) F I G U R E 3.8 Membrane Deformations for Example 1 fa = 1.91, C = \ x 10~ o o o o o 2 - D FEM UEMURA (1972) F I G U R E 3.9 Membrane Deformations for Example 2 fa = 3.40, C = 1 x 10~ Two Dimensional Finite Element Analysis / 3.6.2 60 agreement with Uemura is good. However in Example 2, the deformations are larger and the nonlinearities become significant. Thus the visible difference. The study of the non-interactive analysis presented so far indicates that the results obtained by NFEAM2 agree well with the solutions of perturbation analysis and Uemura's theory. Thus the membrane stress and deformations computed at different wind speeds are accurately predicted. 3.6.2 I N T E R A C T I V E ANALYSIS Having established the accuracy of the finite element formulation on independent structure problems and the accuracy of the boundary element formulation on independent flow problems, these two formulations are coupled together to make them interactive. Hence, the examples presented in this section will exhibit the fluid-structure interaction effects. Since the membrane profile is not known apriori, the fluid loading computed in the flow analysis, is based on an assumed profile. The assumed profile used here is deliberately chosen to be slightly asymmetric, resulting in pressure loads acting on the membrane that are also asymmetric. It is therefore not suprising that the final converged profile obtained in the interactive analysis also exhibits asymmetry. Had a symmetric profile been chosen, the final converged profile would also be symmetric. Such solutions while mathematically correct, are not realistic as one would expect the membrane to be 'lop-sided' when subjected to wind blowing transversely. Another interesting mathematical solution is obtained when the assumed profile is taken to be asymmetric, but in a direction opposing to the flow. The final profile, also asymmetric, is found to converge in that direction which is opposite to the flow direction. Again this mathematical solution is not realistic. The reason for this behavior is that, the pressure loads acting on the membrane are given by the square of the flow velocities. Hence the membrane deformation is indifferent to the flow direction. To select the correct solution from all these mathematical solutions, we must also check that the solution is realistic for the given flow environment. Two Dimensional Finite Element Analysis / 3.6.2 dn 61 ii =0 dn y i\ 1 ^-Flexible Cylinder Rigid Cylinder / \y / / 1 / 1 / >\ \ \ / \) \\ x — F I G U R E 3.10 I*- If \f / • L Flow Domain and Boundary Conditions for Cylindrical mem In all the examples considered, the infinite potential flow region is modeled by a finite flow domain as depicted in Figure 3.10. Note that only Neumann boundary conditions are prescribed for the problem. This implies that the computed potential function is accurate to within an arbitrary constant. It is not necessary to evaluate this constant as the quantities we are interested in, namely, the pressures and velocities, are given by the derivatives of the potential function. The following values are used in all the interactive examples: L = 100 m, H = 100 m, p = 1.25 k g / m , 3 a R = 10 m Q (3.70) The number of boundary elements, singularities and finite elements used in the study are summarized in Table 3.3. They are evenly distributed over each section of the boundary. Two Dimensional Finite Element Analysis / 3.6.2 62 The deformations at each global iteration for two membranes are shown in Figures 3.11 and 3.12. Observe that there is a slight inward deflection at the leeward side of the membrane in Figure 3.12 which is absent in Figure 3.11. This can be explained by the larger used in the latter. Boundary Section No. of Boundary Elements No. of Singularities No. of Finite Elements Vertical 30 x 2 3x2 Horizontal 48 4 - Membrane 50 6 49 Total 158 16 49 Table 3.3 The Discretization Schemes converged pressure coefficients on the membrane are plotted in Fig- ure 3.13. The ratios of the pressure coefficient at the top of the membrane to that at the base are -3.08 and -3.56 for £ = 0.276 and £„ = 1.875 respectively. In contrast with a rigid a circular cylinder which has a pressure ratio of -3.0, we see that a flexible cylinder has a higher ratio. This arises from the fact that a flexible cylinder deforms into a shape that is more distorted than a circle. Another important point that is observed in the plot is that the wind loads based on the rigid model are underestimated at the apex of the cylinder compared to those from the flexible model. Hence, it is important to consider the effects of fluid-structure interactions in a pneumatic membrane structure subjected to wind loads. Figure 3.14 depicts the influence of the loading ratio f , on membrane defleca tions. The results are similar to Kunieda et al. (1981) who have adopted a realflowmodel in their formulation. The small inward deflection in the leeward side of the membrane is also evident in their plot. In fact, the ratio of the two inward deflections 10 The inward deflection at the windward side and the inward deflection at the leeward side in the present Two Dimensional Finite Element Analysis / 63 3.6.2 F I G U R E 3.11 Membrane Deformations at Each Global Iteration (£ = 0.276,C = | x FIGURE Membrane Deformations at Each Global Iteration (£ 3.12 a a = 1CT J. 4 1.875,C = i x i o - ; . 4 Two Dimensional Finite Element Analysis / 3.6.2 64 o F I G U R E 3.13 Pressure Loads acting on Membrane in Equilibrium ( C = | x 10~ J. 4 analysis agrees well with their results. Also, just as in their plot, this inward deflection at the leeward side vanishes at smaller values of the loading ratio. It is thus seen that the deformed shape of the membrane structure changes with £ . 0 The variation of hoop stress , with internal pressure p at different wind a speeds, is sketched in Figure 3.15. The speeds considered are 0, 10, 20, 30, and 40 m/s. The wind pressures shown in the figure are dynamic pressures, corresponding to the wind speeds mentioned. Observe that the additional stress due to wind, defined as the difference Two Dimensional Finite Element Analysis F I G U R E 3.14 / 3.6.2 65 Influence of the Loading Ratio £ , on Membrane Deflections (C = a between the total stress and the stress due only to the internal pressure, changes rapidly at small p . It then decreases with increasing internal pressure. This is attributed to the a fact that the membrane becomes more and more 'rigid' at higher p . Also, recalling that a the perturbation method is applicable to the case of small wind loads and large internal pressure, it will be interesting to compare these interactive results with the solutions of the perturbation analysis in the region of high p . A difficulty that quickly arises in the attempt a to apply the perturbation technique here, is the characterization of the wind pressure on the membrane. Naturally the characteristic wind pressure, p« required in the perturbation analysis, is obtained from the final pressure loads acting on the membrane. From Two Dimensional Finite Element Analysis / 3.6.2 66 Py, = 102 kfi/m p = 57 kg/m 2 P„ = 25kg/m p = 6 kg/m* P = 0 kg/m 8 ! w to w A ! w < ft O o 200 Internal Pressure p , kg/m 250 2 a F I G U R E 3.15 Variation of Membrane Stress Resultant P a (c=\x , with Internal Pressure io~ ;. 4 Figure 3.13, it is found that the final pressure loads on the membrane can be accurately modeled by the following equation, Pw(<t>) =p*(coso^ + a,) (3.70) where p., a and a* are constants chosen to match the final pressure distribution as closely as possible. In the examples given here, we used a = 2.0 and a„ = —0.50. The values of p, range from 12.8 to 200.0 depending on the wind speed. Using this wind pressure distribution, the 0(e ) term for the hoop stress, N 1 l is derived from Eq. (3.67) to be, aZi {1 - cos (A$)} - Z sin (A$) + (a - l) a a , Z 2 2 AT, = 1 (a - I) a {(7A$sin(A$) + (C + ^)[2 - A$sin(A$) - 2cos(A$)]} 2 where Z = cos(a$!) + cos(a$ ) x 3 2 (3.71) Two Dimensional Finite Element Analysis Z = sin (a$ ) 2 — s 2 m / 3.7 67 ( $i) a Z = A$sin(A$) - 2{1 -cos(A$)} 3 The results of the perturbation analysis, accurate to 0(e ), are shown by broken lines in l Figure 3.15. The plot shows that in the region of high internal pressure, the interactive results approach the perturbation solutions. However at higher wind speeds, the two sets of curves differ by up to 15%. This is expected and is attributed to the deterioration of the perturbation solutions. The internal pressure at this wind speed is not high enough and consequently, the perturbation parameter e, denned as the ratio of the wind load to the internal pressure, is almost 1.0 in this region. 3.7 SUMMARY A two-dimensional interactive analysis of cylindrical pneumatic membranes subject to wind loads is presented. Theflow,modeled with an inviscid fluid, is analysed using the boundary element technique. The structural analysis is based on a large deformation but smallstrain model, performed using the finite element method. Two computer programs are developed for the analysis: NFEAM2 for the non-interactive analysis and INFEAM2 for the interactive analysis. To assess the performance of NFEAM2, the following procedure is followed. At low wind speeds, the results of NFEAM2 are compared with the predictions from a perturbation analysis. At moderate to high wind speeds, these results are compared with Uemura's solutions. These comparisons show that accurate membrane stresses and displacements are obtained. Having etablished the accuracy of NFEAM2, the boundary element program is coupled to the finite element program to form INFEAM2. Results displaying fluid- structure interactions are then generated. The investigation shows that the deformed shape of the membrane changes with f . Also the additional hoop stress due to wind a changes rapidly at smaller p , and then decreases with increasing p . a a CHAPTER 4 Three Dimensional Finite Element Analysis 4.1 I N T R O D U C T I O N The interactive analysis of a wind loaded, three-dimensional pneumatic membrane structure is presented here. To study the motion of a body undergoing large deformations, a consistent formulation based on continuum mechanics principles is desired. In the classical approach for finite elasticity analysis, either a Lagrangian or an Eulerian description is used to characterize the deformation. However, in this presentation, a material coordinate description which retains some of the basic features of the Lagrangian description is used. The coordinate system employed here, is embedded in the body. As the body deforms, the coordinate system displaces and deforms with the body such that the coordinates of a point on the body remain unchanged. Hence this material description can be called a convected coordinate description. The incremental equations of equilibrium are derived by applying the principle of virtual work to two neighboring configurations. By considering the configurations to be infinitesimally close, these incremental equilibrium equations are then linearized to permit easy solution. The error arising from the linearization constitutes a measure of the non-satisfaction of equilibrium for the configuration obtained. To evaluate and control this error, an unbalanced load term is derived. The iterative process in the solution algorithm is halted when this unbalanced load satisfies a prescribed tolerance. 68 Three Dimensional Finite Element Analysis / 4.2 69 The nonlinear deformation is analysed using a displacement-based finite element procedure. The incremental continuum mechanics equations are used to develop the governing finite element equations. To evaluate the integral associated with the external load term, a truncated Taylor series expansion of the external virtual work is used. However, due to the follower-load nature of the pressure loading, we have an unsymmetric stiffness contribution to the system tangent stiffness matrix. A three dimensionalfiniteelement program has been developed for the structural analysis. The program, called NFEAM3 uses the warped isoparametric quadrilat1 eral element of Haug and Powell (1972a). This membrane element is warped in space and has 12 degrees of freedom with the nodes connected by straight lines. Newton-Raphson iteration technique is employed in the solution process. To promote convergence, an automatic line search algorithm is incorporated into the Newton-Raphson iteration scheme. A displacement convergence criterion is adopted in the program. By coupling the finite element program for the structural analysis with the boundary element program outlined in Chapter 2 for the flow analysis, the fluid-structure interactions are achieved. The coupled program named INFEAM3, is then used to analyse the air-supported roof of the BC 2 Place Stadium under the action of wind loads. 4.2 M A T H E M A T I C A L PRELIMINARIES Only the relevant mathematical relationships are given here. For a more complete presentation, one could refer to McConnell (1957) and Green and Zerna (1968). Assumptions similar to those stated in the two-dimensional formulation given in Chapter 3 are adopted here. The notations used will be as follows: Greek indices are reserved for surface quantities and take the values 1, 2 while all lower-case Latin indices take values 1, 2, 3. The upper-case Latin indices take values 1, 2, 3, • • •, 12 and the upper-case calligraphic indices Nonlinear Finite Element Analysis of Membranes: 3-D Interactive Nonlinear Finite Element Analysis of Membranes: 3-D Three Dimensional Finite Element Analysis / 4.2 70 a F I G U R E 4.1 2 Coordinate Systems for a Surface in Three Dimensional Space. take the values 1, 2, 3, 4. Some examples depicting the use of these indices are shown in §1.4 Nomenclature of Chapter 1. The summation convention for repeated indices applies unless indicated otherwise. Consider a surface in three dimensional space, as shown in Figure 4.1. Let the surface be denned by general curvilinear coordinates, rectangular cartesian coordinates, x % 3 6 ) with respect to a fixed 2 as follows, * = x '(0 ) < , a (4.1) Denoting the position vector of any point P on the surface by r, with respect to the fixed Since the coordinates x' are in rectangular cartesian system, the location of the index is immaterial. Hence x' or x,- wonld be permissible. Bnt for 0' it is meaningful to write 0' rather than 0;. Three Dimensional Finite Element Analysis / 4.2 71 origin O, we have, r = xg (4.2) i i where g,- are the unit base vectors of the fixed rectangular cartesian coordinates. The covariant base vectors of the surface at P are given by, a ° ~ d9 a8 i = (4.3) where a* = (dz'/<90 ) yields the cartesian components of the base vectors. The contrao a variant surface base vectors, a are given by, a a .ap = 8$ (A A) a where the Kronecker delta, Sjj is defined as, c 5 P = i0 forct^P; ( \ 1 for a = /?, P not summed. 4 . \ ' ) The covariant metric tensor on the surface is given by, (4.6) °a/? = a » / ? a a This yields, «3 fl = a « J 9ij = /?a a a (4.7) Similarly the contravariant metric tensor on the surface can be computed from, a a o A A / ? = $| The normal vector to the surface at P, a (4.8) 3 is given by, a = ai x a 3 2 Substituting Eq. (4.3) into Eq. (4.9) yields, « 3 = A Si = <»i<*2& x 8j (4.9) Three Dimensional Finite Element Analysis / 4.2 72 where a denotes the cartesian components of the normal vector. 3 4® i.e. = <*1 «2 ijk8 e 4 = « i °2 e k = o'i «2 ijkV99 gi e i}k\/99 kt (4.10) H where e.-yj. is the permutation tensor defined as, ijk = >/J7iyjfc e e The permutation symbol e,^ is defined as, 0 when any two of the indices are equal; when t, k is an even permutation of 1, 2, 3; when t, j, A: is an odd permutation of 1, 2, 3. (4.11) Noting that g,- are the unit base vectors in a rectangular cartesian coordinate we have, upon expanding Eq. (4.10), ,2„3 2„3 a 0>\0>22"1 ~ "> a (4.12) a\a\ - a\a\ ,1„2 a,1„2 a - a a x 2 2 1 Dennoting the magnitude of the normal vector by a we get from Eq. (4.9), 3 a (4.13) 3 = l«3l = >J4 4 9ij If e represents the unit normal vector to the surface at P, we have, 3 a 3 a 3 a xa x 2 (4.14) |aj x a | 2 where e' are the cartesian components of the unit normal vector and are given by, 3 4 > V. = ,2 — ° c / 3 3 1 a t „2„3 - a <»3 3 ,2„3 - .1„3 l ° 2 ~~ l 2 a 11 a (4.15) 11 "3 The magnitude of the surface base vector for the covariant components is given by, = (a .a ) / = 1 a a 2 v ^ (no sum) (4.16) Three Dimensional Finite Element Analysis / 4.2 and the contravariant components are given by, \a \ = v/a™ (no sum) a The direction cosine between the surface base vectors is given by, a .a 1 2 = 182)003(8!, a ) 2 or cos(a 12 \/ ii 22 a ) u a 2 a 0 From trigonometry we have, sin(a a ) = yjl - c o s ^ , a ) 2 lf 2 2 • / ^ / n ° 2 2 " 12°12 8111(8!, a ) = W — V <*nL 22 a a 2 a or sm(a u a) = J V n 22 2 a a where a = |°a/3| = °ll 22 - 12 12 a a a The line element of the surface at P, ds is given by, (ds) = dr *dr 2 (ds) = a d9 .a d9^ 2 i.e. a a j3 (is) = a p dO dB a 2 p Q The area of a surface element is found from, dA = \dr x dr \ x 2 dA = \d9 a, x dB a | l 2 2 i.e. <M = dfl <f0 \a \ \a \sin(a , a ) 1 2 x In view of Eqs. (4.16) and (4.19), dA is thus given by, dA = ^d9 d6 l 2 2 1 2 Three Dimensional Finite Element Analysis 4.3 T H E O R E T I C A L / 4.3 74 F O R M U L A T I O N The motion of a surface in a three dimensional space may be described by considering its path of deformation from an undeformed configurations, ° C to a current deformed configuration, C and a neighboring configuration, C. This is shown in Figure 4.2. It is 1 2 useful at this point to define the nomenclature adopted here. All left superscripts indicate the configuration in which a quantity is referred to. Hence the left superscript 0, refers to the original or undeformed configuration while 1,2 pertain to the deformed configurations. All right super-subscripts are indices as used in tensor analysis unless indicated otherwise. All incremental quantities are depicted without left superscripts. / 4.3.1 Three Dimensional Finite Element Analysis 4.3.1 C O N V E C T E D STRAIN 75 T E N S O R A symmetric tensor in a deformed configuration, C can be defined as, V "lap = i -V (4-23) ) where "a j, ° a ^ are the metric tensors in the deformed configuration, C and the undeV a/ a formed configuration, °C respectively. The symmetric tensor, "7 a/3 will be called a con- vected strain tensor because it measures the difference of the squares of the corresponding line elements in a convected coordinate system. That is, (d"S) -(d S) 2 0 = 2''< d9 d9P 2 (4.24) a 1af} A symmetric, incremental convected strain tensor, 7 3 can then be defined from Eq. (4.23) a/ as, la? i.e. = lafi l = ±( a - a ) (4.26) l 2 7 o j 9 (4.25) ~ 1a(i 2 afi a / ) Expanding the covariant metric tensor, a p in Taylor series about the deformed configura2 a tion, C we have, 1 -v * - '*•) + \VM (v - '*•') (V - V) ... + + Neglecting cubic and higher order terms we have, o Q/J - a a/l « - ^ p - x + 2 a l i , a l x i x x> (4.27) where x* = ( x* —x ) represents an incremental displacement term. In view of Eq. 2 1 t (4.27), Eq. (4.26) can be recast in the following form, (4.28) 7o/3 = Cap + Vafi where c a(3 is the linear incremental strain tensor given by, (<•»> and r) is the nonlinear incremental strain tensor given by, a/3 _ ! d (a ) 2 l afi ^ ^ W ? 1 • ^ • ( 4 - 3 0 ) Three Dimensional Finite Element Analysis 4.3.2 / C O N V E C T E D STRESS 4.3.2 76 T E N S O R Since the convected coordinate description is chosen for the finite deformation analysis, a convected stress tensor will be derived. Just as in the Lagrangian description of stress tensors, this convected stress tensor will be defined with respect to the undeformed configuration ° C and will henceforth be called the convected Piola-Kirchoff (P-K) stress tensor. Because stresses occur naturally in the deformed configuration, this convected P-K stress tensor will not yield directly the physical stresses. A transformation will be necessary to obtain the physical stresses. Following the developments given in Green and Zerna (1968), consider the two dimensional equivalent of the Cauchy's tetrahedron in the current deformed configuration C as shown in Figure 4.3. V F I G U R E 4.3 Cauchy's Triangle. Three Dimensional Finite Element Analysis / 4.3.2 77 This two dimensional equivalent of the Cauchy's wedge, will suffice for the description of a membrane continuum. The length of the sides is d S along the curvilinear coordinates, 6 v a a and simply d S along the third side of the triangle. The orientation of the sides are defined u by unit outward normal vectors " v and "v, as shown in Figure 4.3. When the triangle a shrinks to zero at P, keeping "v fixed, the stress vector at P is realized. The length of the third side is given vectorially by, dS = vd S u v (4.31) v or in terms of the contravariant base vectors, " a 2 vd S v = Y v a we have, va a -7=d S (4.32) v Q a=l where d S v a represents the length of the sides along the 6 coordinate lines given by, a dS v = y/^d9 (no sum) a a (4.33) If v denotes the covariant components of "v with respect to the base vectors a , we u a a have, v v V^d S (4.34) = dS v v a a For convenience, the following covariant quantity is defined, d S„ v (4.35) Eq. (4.34) can thus be written as, v d"S v a = d S* (4.36) ca (4.37) v a Also from Eq. (4.8) we have, v„a\ v„ which upon expansion and inverting yields, 1 "a " a ' "o " a ~ "a 11 12 21 where a — an a 22 — 22 O12 a 2 i "°22 -^21 (4.38) Three Dimensional Finite Element Analysis Let t and "t v a / 4.3.2 78 be the traction vectors acting along the sides of the infinitesimal triangle as shown in Figure (4.3). The equilibrium equation for the triangle is then given by, v td S v = %d S (4.39) v a Body forces acting on the triangle are ignored as they are assumed small compared to the surface forces. From Eqs. (4.35) and (4.39) we find that, 2 "t=^"t; v V''t (4.40) / a a Note in Eq. (4.40), t is an invariant to coordinate transformation for any arbitrarily v chosen "v and "v is a covariant vector. Thus it follows that \/ a v a aa v t a is a vector with a contravariant character, a. We may therefore define the following contravariant quantity, "t * = v^a ""ta a (no sum) 5 (4.41) The above contravariant quantity, t * can be decomposed with respect to either a contrav a variant or covariant basis at P as follows, or "t * = "E% V a v (4.43) E ^ and E^ are tensors by Quotient rule and are called convected-contravariant Euler a V and convected-mixed Euler stress tensors. Since we intend to work with stress tensors associated with the undeformed configuration °C, we will derive an alternative definition of stress tensors. A differential force vector, d "T can be defined from Eq. (4.39) as, dT v = td S v v (4.44) Three Dimensional Finite Element Analysis / 4.3.2 79 If "n denotes the traction vector acting on the deformed element per unit of the undeformed length d°S, we have, d T = td S = nd°S v v v (4.45) u With the aid of Eqs. (4.36), (4.39) and (4.41), we can rewrite Eq. (4.45) as, d T = t * d "S* = n * d °S* V v a v (4.46) a From Eq. (4.46), four sets of stress tensors can be defined by decomposing with respect to the covariant and contravariant basis of the deformed and undeformed configurations. We thus have, v t* = E a v afi "ap = E$ "a* (4.47) a? \ (4.48) V or "t * = F v a = F$ V V and n = n n = 9 a^ = K (4.49J or a^ = a (4.50) M Of the many possible definitions of the stress tensors, two will be of special interest, namely the convected P-K stress tensor, n ^ and the convected Euler stress tensor, E & . The u a v a relationship between these stress tensors will now be derived. From Eqs. (4.38) and (4.35) we have, Substituting Eqs. (4.47) and (4.49) into Eq. (4.46) yields, h ^ apd b = a n ^ apd b a from which "E aP But from Eq. (4.51) we have, d "S* = "n a/? d°S* (4.52) Three Dimensional Finite Element Analysis / 4.3.2 80 The above Eq. (4.53) depicts the relationship between the two stress tensors. From moment equilibrium, the convected Euler stress tensor is symmetric. This implies from Eq. (4.53) that the convected P-K stress tensor is also symmetric. The next point of interest is to establish the physical components of these stress tensors. From Eqs. (4.42), (4.45) and (4.46) we have, t d S = "t * d "5* = E & "ap d "S* v V a v (4.54) a Expressing the base vectors "a^ in terms of the unit vectors ep we have v V % H = —J= v (no sum) (4.55) HP Manipulating with Eqs. (4.35), (4.38) and (4.55) we have, a d Sl^^a%d S v v v p a In view of above, Eq. (4.54) becomes, td"S = "E \fia v d "S e td S u v a =E dS e v ap l, (4.56) v al} 0 (4.57) v a 0 in which E > = \fo E v af v (4.58) afi Since E°P is associated with the physical quantities d S and "e^ as shown in Eq. (4.57), V v a v EP A must then yield the physical components of the convected Euler stress tensor given in Eq. (4.58). In a similar manner we can obtain the physical components of the convected P-K stress tensor as follows. From Eqs. (4.45), (4.46), (4.49) and (4.55) we have, 2 td S v u = nd°S = YJ V*HP" P ° a v E d (- ) S 4 59 P=i Substituting Eqs. (4.35) and (4.51) into Eq. (4.59) we have, "nd°S = "n d °5 = ^ where n v afi n ^d S e v afi v v a d "S e =V ^ V * v a ? ? (4.60) (4.61) Three Dimensional Finite Element Analysis / 4.3.3 81 represents the physical components of the convected P-K stress tensor. 4.3.3 I N C R E M E N T A L C O N S T I T U T I V E RELATIONS For any configuration C, the convected P-K stress tensor and the convected strain tensor U are related as follows, where C ^ v a pX is the elastic tensor at configuration "C. For isotropic material, the elastic tensor, as given in Green and Zerna (1968) is, 4 (4.63) 2(1 + /*) where E, \i are the elastic modulus and Poisson's ratio respectively, t is the membrane thickness and a $ is the metric tensor in configuration C. The incremental convected v a V P-K stress tensor can be defined as, n aP 2 a/3pX 2 px + ( C ? 2 l p 1 ^ _ l^pX G = C°^ i.e The elastic tensor C 2 a0pX 2^ = X a - C p x l a f i p x ) l l p X (4.65) is an unknown apriori. However for small strain problems, we can assume negligible change in the elastic tensor in the increment. i.e. 2 ap \ \ app\ C P w C ^ 46 6 j Hence Eq. (4.65) becomes, Further simplification of Eq. (4.67) is necessary as i pX has a nonlinear component. This is done as follows: substituting from Eq. (4.28) we have, 4 cit. i o c , p. 389 Three Dimensional Finite Element Analysis / 4.3.4 82 The nonlinear incremental strain tensor, r\ \ is negligible within an increment and hence p dropped. Thus we have, The above incremental constitutive relation will be applied to linearize the incremental virtual work equations which will be derived next. 4.3.4 I N C R E M E N T A L EQUATIONS O F EQUILIBRIUM A consistent and rigorous formulation of the incremental equilibrium equations based on A denote the area of the deformed virtual work approach is presented as follows. Let V surface "S loaded by a convected Euler stress field, E P. v The principle of virtual dis- a placement for a virtual field of convected strains, 6 "t p requires that, v a jj "E"? 6 (-7^) dA = »W (4.70) ext "A where W v exl is the external virtual work in configuration C. Since the deformed surface V is not known apriori, the above integral may be evaluated only approximately by assuming a trial configuration and iterating. To avoid this kind of approximation, the integral will be transformed into an integral with respect to the undeformed configuration, ° C . The convected Euler stress field will necessarily be replaced by the convected P-K stress field. Substituting Eqs. (4.22) and (4.53) into Eq. (4.70) we have, IJ ^ 5 ( ^ ) d ° A = ^W (4.71) ext °A The incremental virtual work equations are then derived by applying the above virtual work expression to two neighboring configurations, C and C . Linearization of these 1 2 equations results when the configurations are assumed to be infinitesimally close. In configuration C , Eq. (4.71) becomes, 2 / / V^( 7^)d°A = Vext 2 °A (4.72) Three Dimensional Finite Element Analysis / 4.3.4 83 Substituting Eqs. (4.25) and (4.64) into Eq. (4.72) yields, °A That is, 6( ) (V * + Ij d°A = *W lafi (4.73) exl °A II 11 /j In view of Eq. (4.28), Eq. (4.73) becomes, ( n^ l + n^)6(e , °A n"? S^d°A °A + r ,)d A Q 11 JJ Q }a + °A = W 2 ext V » 6 p d°A Va ext °A Using Eq. (4.67) to eliminate n P in Eq. (4.74) we have, C^ l a x 7 a „A 6lafid°A °A n* 6r, p d A = + l °A Q a -Jj JJ = *W - °A V 6 which represent the nonlinear equations for the incremental analysis. The nonlinear term in the above equation is linearized as follows, 7 A hap = (e A + V \) 6{e p + V p) « e A Se P P p a a p (4.76) afi In view of the above, Eq. (4.75) becomes, JJ C^ e 6e d A l x o pX Q0 °A + JJ n^6r, d°A l afl °A elastic stiffness geometric stiffness = *W ext - JJ^6e d A 0 ap °A unbalanced load From Eq. (4.77), the linearized incremental equations of equilibrium are obtained by observing that the equation holds for arbitrary variations of compatible strains. These linearized incremental virtual work equations are of the same form presented by Bathe, Ramm and Wilson (1975) who formulated similar equations but used a Lagrangian description in a rectangular cartesian coordinate system. Observe that the incremental decomposition of the stresses and strains are possible only because they are referred to the same configuration. Solution to these equations is achieved via a finite element modelling. This will be presented next. (4.77) Three Dimensional Finite Element Analysis 4.4 ISOPARAMETRIC / 4.4 84 FINITE E L E M E N T SOLUTION A displacement-basedfiniteelement solution is used as the basis of the numerical modelling. The incremental continuum mechanics equations presented earlier are used to develop the governing finite element equations. Following the work of Haug and Powell (1972a), a warped isoparametric quadrilateral shell (WIQS) element in configuration C is chosen for V the analysis. The element has straight boundaries and 12 degrees of freedom characterized by nodal displacements, TJjy as shown in Figure 4.4. A point P on the middle surface of the shell can be located with respect to a local convected coordinate system, 9 . Global coordinates, x and global displacements a v x components, "J7* of P are given by, V = N "fy (4.78) TP = N Tty (4.79) M M where x* are the nodal coordinates of the ^ node in the I v N tion t h direction at configura- "Ify is nodal displacement defined in a similar manner and is the shape function of the A/ node. The shape functions for WIQS are, th = J(l + 0i)(l + 02) N l J V ^ J (1-0*) (1 + 02) (4.80) \(i-e )(i-e ) N* = l 2 JV = i ( l + 0i)(l-02) 4 Substituting Eq. (4.78) into Eq. (4.2) yields the position vector of the point P, V T = N" x\ (4.81) v g i Differentiating Eq. (4.81) yields the covariant surface base vectors, a ~ d6 a Substituting from Eq. (4.81) we have, dN • M "*a = -gga- *\ gi = X v g i (4.82) Three Dimensional Finite Element Analysis / 4.4 85 «3 K2 Si F I G U R E 4.4 where "a* = Typical Finite Element. dN N (4.83) 1/ » Substituting Eq. (4.83) into Eq. (4.7) yields the surface covariant metric tensor, dN* dN" I> t vt (4.84) Three Dimensional Finite Element Analysis / 4.4 86 The incremental virtual work equations are considered next. Substituting Eqs. (4.29) and (4.30) into Eq. (4.77) yields, +i V « gafipX d xi l °A d V ^Zd9 d9 l (4.85) 2 °A where from Eq. (4.22), d°A = V®ad9 d9 . For external loads acting on a finite element l 2 surface, the external virtual work expression is given by, = JJ V 2 ext 8u V a~d9 d9 i 2 1 2 (4 2A where V is the » component of the externally applied traction vector, 5u* is the * comth t h ponent of the virtual displacement vector and A is the deformed area in configuration C. 2 2 The displacement components can be calculated as, u' = V - V = x (4.87) i Since the external load vector comprises only normal pressure, the traction vector can be written as, 2*» _ 2 „ 2 .« p e t = (4.88) 3 where p is the magnitude of the normal pressure and e is the i 2 2 t n 3 component of the unit normal vector, e in the deformed configuration. Substituting Eqs. (4.87) and (4.88) into 2 3 Eq. (4.86) yields, JJ P 4 6x s/*a'd9 2 ex i 2 = Wt 2 l d9 (4.89) 2 Substituting Eq. (4.89) into Eq. (4.85) we have, JJ 1 [ dx = JJ 4 l d xJ % l 2 [ Jj d x*d xi x x d9 -±JJ 2 [«*•' p e ] v ^ d * 2 2 3 °A 1 ,i % l lap 6x ri ( ap) d la d x* l 0 ad9 d9 1 2 (4.90) Three Dimensional Finite Element Analysis / 4.4 87 To simplify the incremental virtual work equations and thus facilitate computer programming, the following finite element matrices are introduced, A" 0 0 N 0 0 1 [<?] = \P } = 0 a where N N = dN /d9 M a 2 JV 0 0 0 N 0 N 0 0 l 0 0 2 0 0 N •• 3 0 2 • • • N* 0 N, 0 0 N' 0 0 N% N* n 0 0 0 N} 0 N, 0 2 0 a (4.91) a 2 (4.92) Q m—J - l a 0 N 0 \ I — / "V 1 i / 2 i/„ 3 \ " ~ i/_3 t / _ l i/_2 i/_3 2 r f - \l l l 2 2 2 • • •A \ \ I . x x x x x x X x x (4.93) uc\2 (4.94) {%.} = { "a "a , 2 3 In view of Eqs. (4.92), (4.93) and (4.94), Eq. (4.83) becomes, (4.95) {%,} = [P ] TO a Also Eq. (4.84) becomes, "*afi = TO' [^] (4.96) TO where (4.97) [Pa ] = [P.] fa] T fi Substituting Eqs. (4.91) and (4.93) into Eq. (4.78) and the result into Eq. (4.90) we have, 5 ^"a/?) l a^pA (S*) g V //(* r °A ! p Q , / e ] ^adOUQ 2 2 3 2 A II i <fiad9 d9 \ l 2 (4.98) Three Dimensional Finite Element Analysis / 4.4 88 where 6t'=l > l 7 term 1 0, otl otherwise t h a t Since the equation holds for arbitrary variations of the nodal coordinates 6 £ j , we have dC 2 //(' °A ( af)) d l a0p\ l(1 d r a? i\\) de } z 2 d^ 1 J J V de 2 (4.99) °A ext Eq. (4.99) represents the incremental equations of equilibrium for an element. However, it is impossible to evaluate thefirstterm on the l.h.s. of Eq. (4.99) because the quantities associated with the deformed configuration C are not known apriori. The situation is 2 further complicated by the fact that the external loading is deformation dependent. The problem of converting the integrals over A to integrals over A have been discussed by 2 l Oden and Key (1970), Oden (1972) , Larsen and Popov (1974), and Bathe, Ramm and 6 Wilson (1975), among others. To evaluate the integral, the following approximation will be introduced in this thesis. Expanding the external virtual work expression in Taylor series about the deformed configuration C we have, l d( W ),j l W 2 —W 1 ext ext + + (higher order terms) (4.100) ext Examining the second term, we have d( W ) l i.e. ext e J _ \[to K p ) ' ^ 4 . , , ,£fV5) 6 cit. he, p. 242 ^ l c ^ (4.101) Three Dimensional Finite Element Analysis / 4.4 89 Observe that this second term represents a stiffness-like quantity, arising from the followerload character of the pressure loading and hence contributes to the system tangent stiffness matrix. Unfortunately, this is an unsymmetric contribution. This unsymmetric load stiffness term been has identified by several researchers, for example, Oden (1969), Oden and Key (1970), Haug and Powell (1972a), Larsen and Popov (1974), Bathe, Ramm and Wilson (1975), Nagarajan and Popov (1975), and Bathe (1982). Most of them however, neglected this unsymmetric stiffness as the pressure loads in their analysis were small, and the resulting loss of accuracy is within acceptable limits. They also pointed out that the extra computational effort to solve the unsymmetric equations is quite substantial. However, for the case of membrane structures, the wind pressures encountered can be very large. Thus the unsymmetric stiffness matrix is included in the analysis presented here. Eq. (4.100) then becomes, W 2 « / / [ p <?„ V ] v l ext 3 W d9 + I jj % 2 Q i I 1.4 (4.102) + Substituting Eq. (4.102) into the incremental equations of equilibrium depicted in Eq. (4.99) and after some re-arranging, yields, i a/3p\ d( pA) lfl r E -JJ = ff[ pQ» *] l 1 * VA + + P-Q^r y/la+ ^d^d9 -Jf 2 i P ^-dT^r d6 d9 l e ad9 d9 0 d? 1 l 2 2 (4.103) °A Eq. (4.103) depicts the linearized incremental equilibrium equations derived using a virtual work method. They agree with those obtained by Haug and Powell (1972a) whose Three Dimensional Finite Element Analysis / 4.4 90 derivation is based on a force equilibrium method. The unsymmetric stiffness terms were 7 not considered in their analysis and this, they suspected, was the cause of the convergence difficulties they encountered. Eq. (4.103) can be expressed in the following compact form, {K\? K1<?-KW}C> = + P I - (4.104) F I where K.ffi, Kjp , Kjp are the elastic, the geometric and the unsymmetric stiffness matrices respectively and are defined as, (E) K IJ d iff°A u la dH C«p\) d l afi X ( af}) r P 0 1 (G) K IJ 0 ad9 d9 1 ad9 d9 l (4.105) 2 (4.106) 2 °A K IJ -ffo -JJ \ Q i I d ^ W P V V ^ ) * e V a + ,t n ^ ' ^ ^ p ~dqr Va+ I^nV d ^ d9 d9 l 2 (4.107) ^-d^J- P e and Pj, Fj are the load vectors due to external and internal actions respectively. They are defined as, = ff[ pQu <&] l V^Zd9 d9' (4.108) l l A l l„a0 d ( « / ? ) la Fi \fia~d9 d9 1 2 (4.109) °A The difference (Pj — Fj) represents an unbalanced load term. To simplify Eqs. (4.105), (4.106) and (4.109), we observe the following relationships from Eq. (4.96), 7 The force equilibrium method presented by Hang and Powell ( 1972a) b derived for a specialised case. Bathe (1982) has indicated a more general force equilibrium approach to formulate the linearized incremental equilibrium equations (cit. Joe, p. 309-310). / 4.4 Three Dimensional Finite Element Analysis 91 (4.110) and QI^T^J where P pu is the (J, J ) a t h = ( °PKL + P<XPLK) I P 6 term of [P p] • Thus we have, a * ff ^ ^ l Kl M 1 C^ a (>X = \jj u (PafJKL + Pa/UK) (PpXMN P XNM)S}S^V^Ld9 d9 + K (4.111) 6J l (4.112) 2 p (P*f*KL + PaPLK) ti SfV^dd 1 d9 (4.113) 2 °A Pi = i /j ^ (PaPKL + PafiLK) ti (4.114) ^ 0,1 in which C P l ALI X is given by Eq. (4.63) as, IQOPPX and n l _ Et \„ap 2(1 + a/3 1 f)X , l„aX \„fip , n) I (-5L) V> (4.115) is given by Eqs. (4.23) and (4.62) as, (4.116) i.e. To evaluate Eq. (4.107), we make the following observations: the first term is zero if the pressure distribution within an element is assumed constant. The second term is evaluated from Eq. (4.14) as, v -A l " yfa3 3 l *s9u which upon differentiation yields, d£ l J l a e5 3$ a ^ 1 3 e 3 3 e 9kl From Eqs. (4.10) and (4.95) we have, a ^ 1 (4.117) Three Dimensional Finite Element Analysis / 4.4.1 92 where % = e jkl Vgg PiJ + 4 P{j} li (4.H8) l The third term in Eq. (4.107) is evaluated from Eq. (4.20) as, djyfia') 1 = da , p l = d{ a ) l Q afi Substituting from Eq. (4.110) we get, gl£j ~ ~Y ° C [PapKJ + PapJK) (4.119) Thus in view of Eqs. (4.117) and (4.119), Eq. (4.107) becomes, IJ K = JJ + QiJ P l l ^ \ 3 3 S j e e e ^3 1 «/»l^ ( P p a fl K J + e 3 P , a J K 3 ) SjJ V^d9 d9 l 2 (4.120) Having developed the governing finite element equations for the incremental analysis, we shall discuss next, the solution strategy. 4.4.1 EQUILIBRIUM ITERATION To control the error arising from the linearization of Eq. (4.72), it is necessary to iterate within each load increment. Iteration is continued until the difference between the internal virtual work, evaluated at the computed static and kinematic variables, and the external virtual work satisfy a prescribed tolerance. Thus from Eq. (4.72), the error at the fc th iteration is, Error*** = W 2 ext - jj n°^ 2 6 ( $) 2 7 d°A (4.121) °A unbalanced load Comparing the r.h.s. of Eqs. (4.77) and (4.121), we see that the former represents the unbalanced load before the computation of the solution increments while the latter represents the unbalanced load after the solution process. To reduce this second unbalanced Three Dimensional Finite Element Analysis / 4.4.1 93 load, it is necessary to perform an iteration within the load increment. The equation to be solved at the fc iteration is, th Jj i<7 ^ A e J J 6e Q A d°A + JJ n ? 6*$ l afi a d°A = % - JJ °A ezt V**" ) 6 1 Substituting the various finite element arrays and after some manipulations similar to that in deriving Eq. (4.103), we have, {// [i (^ <^^) +l f r $ $ ) ] J * " l d9 de l = JJ[ pQu ei] l l l a 2 + 2 Vhide de* i ff Uap I d( « p) . X ]] \ d*£' d i d^^ I 2 * 1 2 I ad9 d9 0 K °A ' i (4.123) 2 J where the incremental displacements are updated as follows, £j(fc) £J{k-i) = + £J(k) (4.124) A It should be noted that at the iteration k — 1, Eq. (4.123) corresponds to Eq. (4.103). Expressing Eq. (4.123) in a compacted form we have, jrg>} A W {Kff + K<ff where the arrays P , = ( P _ rf-i) . (4 125) K jf , fljj* and A J ^ are given by Eqns. (4.108), (4.112), (4.113) { I} and (4.120) respectively. The array Fj ~^ is obtained by substituting Eqs. (4.110) and k (4.111) into Eq. (4.123) to yield, if _ 1 ) = J / / (V [ ? (P 2 A F I K L + P A F I L K °A + {P pKL + a P pLK)ti^ji \f J a k ) 6} ^ ^~ad9 d9 l 2 (4.126) Three Dimensional Finite Element Analysis / 4.5 94 Observe in Eq. (4.125), the stiffness matrix is not updated at each iteration. Only the unbalanced load is updated. Hence we see that the procedure described so far corresponds to a modified Newton-Raphson iteration. If the stiffness matrix is updated at every iteration in Eq. (4.125), namely as follows, then we have the standard Newton-Raphson iteration. 4.5 C O M P U T E R I M P L E M E N T A T I O N The computer implementation of the incremental finite element equations derived for the structural analysis is discussed here. Only the main programming features of are described. The program NFEAM3, is written in Fortran IV and its program layout is depicted in Figure 4.5. The program exhibits a dynamic allocation of memory. Here all arrays are defined as one-dimensional and stored sequentially in one master array, A . An address system automatically sets up the location of each array in A. Thus each individual array is variably dimensioned to the exact size required for a problem. Also any changes to the overall problem, say in the number of elements, can be effected in NFEAM3 by simply, a single change in A . At the data input phase the user can, for standard shell configurations, invoke the generation of an appropriate finite element discretization scheme from the library. For instance, to discretize a circular cylindrical shell, the user types in the macro command CYLI, followed by the number of nodes, types of boundary conditions etc. CYLI then proceeds to generate the nodal point and element information. In a similar fashion, the macro command SPHE generates the necessary discretization information for a spherical shell. For non-standard shapes, the user inputs his own nodal point and element information following the macro command USER. A summary of the main input data is shown in Tabie 4.1. Except for those variable names that are self-explanatory, the names will be defined as they are introduced in this section. For example, the program execution code MODIN, controls the direction of flow in the program. It is defined as follows, Three Dimensional Finite Element Analysis / 4.5 95 { .EQ. 1 - Data Check Only .EQ. 2 - Program Execution .EQ. 3 - Restart Option All fatal memory storage and data errors are diagnosed by subroutine DOCTOR, which then prints the appropriate diagnostic messages as indicated by the control code, NURSE. Variable Name Description TITLE Title of Job NELEM Number of Elements NNODE Number of Nodes NVAEL Number of Variables per Element NSYMM Symmetry Code for Equation Solver, SOL ISOLN Solution Technique Code ICONV Convergence Criterion Code MODIN Program Execution Code NSTEP Number of Load Steps MITER Maximum Number of Iterations PINT Internal Pressure in Shell TENS Initial Uniform Tension per unit length RADI Principal Undeformed Radius of Shell TOLN Prescribed Solution Tolerance BETA Line Search Parameter Table 4.1 8 Fortran Variable Names used in the Main Data Input The most efficient strategy for the storage of global matrices is the skyline scheme. These matrices are evaluated using a four-point Gauss quadrature. Figure 4.6 shows how the elements of the global stiffness matrix S K are stored in the one-dimensional This terminology is adopted from the frontal solution package given in Irons and Ahmad (1980). Three Dimensional Finite Element Analysis 96 / 4.5 START > begin clock • check RESX\RT code > set up all arrays w.r.t. a l-D master array I READ INPUT DATA • • • • • • DATA ANALYSIS job title & control information nodal coordinates element connectivities boundary conditions material information loading parameters • check memory allocation for data input • diagnosis by DOCTOR • data organisation • printing I I Fatal Errors START LOAD INCREMENT ITERATION CONSTRUCTION TECHNIQUE CONSTRUCT GLOBAL QUANTITIES • compute element stiffness matrix • compute element unbalanced load vector • assemble elemental quantities into global I I > check memory allocation for solution arrays i imposition of boundary conditions during assembly > skyline storage for global stiffness matrix Fatal Errors PREDICTIONS • selection of initial solution SOLUTION OF EQUATIONS CORRECTIONS • choice of solution technique via ISOLN > use of a compacted unsymmetric solver, SOL • improvement algorithm • line search parameter No No < NITER > NITER INVOKE RESTART FACILITY UPDATE LOAD INCREMENT OUTPUT STRAINS & STRESSES print results & stop clock > print program statistics > print errors (if any) 1 ' compute and print strains i compute and print stresses F I G U R E 4.5 Structure of NFEAM3. Three Dimensional Finite Element Analysis / 4.5 97 array A . All zeros inside the skyline are included as they can become non-zero in the solution process. If only symmetric matrices are encountered, then only the diagonal terms and the upper triangular part need to be stored. A precise procedure to address the elements of S K in A is required. This is achieved by establishing the column height array L M and the connectivity array K L D . Thus we see that the efficiency of this skyline 9 scheme lies in the fact that no elements outside the skyline are stored or processed in the computation stage. Also to keep S K as small as possible, the imposition of boundary conditions is carried out during the assembly process. The resulting system of equations is solved by a symmetric-unsymmetric skyline storage solver, SOL. A solution technique code ISOLN, is provided in NFEAM3. Its definition is, .EQ. 1 - Triangularization of Stiffness Matrix ISOLN .EQ. 2 - Modified Newton-Raphson Iteration .EQ. 3 - Standard Newton-Raphson Iteration > .EQ. 4 - Automatic Selection of N-R Iteration At each step of the Newton-Raphson iteration, the unbalanced forces are computed. The stiffness matrix is either kept constant or re-computed, depending on ISOLN. For the case of ISOLN = 4, the program starts with the modified Newton-Raphson iteration and progresses automatically to the standard Newton-Raphson iteration when the results have attained 80% of the prescribed convergence criterion. The external loading is applied incrementally, as defined by the number of steps parameter, NSTEP. If NSTEP is too small (i.e. large load increments) the solution process could diverge. The proposed incremental solution strategy is based on an iterative technique. For this strategy to be effective, proper convergence criteria for the termination of the iteration run should be adopted. As discussed in Chapter 3, Bergan and Clough (1972) have argued that a displacement-based criterion is preferable to the unbalanced load criterion. Another feasible criterion is that based on energy. Here the increment in internal 9 The scheme used here is taken from Dhatt and Tonzot (1984), p. 223-225 Three Dimensional Finite Element Analysis / 4.5 98 Skyline Profile [_^21_ %22^ ^23 0 i I SKD = SKU = SKL = (#11, ^35 0 J -^44 -^45 (#21, 0 0 >. ^ -^53 -^54 -^55 -^56 Actual Stiffness Matrix (Unsymmetric). #221 ^33> #44> ^ 5 5 » ^ 6 6 ) Of °i ^Se) # 3 2 » -Ktl 1 Of -^43J ^531 ^ 5 4 > ^621 0> 0> ^65) (#12, # 2 3 ) ^14> 0> ^ 3 4 ) ^35> ^45> ^26> (b) F I G U R E 4.6 -^43 0 j (a) Q. J ^26 ^ < 0 1 K33 K$2 1 1^41 0 Arrays for Storing Elements of S K . Storage Scheme used for Global Stiffness Matrix. Three Dimensional Finite Element Analysis / 4.5 99 energy during each iteration is compared to the initial internal energy increment. Bathe 10 and Cimento (1980) have reported their experiences with this type of energy convergence criterion. In this thesis, the displacement convergence criterion is used. A convergence criterion code ICONV, sets up the appropriate vector norm to be used in the convergence test. It is defined as follows, { .EQ. 1 - Maximum Norm Convergence Criterion .EQ. 2 - Euclidean Norm Convergence Criterion A line search parameter BETA, is also provided in NFEAM3 to assist solution convergence. As the idea has been explained quite extensively in Appendix D, no further elaboration will be pursued here. Eq. (4.124) is thus modified to ZW = £ * J J ( _ 1 ) + BETA<*) A f ' W (4.128) Also, incorporated in NFEAM3 is a RESTART facility. If for some reasons, the solutions are not converging (for example, the load increment is too large) then the RESTART option is automatically invoked through MODIN. The solving process is halted and all essential information pertaining to the previously computed equilibrium configuration like the global stiffness matrix and the unbalanced load vector, is stored in a preassigned file. An error message is then printed by DOCTOR and the program execution terminated. The user can now examine his output and make necessary adjustments to enhance solution convergence. When the program execution is restarted by the user, the program begins iteration from the previously obtained equilibrium configuration rather than from the undeformed configuration. Finally, the coupled program, INFEAM3 is obtained through the synthesis of the 3-D boundary element program and the 3-D finite element program. We shall now present some of the results obtained in the uncoupled finite element analysis as well as the coupled boundary element-finite element analysis. This is defined as the amount of work done by the unbalanced loads on the displacement increments. Three Dimensional Finite Element Analysis 4.6 N U M E R I C A L / 4.6.1 100 RESULTS As in the two-dimensional analysis, this section is presented in two parts: the noninteractive analysis where the results of the finite element program, NFEAM3 are discussed, and the interactive analysis in which results of the coupled boundary element-finite element program, INFEAM3 are reported. 4.6.1 NON-INTERACTIVE ANALYSIS We shall now examine the performance of the proposed finite element formulation. Two membranes loaded laterally by pressure are considered in this exercise. The first example consists of a pressure-loaded, pretensioned rectangular membrane while the second example involves a long cylindrical membrane. To avoid unnecessary complications, the fluid-structure interactions have been ignored in the analysis. The computed membrane deformations and stresses are compared with published results. This then enables us to gauge the accuracy of the proposed structural analysis. Example 1: Pressure-Loaded, Pretensioned Square Membrane An investigation into the inflation of a square membrane with dimensions axaxt and fixed at the edges is described here. The membrane, held initially taut and flat by a uniform pretension per unit length T, is loaded by internal pressure, p. The deflection response and incremental tension are then computed and compared with the analytical solutions of Irvine (1976, 1981). He analysed a membrane with the following parameters: o = 60 in t = 0.575 x I O in - 3 E = 15 x 10 psi 6 (4.129) r = 10lb/in p = 0.035 psi In the analytical solutions, the membrane is subjected to vertical loads per unit area while the membrane studied here is loaded by normal pressure. This discrepancy can be ignored as the curvature of the membrane under load is small. Three Dimensional Finite Element Analysis 1 / 4.6.1 101 6 C 4 5 B 2 3 A a a F I G U R E 4.7 Discretization Scheme for the Square Membrane. The membrane is discretized with 36 elements as shown in Figure 4.7. To be consistent with the analytical solutions, a slightly non-uniform mesh is used. The loads are applied in 10 steps and a convergence tolerance of 1% on displacements is employed. To start the iteration, an assumed initial profile is prescribed. The line search scheme is also invoked in the solution iteration to assist convergence. The average line search parameter, Pav P e r l ° d increment for convergence is found to be 0.459 while the average number of a iterations required to achieve convergence at each load step is 5.3. Three Dimensional Finite Element Analysis / 4.6.1 Points 2 102 Deflection in inches 1 Linear Theory 0.164 0.323 0.366 Irvine (1976) 0.135 0.264 3-D F E M 0.142 0.268 Table 4.2 on Membrane 3 4 5 6 0.691 0.798 0.928 0.299 0.565 0.653 0.759 0.296 0.572 0.656 0.763 Comparison of Deflection Response for the Square Membrane A comparison of the deflection response is presented in Table 4.2. For convenience, the units used by Irvine is adopted here. The deflected profile is sketched in Figure 4.8. It is apparent that good agreement exist between the numerical and the analytical solutions. Also shown in Tabie 4.2, are solutions obtained from the linear analysis. It is observed that the deflections predicted by the linear analysis are always larger than those obtained by the nonlinear analysis. The centerline deflections (i.e. along the planes located at x — 30 in. or y = 30 in.) for the various stages of load increments are depicted in Figure 4.9. A nonlinear stiffening effect is observed as the load increment approaches the 100% pressure. Tabie 4.3 compares the theoretical predictions of the incremental tension in the membrane obtained by Irvine (1976, 1981) with the results from the finite element analysis. In the theoretical analysis, Irvine has assumed that the incremental tension say, in the i-direction, is a function of the y-coordinates only. That is to say, the incremental tension is constant along any line parallel to the z-axis. To be consistent with this assumption, the finite element results listed in Table 4.3 represent the averaged values. For example, to obtain the incremental tension along a line parallel to, say, the z-axis, we take the average of the incremental tension evaluated at several points along that line. No distinction between the horizontal components of the incremental tension and the true incremental tension is made as the curvature of the membrane under load is small. Three Dimensional Finite Element Analysis F I G U R E 4.8 / 4.6.1 103 Defected Profile for the Square Membrane. DISTANCE IN INCHES F I G U R E 4.9 Defection Response of the Centerline with Increasing Load Increments: Square Membrane. Three Dimensional Finite Element Analysis / Incremental Locations Tension in lb/in Table 4.3 4.6.1 104 in Membrane A B C Linear Theory 0.0 0.0 0.0 Irvine (1976) 0.210 1.878 3.667 3-D F E M 0.242 1.955 3.770 Comparison of Tension Increments for the Square Membran The relationship between the central deflection, 5 m a x with the internal pres- sure, p for different Poisson's ratios is shown in Figure 4.10. Like a cable net structure, a nonlinear stiffening of the membrane is observed with the increasing pressures. The possibility of solution divergence for such hardening structures has been mentioned by Cook (1981). He suggested the use of a line search parameter, (3 < 1.0, to assist convergence. Pursuing in this manner, it has been found for the examples tested that, 0.35 <j3 < 0.50 (4.130) To check that the small-strain assumption is not violated, we monitor the largest principal strains in the membrane. The results are displayed in Figure 4.10 for varying levels of internal pressures. The table in the graph shows that the membrane is still exhibiting small-strain deformations. The variation of membrane tension increments with internal pressure is shown in Figure 4.11. The tension increments here are the maximum x - direction tension. It can be seen that the tension increases with the pressure at an increasing rate. This is consistent with the observation that the membrane is stiffening as the pressure increases. Also at larger Poisson's ratios, the tensions at corresponding pressures, are larger. Three Dimensional Finite Element Analysis FIGURE 4.11 / 4.6.1 105 Variation of Maximum Tension Increments with Internal Pressures for the Square Membrane. Three Dimensional Finite Element Analysis / 4 106 Example 2: Long Cylindrical Membrane The long cylindrical pneumatic membrane introduced in Chapter 3 will be resolved here using the three-dimensional finite element representation. The cylinder is initially circular and is subjected to wind blowing from the broadside. This problem was solved in Chapter 3 by assuming the cylinder to be sufficiently long so that it reduces to a two-dimensional problem. We will again assume the cylinder to be sufficiently long. But instead of treating the problem as two-dimensional, we will retain the three-dimensionality by considering a slice of the cylinder of unit width and employing the appropriate boundary conditions to simulate the state of plane-strain. The following parameters are prescribed in the analysis: Ro = 10 m, p = 1.25 kg/m , 3 a p = 0.20 (4.131) The discretization scheme used is depicted in Figure 4.12. The loads are applied in 10 steps and a convergence tolerance of 0.1% on displacements is sought. No line searches are required in the solution iteration. The average number of iterations required to achieve convergence per load step is 5.8. F I G U R E 4.12 Discretization Scheme for the Cylindrical Membrane. Three Dimensional Finite Element Analysis / 4 107 Examples 1 and 2 of Chapter 3 are used in this investigation. The results obtained are compared with solutions from the two-dimensional finite element analysis and from Uemura (1972). For Example 1, the deformation profiles and the hoop stress resultants are shown in Figure 4.13 and TaWe 4.4 respectively. Reasonably good agreement is observed. Again, for convenience, the metric units used by Uemura are retained here. U E M U R A (1972) & 2 - D F E M 3-D FEM F I G U R E 4.13 Membrane Deformations for Example 1 (£, = 1.91, C = | x IO ,). -4 Hoop Stress Resultants, Table 4.4 in Kg/m Uemura (1972) 2-D F E M 3-D F E M 818.7 864.1 886.3 Comparison of Hoop Stress Resultants for Example 1 Three Dimensional Finite Element Analysis / 4 108 For Example 2, the deformation profiles are sketched in Figure 4.14 and the hoop stress resultants listed in Table 4.5. Reasonably good agreement is seen between the twofiniteelement solutions. This suggests that the approximations introduced in the strain-displacement relations for the two-dimensional analysis are valid, as such approximations are not used in the three-dimensional analysis. Also, for reasons mentioned in Chapter 3, the agreement with Uemura's solutions is poor. UEMURA ° ° 0 0 F I G U R E 4.14 (1972) 2-D FEM Membrane Deformations for Example 2 Hoop Stress Resultants, Table 4.5 (£ = 3.40, C = l x 10 j. _4 a in Kg/m Uemura (1972) 2-D FEM 3-D FEM 1227.3 1271.9 1299.3 Comparison of Hoop Stress Resultants for Example 2 Three Dimensional Finite Element Analysis 4.6.2 I N T E R A C T I V E / 4.6.2 109 ANALYSIS Having studied the performance of the three-dimensional finite element program, we consider next, the interactive analysis of wind-loaded pneumatic membrane structures. The structure analysed is the membrane roof of the BC Place Stadium. BC Place Stadium Roof The stadium, located in Vancouver, Canada, was completed in 1983. It has an air-supported roof and covers an area of 4.05 hectares. The roof is a composite structure made up of orthotropic membranes and cables. In this analysis, we will assume homogeneous and isotropic membranes without any cable reinforcements. A schematic depiction of the stadium together with its flow domain is shown in figure 4.15. Note that the roof's minor axis is in the i-direction which is the wind direction, while the major axis is in the y-direction. The following parameters are used in the flow analysis: L = B = 1600 m, H = 800 m, p = 1.20 kg/m 3 a (4.132) and for the roof, Li = 182.9 m, B = 224.3 m, x H = 26.8 m, x H = 33.5 m 2 (4.133) E = 0.690 GPa, fi = 0.30, t = 7.62 mm, p = 345 Pa a Table 4.6 summarizes the discretization schemes used in the analysis. The boundary element grid adopted for the flow domain is the same as the 'cube' used for the sphere (see Figure 2.11) and for the stadium, is shown in Figure 4.16. Note that symmetry about the x-z plane is exploited. The finite element grid for the roof is identical to the boundary element grid shown in Figure 4.16(B). The following specifications are common for the analyses presented here: all pressure loads are applied incrementally in 10 steps, a convergence criterion of 1% is imposed on the displacements and the line search scheme is invoked to assist convergence in the finite element iteration. F I G U R E 4.15 Stadium and Flow Domain. Three Dimensional Finite Element Analysis / 4.6.2 112 Boundary No. o f B o u n d a r y No. o f No. o f F i n i t e Section Nodes Singularities Elements Roof 50 6 38 Wall 60 4 Flow Domain 171 14 - Total 281 24 Table 4.6 - 38 Discretization Schemes F i n a l l y , t h i s r e p o r t i s d i v i d e d i n t o three stages: (a) Inflation (a) Inflation analysis (b) Flow model investigation (c) Interactive analysis Analysis To initiate the solution process f o r m e m b r a n e deformations due t o wind, t h e i n f l a t i o n p r o b l e m i s first s o l v e d . T h i s i n v o l v e s c o m p u t i n g t h e s t r u c t u r a l r e s p o n s e o f t h e r o o f d u e only t o i n t e r n a l p r e s s u r e , p . T h e u n d e f o r m e d r o o f g e o m e t r y w h i c h i s s t r e s s - f r e e , a is s h o w n i n Figure 4.17. 11 T h e d a t a set p e r t a i n i n g t o t h i s u n d e f o r m e d r o o f is t a b u l a t e d i n Appendix E. I t s h o u l d b e e m p h a s i z e d h e r e t h a t t h e f o l l o w i n g t e r m i n o l o g y i s a d o p t e d : undeformed p r o f i l e refers t o t h e s t a t e p r i o r t o i n f l a t i n g t h e r o o f , w h i l e initial p r o f i l e r e f e r s t o t h e s t a t e a f t e r t h e r o o f is f u l l y i n f l a t e d . T h e i n t e r n a l pressure is a p p l i e d i n c r e m e n t a l l y u n t i l t h e p r e s c r i b e d f u l l l o a d is r e a c h e d . A l s o a t each load increment, a 'local' iteration is set u p t o c o m p u t e t h e corresponding converged profile. T h i s converged profile t h e n becomes the starting profile f o r t h e n e x t l o a d i n c r e m e n t i t e r a t i o n . T h e a v e r a g e n u m b e r o f ' l o c a l ' i t e r a t i o n s i s 3.1 w h i l e t h e a v e r a g e l i n e s e a r c h p a r a m e t e r e n c o u n t e r e d f o r t h e o v e r a l l s o l u t i o n is 0.873. Recall that self-weight of the membrane has been neglected and hence the 'np' position of the profile. Three Dimensional Finite Element Analysis F I G U R E 4.17 / 4.6.2 113 Undeformed and Initial Roof Profiles due to Internal Pressure. Figure 4.17 depicts the deflected shape due to the full pressurization. The maximum height of the roof increased by 6.8% or 1.98m. The distribution of membrane stress resultants along the minor axis of the roof is shown in Figure 4.18. As expected, 12 N x appeared to be fairly constant while N the apex, the shear stress resultant, N xy y increases towards the apex of the roof. At is zero which is also expected, as it is the point of double symmetry. A plot of the distribution along the major axis is not shown as it is felt that no new information is conveyed. But it suffices to point out that the N x roles are interchanged. Hence, along this axis, N y is fairly constant while N x and N y increases towards the apex. The stress resultants at the apex are also tabulated in Table 4.7. Also listed are the stress resultants computed by the elliptic-paraboloid shell theory for the same pressure load. We note that the roof deformed symmetrically about both the minor and major axes but with very different radii of curvature in these two axes. Because of the latter, the top of the roof is best modeled as an elliptic-paraboloid shell. This analytical approach provides an immediate check of the finite element solutions shown in Table 4.7. Good agreement is observed. The stress resultants are actually evaluated at the gauss points but this difference is small and thus neglected. Three Dimensional Finite Element Analysis / 4.6.2 114 o to F I G U R E 4.18 Stress Resultants along the Minor Axis for the Inflation Analysis. Stress Resultants Elliptic-Paraboloid Shell Theory kN/m Finite Element Solution kN/m N 44.10 47.01 N 33.22 36.60 0.0 0.0 x y Table 4.7 Stress Resultants at the Apex of the Roof, for the Inflation Problem Three Dimensional Finite Element Analysis / 4.6.2 115 (b) Flow Model Investigations The wind, as mentioned previously, is modeled by potential flow theory. To see how closely this approach models real flow, the following non-interactive analysis is performed. The roof is subjected to two types of wind pressures: wind pressures computed by potential flow theory and wind pressures experimentally measured by Parkinson (1980). In the latter, a rigid-roof model of the stadium was built and the pressure measurements were carried out in the boundary layer wind tunnel at the University of British Columbia. The set-up is shown in Figure 4.19. For the purpose of this study, a wind blowing at 25 m/s in the direction of the minor axis of the roof is considered. Figure 4.20(A) compares the two resulting wind pressures along the minor axis of the roof. We see in the case of the model roof, that the wind pressure comprises only the suction type. This implies the occurrence of flow separation on the roof. 13 However, in the case of the potential flow model, the wind pressure on the roof consists of both positive and negative (or suction) pressures. The occurrence of these conflicting pressure predictions is attributed to the built-in limitations of the potential flow theory. But towards the top of the roof, the pressures obtained by these two approaches do agree in sign and differ at most by 10%. The wind pressure along the major axis which is perpendicular to the flow direction, is shown in Figure 4.21(A). Here, as expected, both flow models predicted suction pressures over the entire length of the axis. The deformed profiles computed by the two flow models, along the minor axis are shown in Figure 4.20(B). Note that, in the case of the potential flow model, there is inward deflection of the roof at two locations: the windward and the leeward side of the roof. This is not observed in the deformed profile predicted by the real flow model. Again, this contradiction in the deflected shape is attributed to the inherent limitations of the potential flow theory. However, as expected from the corresponding pressure plot, As reported in Parkinson (1980),flowseparation over a large portion of the roof was encountered. Three Dimensional Finite Element Analysis F I G U R E 4.19 / 4.6.2 Model Test of BC Place Stadium (Parkinson (1980)). 116 Three Dimensional Finite Element Analysis / 4.6.2 117 POTENTIAL FLOW REAL FLOW (PARKINSON (1980)) (A) INITIAL PROFILE POTENTIAL FLOW REAL FLOW (PARKINSON (1980)) o. S o EH K m — •« « * i ft n r; n r i Ti .t .-; »• - •« „ o I—I w a (B) F I G U R E 4.20 Pressure and JZoof Profiie for the Minor Axis. the agreement in the vicinity of the roof's top is very good. The maximum elevation of the roof predicted by the two flow models differ by only 2.4%. Figure 4.21(B) shows the deformed profiles for the major axis. Here, the deformed profiles agree very well over the entire length of the axis. The corresponding distributions of membrane stress resultants along the minor axis are shown in Figure 4.22. Again, we see for the case of the potential flow model, that N x is fairly constant while N y increases towards the apex. Similarly, N xy is zero at the apex. For the real flow model, the distribution is not much different from that of the potential flow model. But, because the symmetry in the flow direction is destroyed, N xy Three Dimensional Finite Element Analysis CM « ID W Ui w / 4.6.2 118 POTENTIAL FLOW REAL FLOW (PARKINSON (1980)) CO I" i Pi «r> cu i (A) O. OO INITIAL o CO PROFILE POTENTIAL FLOW REAL FLOW (PARKINSON (1980)) o • .*» ."T . T . T .T.T.t VT . T . ^ . ^ , ^ ^ . t (B) F I G U R E 4.21 Pressure and Roof ProGle for the Major Axis. is now no longer zero at the apex. The stress resultants at the roof's apex are tabulated in Table 4.8. Also shown are stress resultants computed from the elliptic-paraboloid shell theory for a pressure loading corresponding to the potential flow model. We see from the table, that the stress resultants predicted by the two flow models differ by at most 16% which is within acceptable engineering practice for these structures. Thus, despite the inherent limitations of the potential flow theory, the controlling stresses and deformations are obtained to acceptable accuracy for design. It is on this ground that we feel it is justified to use the much simpler potential flow theory to model the wind. Three Dimensional Finite Element Analysis / 4.6.2 119 s F I G U R E 4.22 Stress Resultants along the Minor Axis for Flow Model Investigations. Stress Resultants *y N xy Table 4.8 Elliptic-Paraboloid Shell Theory kN/m Potential Flow Model kN/m Real Flow Model kN/m 74.17 91.53 78.91 50.52 58.87 69.49 0.0 0.0 4.37 Stress Resultants at the Apex of the Rigid Roof (V = 25m/s) Three Dimensional Finite Element Analysis / 4.6.2 120 (e) Interactive Analysis With the determination of the converged roof profile due to internal pressure in part (a), the wind is now 'turned on' the structure. The deformed profiles are solved initially, at low wind speeds and then at higher speeds until the desired speed is attained. In this way, the wind loads are applied incrementally, with the converged profile determined at each load step until the desired wind loads are reached. Solution is completed when the converged profile corresponding to this desired wind load is obtained. (a) Wind Speed, m/s (b) Height of Apex, m (c) Displacements in m, (%) (d) Displacement Increments, m 15 31.176 2.366 (8.2) 2.366 25 32.214 3.404 (11.8) 1.038 35 32.727 3.917 (13.6) 0.513 Table 4.9 Maximum Roof Elevations (Initial Height — 28.81m) Some typical results of the interactive analysis for wind speeds of 15, 25 and 35 m/s in the direction of the minor axis are presented in Figures 4.23 - 4.25. Plots (A) and (C) show the roof pressure along the minor and major axes of both the initial and deformed roof while plots (B) and (D) show the corresponding roof profiles. We see that the wind loads corresponding to the flexible roof are larger than those from the rigid roof. Also, from these figures, one could see that the membrane deformations are not very sensitive to the interactive effects in which the flow around the deformed structure is computed. Table 4.9 shows in column (b) the height of the apex, (c) the apex displacement away from the initial state, and (d) the increment in this displacement. Clearly, we observe a sharp decrease in the displacement increments with increasing wind speeds, due to the nonlinear stiffening effects of the membrane. Three Dimensional Finite Element Analysis / 4.6.2 121 RIGID MODEL FLEXIBLE MODEL w OT OT o. INITIAL PROFILE DEFORMED PROFILE o (B) o o. RIGID MODEL FLEXIBLE MODEL OT OT W (C) o. S .S 00 to INITIAL PROFILE DEFORMED PROFILE CM (D) F I G U R E 4.23 Pressure and Roof Profile: (A) and (B) refer to the Minor Axis, (C) and (D) refer to the Major Axis ( V = 15 m/s). Three Dimensional Finite Element Analysis / 4.6.2 122 INITIAL PROFILE DEFORMED PROFILE OP) F I G U R E 4.24 Pressure and Roof Profile: (A) and (B) refer to the Minor Axis, (C) and (D) refer to the Major Axis (V — 25 m/s). Three Dimensional Finite Element Analysis / 4.6.2 123 FLEXIBLE MODEL (D) FIGURE 4.25 Pressure and Roof Profile: (A) and (B) refer to the Minor Axis, (C) and (D) refer to the Major Axis (V = 35 m/s). Three Dimensional Finite Element Analysis F I G U R E 4.26 / 4.6.2 124 Membrane Deformations along the Minor Axis for Varying Loading Ratios. The net membrane deformations away from the initial state are shown in Figure 4.26. We see that the roof deforms almost symmetrically about the major axis. This pattern of deformation appears to contradict similar results obtained in the twodimensional interactive analysis in which the flow past an infinitely long cylinder was studied. Recall in those results the symmetry in the deformed profile was destroyed. This 14 difference in the deformation pattern is due to pressure loads being more symmetrical for the three-dimensional structure than the two-dimensional one. Note that in the case of the two-dimensional structure, the flow goes up and over the structure whereas for the threedimensional structure, the flow goes up and over, and as well as around the structure. Thus 1 4 see Figures 3.11 and 3.12 Three Dimensional Finite Element Analysis / 4.6.2 125 this 'extra dimension' creates less disturbance to the flow and therefore less disturbance to the pressure field. That is why, the pressure coefficient for a circular cylinder is -3.0 15 while for a sphere of the same radius, it is only -1.25. With less disturbance, the roof has a greater tendency to remain in the symmetrical mode. The foregoing discussion is not applicable to the case of real flow. Recall in §(b) that because of flow separation in the real flow model, the profile becomes asymmetric. Instead of looking at the deformed roof from its minor and major axes plots, it is useful to present some three-dimensional views. These are given in Figure 4.27 for the three selected wind speeds. The initial roof is shown in red while the deformed roof is plotted green. The distributions of membrane stress resultants along the minor axis are shown in Figure 4.28. We see that N the apex. At that point, N xy x is generally constant while N y increases towards is almost zero (see Table 4.10). The stress resultants at the Apex of the roof are also tabulated in Table 4.10. The maximum shear stress resultant for the corresponding wind speeds are listed in Table 4.11. Note that they occur towards the edges of the roof. Stress Resultants kN/m N x N y Table 4.10 Wind Speed V = 15 m/s Wind Speed V = 25 m/s Wind Speed V = 35 m/s 54.82 81.25 110.35 41.35 53.78 70.84 0.10 0.13 0.28 Stress Resultants at the Apex of the Flexible Roof defined as the ratio of the maximum pressure to the dynamic pressure of theflow,j p U 2 a Three Dimensional Finite Element Analysis / 4.6.2 F I G U R E 4.28 127 Stress JResuitants along the Minor Axis: (A) for V=15m/s, (B) for V=25m/s, (C) for V=35m/s. Three Dimensional Finite Element Analysis Wind Speed, m/s 15 25 35 Table 4.11 4.7 / 4.7 Maximum N , xy kN/m 11.40 20.48 31.65 128 Locations in Roof X Y Z -41.0 -65.1 -65.1 71.1 82.3 82.3 17.0 4.62 4.66 Maximum Shear Stress Resultants and their Locations S U M M A R Y A three-dimensional interactive analysis of pneumatic membrane structures is presented. Theflow,modeled as an inviscid fluid, is analysed using a boundary element approach. The structure is analysed with a finite element technique. The finite element method is based on a continuum mechanics-based incremental formulation using a convected coordinate description. Two computer programs are developed for the analysis: NFEAM3 for the non-interactive analysis and INFEAM3 for the interactive analysis. The performance of NFEAM3 are assessed on two pneumatic membranes: a pressure loaded, pretensioned square membrane and an infinitely long cylindrical membrane. Three-dimensional representation is used in the latter. Comparisons with published results are made and good agreement is attained for the deformations and stresses. As an example of a three-dimensional pneumatic membrane structure, the air-supported roof of the BC Place Stadium is analysed. To justify modeling the wind by potential flow theory, the roof deformations are solved with wind pressures computed using this theory and from a model test in a wind tunnel. The results show that, despite the built-in limitations of the potential flow model, the critical deformations and stresses are obtained to reasonably good accuracies. However, it is felt that the inclusion of real flow effects such as flow separation, would be desirable before the procedure is put to engineering use. Three Dimensional Finite Element Analysis / 4.7 129 Interactive results for the deformed roof and stress resultants for wind blowing at the three selected wind speeds are also presented. It is shown that if the effects offluid-structureinteractions are ignored, the wind loads computed can err on the nonconservative side. CHAPTER 5 Summary and Conclusions The aim of this thesis is to develop a method for the analysis of wind-loaded pneumatic membrane structures, taking into accountfluid-structureinteraction. The wind, assumed to be inviscid, is modeled using a boundary element technique with moving singularities. The membrane structure, assumed to exhibit large deformations but small-strains, is analysed using a finite element scheme. For two-dimensional structures, the finite element method used is based on the Galerkin formulation, while for three-dimensional structures, a continuum mechanics-based incremental formulation using a convected coordinate description is employed. The accuracies of the foregoing numerical methods are assessed separately on flow and structural problems with known solutions. To make the analysis interactive, the boundary element program is coupled to the finite element program. This is done in the following manner: given an assumed membrane profile, the wind loads are calculated from the flow analysis using the boundary element program. These wind loads are then fed into the finite element program which performs the structural analysis for a new membrane profile. Convergence is attained when an equilibrium profile consistent with the wind loads results. Interactive analysis is performed on the these pneumatic membrane structures: an infinitely long cylindrical membrane as an example of two-dimensional and the air-supported roof of BC Place Stadium, as a three-dimensional structure. It is shown that the wind loads predicted by the flexible model can be larger than those of the rigid model. 130 Summary and Conclusions / 5.1 131 Hence, a pneumatic membrane structure analysed using a rigid model approach which is the current engineering practice, may not yield a conservative design. The conclusions regarding this thesis can be conveniently grouped into three categories. These are, the flow analysis, the structural analysis and the interactive analysis. Each will be discussed in turn. 5.1 T H E F L O W ANALYSIS In the flow analysis, we are interested in obtaining the wind loads acting on the membrane structure. This is done via a boundary element technique. The singularities in this formulation are placed on an auxiliary boundary located outside, encompassing the domain of the problem. Also, because they are allowed to move, the number of singularities used can be significantly less than the number of boundary elements employed in the discretization without any apparent loss of accuracy. The auxiliary boundary is determined using a minimization criterion. Two criteria have been used: the least squares criterion and the Galerkin criterion. In both criteria, the usual weighting functions have been augmented with functions defined as derivatives with respect to the singularity locations. By allowing the singularities to move freely until they attain optimum positions, a highly accurate and yet adaptive method for computing wind loads is realized. However, one pays a price for this feature. The resulting boundary element equations are nonlinear, with the nonlinearity occuring in the determination of the coordinates of the singularities. From the plot of the loci of the singularity movement, we see that the singularities generally tend to move away from the domain of the problem. Also, as the solution approaches convergence, minimum movement is observed. The method is illustrated with several examples involving the Laplace's equation in two and three-dimensions. They are, deflection of a stretched membrane,flowpast a rigid circular cylinder, heat conduction in a rectangular prism and flow past a rigid sphere. To accurately gauge the performance of this boundary-based numerical method, 1 as opposed to thefiniteelement method which is a domain-based numerical technique Summary and Conclusions / 5.2.1 132 the following boundary conditions are prescribed in the examples: Dirichlet boundary conditions, Neumann boundary conditions, and mixed boundary conditions. Comparisons with known solutions of these examples show the computed values on the boundary as well as in the domain of the problem to be in excellent agreement. 5.2 T H E S T R U C T U R A L ANALYSIS In the structural analysis we are interested in obtaining the membrane deformations, the wind loads acting on the membrane being assumed given. The membrane is assumed to be perfectly flexible in bending but possessing finite extensional rigidity. The self-weight of the membrane has also been ignored in comparison with the inflation pressure. 5.2.1 TWO-DIMENSIONAL P N E U M A T I C M E M B R A N E S T R U C T U R E S The structural analysis of two-dimensional pneumatic membrane structures is formulated on an infinitely long, wind-loaded cylindrical membrane. A small strain-finite deformation membrane theory is derived. The resulting second-order integrodifferential equation for the membrane profile is solved via a Galerkin finite element procedure using Newton-Raphson iteration with line^search. The line-search parameter encountered in the analysis, ranges between 0.15 to 0.30. To check the accuracy of the finite element results at low wind speeds, a perturbation technique is developed. The perturbation parameter, e used here, is defined as the ratio of a characteristic membrane displacement to the diameter of the cylinder. The foregoing theory is then applied to solve a cylindrical membrane subjected to a cosine wind load distribution. Good agreement in the deflected shape and stresses between the pertubation solutions and the finite element solutions are obtained. For moderate to high wind speeds, the finite element results are compared with published soutions. Again, good agreement is observed. Summary and Conclusions 5.2.2 133 / 5.3 THREE-DIMENSIONAL P N E U M A T I C M E M B R A N E S T R U C T U R E S An incremental formulation using a convected coordinate description is developed for the analysis of three-dimensional pneumatic membrane structures. The metric tensor for any deformed configuration is derived and from this, a consistent measure of deformation is obtained. The incremental continuum mechanics equations are then used to develop the governingfiniteelement equations. Thesefiniteelement equations are solved using NewtonRaphson iteration. An 'automatic' line-search algorithm is also incorporated into the iterative process. This algorithm using either a quadratic or a cubic backtracking strategy, determines an appropriate step-size increment along a given Newton direction so as to yield an acceptable next iterate. The average line-search parameter experienced is 0.42. To examine the performance of the proposed formulation, a pressure-loaded, pretensioned square membrane and an infinitely long, wind-loaded cylindrical membrane are analysed. In the latter, three-dimensional representation is used in the analysis. Comparison with published results is made, and good agreement is attained in both the examples. For the square membrane, the variations of central deflections and maximum tension increments, with internal pressures are presented. The results show a nonlinear stiffening of the membrane with the increasing pressure, just as in a cable-net structure. 5.3 T H E I N T E R A C T I V E ANALYSIS Having discussed the performance of the boundary element and the finite element programs, interactive results for wind-loaded pneumatic membrane structures are then generated by coupling these two programs. As a measure of loading on the membrane, a loading parameter £„ is presented. It is defined as the ratio of the dynamic pressure to the internal pressure. Thus the membrane experiences larger deformations with increasing £ . a For two-dimensional structures, an infinitely long cylindrical membrane is considered. Interactive solutions pertaining to the influence of the loading ratio, membrane deformations and the variation of membrane stress resultants, on with the Summary and Conclusions / 5.4 134 internal pressure, p are presented. It is seen that the additional stress due to wind blowing 0 at a specific velocity, changes rapidly at smaller p , and then decreases at increasing a p . This is because the membrane structure has becomes so 'rigid' that there is little 0 fluid-structure interactions. Consequently, the stresses also converge to the perturbation solutions. Also presented is a comparison of the final wind pressure coefficients predicted by the rigid membrane model and the flexible membrane model. The result shows that the coefficients from the flexible model can be larger than those from the rigid model. As an example of a three-dimensional pneumatic membrane structure, the air-supported roof of the BC Place Stadium is analysed. Three wind speeds are considered: 15, 25 and 35 m/s. The wind as mentioned, is modeled as an inviscid fluid. To justify this approach, the roof is solved for two sets of wind pressures. The first set are those computed by the potential flow theory and the second, from a model test measurement in a wind tunnel. The critical stresses and deformations in the roof predicted by these two flow models agree very closely. Some typical interactive results of the deformed roof and stress resultants are also presented. Just as in the two-dimensional analysis, the results show that the wind loads of the flexible roof can be larger than those of a rigid roof. Thus it is important to consider thefluid-structureinteraction effects. 5.4 R E C O M M E N D A T I O N S F O R F U R T H E R S T U D Y Some specific recommendations for further studies are: a) The handling of anisotropic membranes and the inclusion of cables in the case of the three-dimensional analysis. b) A turbulent flow model for the wind to account for Reynolds number effects such as flow separations. c) Dynamic response of pneumatic membrane structures. d) Dynamic stability analysis and static stability analysis to account forflutteringand wrinkling of membranes respectively. e) Experimental investigations. 135 References ADKINS, J.E. and RIVLIN, R.S. (1952) Large Elastic Deformations of Isotropic Materials Part IX, Phil. Trans. A, Vol. 244, No. 505 AHMADI, G. and GLOCKNER, P.G. (1983) Ponding Instability of Insulated Imperfect Cylindrical Membranes J. Struct. Div., ASCE, Vol. 109, p. 297-313 AHMADI, G. and GLOCKNER, P.G. (1984) Effect of Imperfections on Ponding Instability Eng. Mech. Div., ASCE, J. Vol. 110, p. 1167-1173 ALLGOWER, E. and GEORG, K . (1980) Simplicial and Continuation Methods for Approximating Fixed Points and Solutions to Systems of Equations SIAM Review Vol. 22, p. 28-85 ARMIJO, L. (1966) Minimization of Functions having Lipschitz-Continuous First Partial Derivatives ASCE (1979) Pacific J. Math. Vol. 16, p. 1-3 State-of-the-Art report on Air Supported Structures Task Comm. on Air Supported Struct, of the Comm. on Metals, of the Struct. Div. 95 p. Finite Element Procedures in Engineering Analysis BATHE, K . J . (1982) Prentice Hall, New Jersey, 735 p. BATHE, K.J., and CIMENTO, A.P. (1980) Some Practical Procedures for the Solution of Nonlinear Finite Element Equations J. Comp. Mthds. App. Mech. Engrg., Vol. 22, p. 59-85 BATHE, K.J., R A M M , E., WILSON, E.L. (1975) Finite Element Formulations for Large Deformation Dynamic Analysis Int. J . JVum. Meth. Eng., Vol. 9, p. 353-386 B E C K E R , E.B. (1966) A Numerical Solution of a Class of Problems of Finite Elastic Deformations Ph.D. Thesis, Univ. Calif., Berkeley BENZLEY, S.E. and K E Y , S. W. (1976) Dynamic Response of Membranes with Finite Elements J. Eng. Mech. Div., ASCE, p. 447-460 BERGAN, P.G. and CLOUGH, R.W. (1972) Convergence Criteria for Iterative Process Journal, Vol. 10, No. 8, p. 1107-1108 BIRD, R.B., STEWART, W.E. and LIGHTFOOT, E.N. (1960) Transport Phenomena AIAA John Wi- ley and Sons, New York, 780 p. BIRD, W.W. (1967) The Developments of Pneumatic Structures, Past, Present and Future Proc. IASS, 1st Int. Coll. Pneumatic Structures, Stuttgart, p. 1-9 BLUE, J.L. (1976) Solving Systems of Nonlinear Equations Computing Science Technical Report #50 Bell Lab., Murray Hill, N.J., 14 p. BLUE, J.L. (1978) Robust Algorithms for Solving Systems of Nonlinear Equations puting Science Technical Report #67 BREBBIA, C A . and DOMINGUEZ, J . (1977) lems Com- Bell Lab., Murray Hill, N.J., 10 p. Boundary Element Methods for Potential Prob- Appl. Math. Modeling Vol. 1, p. 372-378 136 On Kupradze's Functional Equations for Plane Harmonic Prob- C H R I S T I A N S E N , S. (1967) Function Theoretic Methods in Differential Equations lems (Eds. R.P. Gilbert and R.J. Weinacht), Pitman, London, p. 205-243 C O O K , R.D. (1981) Concepts and Applications of Finite Element Analysis 2nd Ed., John Wiley k Sons, N.Y., 537 p. A Faster Modified Newton-Raphson Iteration C R I S F I E L D , M.A. (1979) ods in Applied Mechanics and Engineering, Computer MethVol. 20, No. 3, p. 267-278 An Arc-Length Method including Line Searches and Accelerations C R I S F I E L D , M A . (1983) Int. J. Num. Meth. Engrg., Vol. 19, p. 1269-1289 .Numerical Methods for Unconstrained Optimization and Nonlinear Equations Prentice Hall, New Jersey, 378 p. D E N N I S , J.E. Jr. and S C H N A B E L , R.B. (1983) D E N T , R.N. (1971) Principles of Pneumatic Architecture The Architectural Press, Lon- don 236 p. D H A T T , G. and T O U Z O T , G. (1984) G.) Beams, Plates and Shells D O N N E L L , L.H. (1976) F L U G G E , W. (1973) The Finite Element Method Displayed (transl. CANTIN, John Wiley & Sons, 509 p. Stresses in Shells McGraw Hill, New York 453 p. 2nd Ed., Springer Verlag 525 p. U.S. Pavillion at Expo 70 features Air Supported Cable Roof Eng., ASCE, p. 48-50 G E I G E R , D.H. (1970) G O L D S T E I N , A.A. (1967) Constructive Real Analysis Harper ic Row, N.Y., 178 p. Large Elastic Deformations G R E E N , A.E. and A D K I N S , J.E. (1970) Civil 2nd Ed., Oxford Univ. Press, London 324 p. G R E E N , A.E. and Z E R N A , W. (1968) Theoretical Elasticity 2nd Ed., Oxford Univ. Press, London 457 p. H A N , P.S. and O L S O N , M.D. (1986) Static Analysis of Pneumatic Structures Loaded by Wind submitted to IASS Int. Symp. Membrane Structures and Space Frames Osaka, Japan H A N , P.S., O L S O N , M.D. and J O H N S T O N , R.L. (1984) lation with Moving Singularities A Galerkin Boundary Element FormuJ. Eng. Comput. Vol. 1, p. 232-236 Large Elastic Deformations of Thin Rubber Int. J. Eng. Sc., Vol. 5, p. 1-27 H A R T - S M I T H , L.J. and CRISP, J.D.C. (1967) Membranes Finite Element Analysis of Nonlinear Membrane StrucRep. No. UC-SESM 72-7, Univ. Calif., Berkeley H A U G , E. and P O W E L L , G.H. (1972a) tures H A U G , E. and P O W E L L , G.H. (1972b) tures Finite Element Analysis of Nonlinear Membrane Struc- Proc. IASS, Pacific Symp. Tension Struct. Space Frames, Tokyo and Kyoto, p. 165-175 Numerical Properties of Integral Equations in which the Given Boundary Values and the Sought Solutions are defined on Different Curves Comput. Struct. Vol. 8, p. 199-205 H E I S E , U. (1978) 137 HOPPER, M.J. (ed.) (1973) Harwell Subroutine Library Catalogue Division, A.E.R.E. Harwell, England Theoretical Physics HO-TAI, S., JOHNSTON, R.L. and MATHON, R. (1979) Software for Solving Boundary-value Problems for Laplace's Equation using Fundamental Solutions Tech. Rep. 136/79, Dept. Comp. Sc., Univ. Toronto IRONS, B . M . and A H M A D , S. (1980) Techniques of Finite Elements Ellis Horwood Ltd., Chichester, 529 p. IRVINE, H.M. (1976) Analytical Solutions for Pretensioned Cable Nets Div., ASCE, Vol. 102, p. 43-57 IRVINE, H . M . (1981) Cable Structures J. Eng. Mech. MIT Press, Cambridge, Mass., 259 p. JOHNSTON, R.L. and FAIRWEATHER, G. (1984) The Method of Fundamental Solutions for Problems in Potential Flow Appl. Math. Modeling Vol.8, p. 265-270 JOHNSTON, R.L., FAIRWEATHER, G. and H A N , P.S. (1984) The Method of Fundamental So- lutions, an Adaptive Boundary Element Method, for Problems in Potential Flow and Solid Mechanics Div., Proc. Fifth ASCE Specialty Conf., Eng. Mech. Laramie, Wyoming, p. 500-504 KLINGBEIL, W.W. and SHIELD, R.T. (1964) Some Numerical Investigations on Empirical Strain Energy Functions in Large Axisymmetric Extensions of Rubber Membranes Z. angew. Math. Phys., Vol. 15, No. 608, p. 608-628 KUNIEDA, H., Y O K O Y A M A , Y . and A R A K A W A , M . (1981) Structures Subject to Wind 867 Cylindrical Pneumatic Membrane J. Eng. Mech. Div., ASCE, Vol. 107, p. 851- K U P R A D Z E , V. (1967) On the Approximate Solution of Problems in Mathematical Physics Russ. Math. Surveys Vol. 22, p. 59-107 [Uspehi Mat. Nauk, Vol. 22, p. 59-107] LANCHESTER, F.W. (1917) Patent 119,339 LANCHESTER, F.W. (1938) Span Trans. Sessions 1938-1939, Manchester Assoc. Engineers, Machester, England p. 59-97 LARSEN, P.K. and POPOV, E.P. (1974) Large Displacement Analysis of Viscoelastic Shells of Revolution Computer Methods in Applied Mechanics and Engineering, Vol. 3, p. 237-253 LUKASIEWICZ, S.A. and GLOCKNER, P.G. (1983) Ponding Instability of Air Supported Spherical Membranes with Initial Imperfections J. Nonlinear Mech., Vol. 18, p. 1-9 McCONNELL, A . J . (1957) Applications of Tensor Analysis Dover Publications, N.Y., 318 p. M A L C O L M , D.J. and GLOCKNER, P.G. (1978) Collapse by Ponding of Air Supported Membranes J. Struct. Div., ASCE, Vol. 104, p. 1525-1532 M A L C O L M , D.J. and GLOCKNER, P.G. (1981) Collapse by Ponding of Air Supported Spherical Caps J. Struct. Div., ASCE, Vol. 107, p. 1731-1742 138 MATHON, R. and JOHNSTON, R.L. (1977) The Approximate Solution of Elliptic Boundary Value Problems by Fundamental Solutions SIAM J. Num. Anal. Vol. 14, p. 638-650 MATTHEIS H. and STRANG, G. (1979) The Solution of Nonlinear Finite Element Equations Int. J. Num. Meth. Engrg. Vol. 14, p. 1613-1626 MORE, J.J. (1977) The Levenberg-Marquardt Algorithm, Implementation and Theory Numerical Analysis, (ed. G.A. Watson ), Lecture Notes in Mathematics, Springer-Verlag MORRIS, N.F. (1975) Cable Syatems Shock and Vibration Computer Programs, Review and Summaries, SVM-10, tion Monograph Series (eds. W. and B. Pilkey), The Shock and Vibra- MORRIS, N.F.(1981) Wind-Structure Interaction of Flexible Roofs Proc. Symp. Long Span Roof Structures, Struct. Div., ASCE, p. 297-312 N A G A R A J A N , S. and POPOV, E.P. (1975) Non-Linear Dynamic Analysis of Axisymmetric Shells Int. J. Num. Meth. Eng., Vol. 9, p. 535-550 ODEN, J.T. (1966) Analysis of Large Deformations of Elastic Membranes by the Finite Element Method Proc. IASS, Int. Congr. Large-span Shells, Leningrad ODEN, J.T. (1967) Numerical Formulation of Nonlinear Elasticity Problems Div., ASCE, Vol. 93, p. 235-255 J . Struct. ODEN, J.T. (1969) Discussion of "Finite Element Analysis of Non-Linear Structures" Struct. Div., ASCE, Vol. 95, p. 1379-1381 Finite Elements of Nonlinear Continua ODEN, J.T. (1972) J. McGraw Hill 432 p. ODEN, J.T., K E Y , J.E. and FROST, R.B. (1974) A Note on the Analysis of Nonlinear Dynamics of Elastic Membranes by the Finite Element Method Comp. <fe Struct., Vol. 4, p. 455-452 ODEN, J.T. and KUBITZA, W.K. (1967) Numerical Analysis of Nonlinear Pneumatic Structures Proc. IASS, 1st Int. Colloq. Pneumatic Structures, Stuttgart, p. 82-107 ODEN, J.T. and SATO, T. (1967) Finite Strains and Displacements of Elastic Membranes by the Finite Element Method Int. J. Solids Struct., Vol. 1, p. 471-488 OLIVEIRA, E.R. (1968) Plane Stress Analysis by a General Integral Method Div., ASCE Vol. 94, p. 79-101 OLSON, M.D. (1982) Notes CIVL 538: J. Eng. Mech. Advanced Topics in the Finite Element Method, Class University of British Columbia OLSON, M.D. and H A N , P.S. (1984) Interactive Analysis of Wind-Loaded Membranes Proc. 2nd Int. Conf. Num. Methods Nonlinear Problems, Barcelona, p. 501-512 Iterative Solution of Nonlinear Equations in Academic Press, N.Y., 572 p. ORTEGA, J.M. and RHEINBOLT, W.C. (1970) Several Variables OTTO, F. (1962) Zugbeanspruchte Konstruktionen and Published as Tensile Structures Ullstein Fachverlag, Berlin. Translated MIT Press, Cambridge, Mass. 320 p. 139 P A R K I N S O N , G.V. (1980) Wind Pressure Distributions on a Model ofB. C. Place Amphitheatre for Phillips Barratt Engineers <fe Architects, 23 p. P A T T E R S O N , C. and SHEIK, M.A. (1981) Regular Boundary Integral Equations for Stress Proc. Third Int. Sem. Boundary Element Methods Analysis Irvine, Ca. p. 85-104 On the Use of Fundamental Solutions in Trefftz Method for Potential and Elasticity Problems Proc. Fourth Int. Sem. Boundary Element Methods in Eng. Southampton, p. 43-57 P A T T E R S O N , C. and SHEIK, M . A . (1982) A Hybrid Method for Nonlinear Equations Numerical Methods for Nonlinear Equations, (eds. Rabinowitz, Gordon and Breach), p. 87-114 P O W E L L , M.J.D. (1970) ROSS, E.W. Jr. (1970) Large Deflection of an Inflated Cylindrical Tent ASME, No. 69-WA/APM-9, SIMONS, J.W. and P O W E L L , G.H. (1982) Structures S T E V E N S , H.H. (1942) U E M U R A , M . (1972) p. 845-851 Solution Strategies for Statically Loaded Nonlinear UCB/EERC-82/22, Univ. Calif. Berkeley Air Supported Roofs for Factories Architectural Record, Membrane Tension and Deformation in Cylindrical Pneumatic Struc- tures Subject to Wind Loading Space Frames, ZIENKIEWICZ.O.C. (1977) 787 p. J. Applied Mech., Proc. IASS, Pacific Symp. Tension Struct. Tokyo and Kyoto, p. 199-210 The Finite Element Method 3rd Edit., McGraw Hill, London, APPENDIX A Boundary Element Equations The governing boundary element equations for the Laplace's equation are presented in this appendix. The presentation is divided into two sections: results for 2-dimensional space and results for /-dimensional space. These equations are derived for the most general boundary condition, namely, the mixed boundary condition. By the appropriate selection of P and 7 in these results, the equations for Dirichlet and Neumann boundary conditions 1 are realized. The Laplace's equation in /-dimensional space is given by, du 2 du du 2 du 2 2 where u, in potential flow theory, represents the velocity potential. Eq. (A.l) is solved numerically. This is because the boundary enclosing the domain where this equation is defined, is highly irregular. Substituting Eqns. (2.3c) and (2.9) into Eq. (2.8) yields, M N + i w or £{ ^ " })],. ' c ( t x ) < r = o S/L fil(' *ww ri w+ r, 1 See E q . (2.3) 140 Boundary Element Equations / A 141 Boundary Element Equations for 2-Dimensional Space The normalized Green's function for the 2-dimensional Laplace's equation is given by Eq. (2.9a). Substituting this equation into Eq. (.4.2) results in, M N EE/ .=i j=i C d(/dt k dr = 0 (a(x) + tf(t ,x)) lk y \<k<N (A.Z) . r where x — (xi, x ) are the coordinates of the boundary node t = (t , t ) are the coordinates of the singularity 2 lk k ^,x) = 2k -log V ^ - + ^ ^ - x ^ 2 R + >/'?* + <2* 2(<ijfc - x ) d$ dt 2(r [t 2k lk - x) 2 2Jt - x f + (t x 2k - x )2 k k (AA) 2k 2t2k + R^Tq + \ +q k 2 2 (x)C H(t ,x) = /3(x)C $(t ,x) + k 2*life + x 1 t k . 3xi . k (hk- l)-g^ X (<lJk-Xl) + (*2Jfc-^2) 2 2 k , . + (t2k-X ) dx 2 ^ 2 The expressions (dx /dn) and (dx /dn) are the direction cosines with the x and x axes l 2 x 2 respectively. Boundary Element Equations for l-Dimensional Space; I > 3 The normalized Green's function for the /-dimensional Laplace's equation is given by Eq. (2.9b). Substituting this equation into Eq. (A.2) results in, f M f(t*,x) ' C ds/dt N EE/ k lk dT = 0 Ckd$/dt 2k K C d$/dt , k lk \{a(x) + n{t x)) j) l<k<N, l>3 (A.5) Boundary Element Equations / A 142 where x = (s l f x , 1 3 , . . . , X/) denotes the b o u n d a r y nodes 2 (*iJti hk> hki • • • > '/it) = denotes the singularities (2-0 ?(t ,x) = [(t fc - X ) + (< * " * ) + • • • + ft* " lfc I/) ] 2 2 2 t 2 2 2 2-J -(/-2)(t 1 J t -z ) 1 //2 [(<lJfc - l l ) + (*2Jb " *2? + ' •' + ftfc 2 I/) ] 2 ^l-l -(* - 2)(* - x ) 2 2jk 3/ 1/2 [ ( * i * - x , ) + (« ib - * 2 ) + • • • + ft* ~ x / ) ] 2Jk 2 2 (A.6) 2 2 (' ~ 2)* Jfc + 2 /-1 +4 + -- + 'f*[* + >/'?* 4 + - + 4] + -(/-2)(t / J t -x,) 1/2 (hk ~ * i ) + ihk ~ * ) + • • • + ft* - * / ) ] 2 2 2 2 + l-l (/ - 2) (x)C 7 ^(ft x) = /?(x)C (t x) + Al ibf ft* X Jfc " *i) 2 + ihk - x . ) ^ - + ft* - x )-^ 2 (hk ~ * ) + ••• + 2 2 2 l 2 , ..ii t - ft* dxi + • • • + ft* - x ) As before, the expressions (dxt/dn), (3x /9n),..., (dxi/dri) the x , i -.f/2 jfc> Q x,) 2 n are the direction cosines with axes respectively. Eqns. (A.3) and (A.5) constitute the necessary (/ -I- 1)7Y boundary element equations for the solution of the unknown coefficients and singularity locations. These Boundary Element Equations / A 143 equations are however, non-linear and are solved using an alternating algorithm with line search. The integration is performed numerically, using the Gauss quadrature rule. In the program, the computation of the jacobian is obtained from exact expressions. Some simplifications are achieved when the problem exhibits symmetry. Having obtained the coefficients and the singularity locations, the solution is then given from Eq. (2.4) as, N u(x)=£C y ? (A.7) (ty,x) 3=1 Since the boundary element method presented here is used for the flow analysis, it will be useful to give the following expressions for flow velocity and pressure. The flow velocities in x lf x , x , . . . directions, namely, v , v , v , . . . are given by, 2 z Xl x% X3 du dXi That is, 2fty-x.) • = 1,2 (l-2)(t -x ) i] 3=1 i 772 * = 1> 2,3, />3 The resultant velocity, V(x) is given by, V(x) = dxi dx (A.8) dxi 2 The Pressure, P(x) is given by Bernoulli's equation which is, P(x) = P (x) + ip{(Vb(x)) - (K(x)) } 2 Q 2 (A.9) where PQ(X) and K (x) are the datum pressure and free stream velocity respectively. 0 APPENDIX B Stiffness Matrix G tJ The stiffness matrix G,y is defined in Eq. (3.40) as, r > (B.l) dP OWj where P, is given by Eq. (3.53) as, P =P l = R A<p j l + C{N+ - R ) - 2 2 Pa R[Va ~^ P w ) | (B.2) Since Pj = P , the subscript can be dropped. That is, 2 From Eq. (-B.2) we have, dPi dP dwj du)j + = RA<j>^R^C + 2+ R(Pa-PwY C{2N -3R )<j> J N} dtVj > -p„)1 dR\ 3R(Pa " (5.3) g Pa N+ J a»j J where from Eq. (3.46) we have, dw 1 2C($ 2 L, \ Jn K n=l dR dwi - dR [9 (w\ + tof) + &«>i«>2 " 4CR p &<}>] + R A<t>} z x N t 2 a { E\i [^i^i + ^ 1 \ + &«>i«>2] + Atfun + u; + 2CJ2 p ) J J (B.4) 2 2 a 144 Stiffness Matrix / B 145 dN. dw 2 2C7($2 f» = l dR dw 2 _ dR * t in which Q± and Q are defined as, d w 2 3 + A<f> 2 A<f> -6 2 92 = 3A<(> 3&<p From Eq. (3.47) we have, dR dwi (2fl I y^ (2& 2 2 + je ) 2 2 + k\ ' (B.6) +zy 2 dR dw \(2Q 2 2 +2Q + 4Q R Z 2 2 + Soc ) 2 2 QZJ + Zl< (2fi2 £2)2 + in which R\ and £ are defined as, 2 Q = * l = (*b + « l ) ) (B.7) *2 = (*0 + «>2) The stiffness matrix G,-y is obtained by substituting Eq. (B.Z) into Eq. (5.1) which yields, ON* 3ici dN* "I 5 tea J a« dR ai? n (B.8) dR where ? = R &<p{C + 2 t 7 = RA+ | . + C 2 2 } - 3R ) Pa - 3RiPa N ~ } Pw) APPENDIX C Perturbation Analysis For the case of small wind loads, the finite element solutions will be checked with results from a perturbation analysis. The formulation for the perturbation analysis will be given here. IntegrodiSerential Equation From Eqs. (3.22) and (3.24), the integrodifferential equation for the structural analysis of a cylindrical pneumatic membrane subjected to arbitrary wind is given by, 2Rw + (wl + w 2 H> Iw+Rt) + 2Rw = 2R 2 1+ R(Pa~Pu>Y ( C . l ) C{N+-Rp )a where f{^(rvl = ~ + w )+w 2 + CR }d<p 2 Pa 5 (C2) CfR(<f>)d<j> *i For small wind loads, the radius of curvature does not change appreciably. Hence we will assume that the radius of curvature remains constant. i.e. R(<f>) = constant = R (C.3) Eqs. (C.l) and (C.2) then become, + + ™ ) + » = R [ l + C[N+ - R ) 2 Pa - ~^ R{Pa P w ) (CA) 146 Perturbation Analysis / C 147 where *2 Boundary Conditions The cylinder is held down at the two edges and this yields, xv(i ) = w($ ) = 0 l (C.6) 2 Non-Dimensional Form The integrodifferential equation in Eq. (CA) can be recast in a non-dimensional form as g(<j>), by defining an unknown characteristic displacement w , given by, t w(<f>) = tw, g(<*>) (C.7) Also, the wind pressure p (<f>) can be non-dimensionalized as $(<f>), by defining a characw teristic wind load p, as follows, PM) (C.8) = P.$(<f>) When a wind is prescribed, it fixes the value of p, and is thus a known quantity. Substituting Eqs. (C.7) and (C.8) into Eqs. (C.4) and (C.5) we get, 1 + C(N and *1 <j> - R ) Pa R(p -P.t(<l>)) N. a (C.9) Perturbation Analysis / C 148 Perturbation Parameters At low wind speeds, the deformation of a cylindrical membrane is small compared to its diameter. Hence a perturbation parameter e can be defined as follows, •=£ < > cn The above definition is obvious as it multiplies the nonlinear terms depicted in Eq. (C.9). Also another small parameter, TJ can be defined with respect to the pressure loads. At low wind speeds, the wind pressure is small compared to the internal pressure. Hence we have, s = (C.12) Pa It will be shown later that n = e in order for terms of the same order of magnitude to balance. This would result in the following expression, (C.13) Perturbation Solutions In view of Eqs. (C.ll) and (C.12), the integrodifferential equation becomes, 9" + e[9 ,2 1 + C{N+ - R ) + 9 ]+9=lr 2 Pa Rp (l-V$(<t>)) - a (C.H) where 1 (C.15) C($ - $ o 2 Perturbing g(<f>) and as follows, 9{<t>) = 9o(<f>) + Witt) + e 92{4>) + 0{e ) 2 N <j> = N + N Q e + eN 2 l + 0(e ) 3 2 3 (C.16) (C.17) Perturbation Analysis / C 149 Substituting Eqs. (C.16) and (C.17) into Eqs. (C.14) and (C.15), and collecting terms yields, (So + %) + <dl + 9i + 9? + <7o) + e {92 + 92 + HA 2 = ^[{i + W - J +*> {CN 2 R p J - * } + - O e{cjv + ^ ( ^ + ^ ) ) } 1 ((ft) " ft 2 [ f t ft - ( f t ) - ft + ( ( f t ) 2 + 2<7 0i) + • • • 2 3 2 } + *>{CN fe^w) - 3 ft) (C.18) } + •••] and N + 0 + e N + • • • = Rp + e < 2 2 a 9o + 9 )d<f>}+--- (C19) 2 2 We can set the first term on the r.h.s. of Eq. (C.18) to zero and this yields, 1 {l + C(N -R )-!fc}=0 0 (C + ^)(N 0 Since Pa - R ) =0 Pa (C + ^ ) ^ 0 (N - Rp ) = 0 then N = Rp i.e. 0 Q a (C.20) a Observe that Eq. (C.20) corresponds to the linear solution for the hoop stress of a pressureloaded cylinder. We can thus expect to get a very good approximation to the hoop stress for the nonlinear problem by taking only a few terms. Substituting Eq. (C.20) into Eq. (C.18) we have, (go + go) + + 9i + g? + go ) + e {g$ + g + 2 2 2 + 2<7ffi) + • • • 0 = 3 { C » i + ft + J f ( « } + | [CN + ft - (ft)' - ft^(*)} 9 +J 1 {^3 + ft + (ft) 3 - 2ft ft Similar conclusion can also be arrived from Eq. ((7.19) + [(ft) - ft] J j W ) } 2 + ••• (C.21) Perturbation Analysis / C 150 Note that on the r.h.s of Eq. (c.21), to balance terms of the same order of magnitude, we must have, »? = e (C.22) Eq. (C.21) then becomes, (flfo + 9o) + e(9i +9i + 9? + 9o ) +s {92 + 92 + ^gWi + 2<7o0i) + • • • 2 = h {CN, + ft + +T H + ft + 2 + l{cN 2 + %- (ft)' - ( * ) " f t ft+ [ ( f t ) ' " ft] rf'>} 3 2 ft^)} + ••• Substituting Eq. (C.16) into Eq. (C.7) and the result into the boundary conditions depicted in Eq. (C.6) yields, 0o(*i) + «Sfi(*i) + e ? ($i) + ° ( 2 2 g ($ ) + e (9 ) 0 gi 2 £ 3 ) = 0 + e g ($ ) + 0(e ) = 0 2 2 (CM) 3 2 2 Assuming the wind pressure distribution on the cylindrical membrane obeys an arbitrary cosine law, we have, C((f>) = cos(or^) (C.25) where a is an arbitrary parameter. 0(e°) Solution From Eqs. (C.19), (C.23) and (C.24) we have the following results, 9o + 9o=' (cNi + ft) + | cos(c^) 2 (C.26) where N = Rp 0 a and Sb(*i) = 0o(* ) = 0 2 (C.27) Perturbation Analysis / C 151 Solution to Eq. (C.26) is given by, g {<p) = A cos <f> + B sin <f> + ^{C + ^) Q 0 cos(c^) J_ Q (C.28) where \B S Q 2sin($ -* ) l/ (S )cos$ -/ ($ )cos$ 2 1 l 1 2 l 2 l J l " J with Observe that Ni is unknown at this stage and will be determined in the 0[e ) l solution. 0(e ) Solution 1 The governing equations from Eqs. (C.19), (C.23) and (C.24) are as follows, + 9i = J (CW + 2 ft " f t ) ~ ^ft ° ( )2 C S ( ^ } " ^ " ^° ( - °) 2 C 3 where and <7i(*i) = <7i(* )=0 (C.32) 2 Substituting Eq. (C.28) into Eq. (C.31) yields the first order approximation to the hoop stress, namely, _ 2 Integrating yields "i — a f c o s a ^ + cosa$ }(l — cos(A$)) — {sina$ p 2 — 2 sina$ }sin(A$) x v- (a - l)a{CA$sin(A$) + ( C + ^)[2 - A$sin(A$) - 2cos(A$)]J 2 where (C33) Perturbation Analysis / C 152 Solution to Eq. (C.30) is given by, gi(<t>) = Ai cos<f> + Bi sin <f> + a cos(2a<^) + a cos(a + 1)^ + a cos (a — 1)<£ x 2 3 + a sin(a + l)<f> + a$ sin(a — 1)^ + a ^ s i n ^ + ai<f>cos<f> + a% cosa<£ + a 4 9 (C.34) where - l °1 a 8(a'^-l)(4o2-l) _ 3 = °2 - °4 2(a-2)to-l)a Bo ° S _ a = 2(a-2)[a-l)a 6 a = iJV B (C+^) 7 1 - °8 - 0 The integration constants A Bi u 2a(a+l)(a+2) 2a(o+l)(a+2) -i^ (C 0 + 3 r 2 0 (C.35) 2(a2-l)2 2A (a -l) r }-) in Eq. (C.34) are determined from the boundary condi- tions shown in Eq. (C.32) and this gives, / ( * ) a n * i ~/2(^i)sin$ t / ^ ^L ) C O S $ - / 2 ( 2 ) l 2 {Bi J sin(* -*i) 2 2 (C.36) 2 $ c 2 o s $ where / (^) = ox cos(2or0) + a cos(a + l)<f> + o cos(a — l)<f> + o sin(a + l)<f> 2 2 3 4 + a sin(a — 1)^ + a 4> sin ^ + 07 <f> cos ^ + ag cos(ar<£) + a 5 Again N 0(s ) 2 6 (C37) 9 is unknown and will be determined in the 0(e ) solution. 2 2 Solution Since it is intended that the results of the perturbation analysis be evaluated up to 0(e ) l accuracy, the governing equations for this section will not be produced here. The purpose of this section is to determine the second order approximation to the hoop stress N . This 2 information is necessary for the complete solution of the O^e ) problem. Hence, we have 1 from Eq. (C.19), Perturbation Analysis / C 153 N Substituting Eqs. * C{* -* )I 'o i9l+9 2 2 = 2 l (C.28) and (C.34) into Eq. + 92 ° )d4, ( a 3 8 ) (C.38) we have, Vf 2 •^2 = "FHTk *~T / S - ^ l c o s ^ + B\ $ + i cos(2a^) + a cos(a + 1)# + a cos(a - sm 1)<£ a 2 3 + a sin(a + l)<f> + o sin(a - 1)<^ + a ^ sin <f> + a #cos<f> + a cos(a^) + a 4 5 + 42 , 2 , 0L( . i\ + 2(^1) K - 0 a , 2 r B 6 c o s <* + 2 ( + W ~( + ) a a l 7 8 9 cog(2afl 1 C 0 8 ( - *M a + 2 ( ^ 1 ) K " ) ( + )+ + ( + !) (« " a J s i n a l a sin ( C + jfe) cos^ + iV^o + iV^o ( c + jfe) sin* - ^ t ^ f cos(a^) J <ty Integrating yields, 2/t (sin $ - sin $x) + 2fe (cos $ ~ t 2 ~ CA$ sin( A$) 2 2 c o s 2 + (C + ^){2 - A $ sin( A$) $ 1 ) + 2fe sin(A$) (C39) 3 - 2 sin $ sin $ - 2 cos $j cos $ } x 2 ' 2 in which A$ = $ - $1 2 ^1 = /3(^2) sin - / ( $ i ) sin $ 3 2 ^2 = /3(^2) cos $! - / ( $ i ) cos $ 3 2 / i = 6 (sin 2a$ — sin 2a$ ) + 6 ( 1 3 2 x + 6 {sin(a + 1)$ 3 _ 2 s m a 2 $2 ~ s m sin(a + l)$i} + fc {sin(a — 1)$ - sin(a 4 + 6 {cos(a + 1)$2 - cos(a + 2 l)$i} + 6e{cos(a - 1)$ - cos(a - l ) $ i } + b 5 2 7 (C.40) where /(</>) and the constants bi... bj are defined as follows, 3 h{4>) = i cos(2a^) + a cos(a + 1)* + a cos(a — 1)# + 04 sin(ar + 1)* + 05 sin(a — l)<f> a 2 3 + a <f> sin <f> + a <f> cos<j> + a cos(or^) - | (ft) 6 7 ^ + B ?+ £ ( c+ 8 A ) J + i ^±L] Perturbation Analysis / 154 C and *1 = 63 —a 4(a -l)(4o a Ao = ; [, 2(a+l)» V h= Bo f. 2(a+l)» V '. = ^ { ' + 1^} b = 7 - 5 (COS$ 0 2 C03$! - $ 2 8 m $ 2 + ^1 Sin - g(^) A$ 2 <- > C 42 APPENDIX D A Globally Convergent Solution Algorithm for Newton-Raphson Iterations Introduction It is a well known fact that Newton-Raphson iteration scheme converges q-quadratically if it is given a sufficiently good starting guess. However, the method may not converge at all if the guess is poor. Hence a globally convergent solution algorithm incorporated 1 into the Newton-Raphson iteration scheme would be highly desirable. Several globally convergent methods with varying degrees of sophistication are available (see Ortega and Rheinbolt (1970) and Dennis and Schnabel (1983)). However, Allgower and Georg (1980) has pointed out that any method that guarantees convergence from any starting point is probably too inefficient for general use. The globally convergent algorithm presented here is the method of line search. The basic concept is simple: we move along a given descent direction, d k by such an amount so as to yield an acceptable next iterate. That is, at iteration k: compute Pi. > 0 such that becomes an acceptable next iterate. By a globally convergent algorithm, we mean a scheme that is designed to converge to some solution from any starting point Mhnost 155 A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 156 The line search parameter, fi in Eq. (D.l) exhibits these features in the solution iteration: k ' < 1 decelerating N-R iteration = 1 standard N-R iteration > 1 accelerating N-R iteration Pk (£-2) k In well-behaved problems, (3 varies with the iteration cycle in this manner: usually, at k the starting guess, {3 < 1. As the solution approaches convergence, f3 —• 1 and when the k k solution becomes very close to convergence,fi > 1. The implementation of the acceleration k feature for general cases can be complicated and will not be attempted here. Following the developments given in Dennis and Schnabel (1983), we shall discuss next the conditions governing global convergence of a solution method. Conditions of Global Convergence For simplicity, the following discussion is restricted to one-dimensional nonlinear problems. This basic framework will then be extended to handle multi-dimensional nonlinear problems. Consider a nonlinear problem defined below, of which solution is desired: F{x)=0 (D.S) Since a globally convergent scheme for an associated minimization problem is easily formulated, we shall convert the problem of solving a nonlinear equation to one of minimizing an associated function. This is achieved by forming the l norm of Eq. (D.3) which is, 2 /(x) = i | F | A new solution x i k+ f(x ). k t {DA) 2 is obtained from the current solution x kt This is achieved by computing x by requiring that /(xjt+i) < in a descent direction d k+l ki from x . Mathek matically, d is a descent direction from x if the directional derivative of / at x in the k k k direction d is negative, namely, k V/(x )4 < 0 fc (D.5) A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 157 If Eq. {D.5) holds, then the following is ensured for sufficiently small fi > 0: k f{x +/3 d )<f{x ) k k k {D.6) k It should be pointed out that satisfaction of Eq. {D.6) does not guarantee that the solution will converge. 2 A more refined criteria is required. Armijo (1966) and Goldstein (1967) have proposed that the line search parameter to be selected from among those /3 > 0 k namely, (3 should satisfy these two relationships: k f{x + (3 d ) < f{x ) + a(3 Vf{x )d k k k k k k {D.l) k and V / ( x ) t i K « V/Cxjfe + P d )d > XVf{x )d + k k k k {D.S) k where a G (0,1) and A G (a, 1) are constants. Eqs. {D.l) and {D.S) are satisfied simultaneously if A > a. Note that Eq. {D.l) ensures that the decrease in f{x) relative to the step length is not too small while Eq. {D.S) ensures that sufficiently large steps are used. A plot of these two line search conditions is given in Figure 1. Having presented the two line search conditions, we can now state the following powerful result: that there exist j3 > 0 satisfying Eqs. {D.l) and {D.S) for any k given descent direction d Vf( k)( k x x - k+i) < 0. x a Q k and that any algorithm that generates solution x obeying k d Eqs. {D.l) and {D.S) at each iteration is essentially globally convergent. A systematic procedure for computing the line search parameter fi is given 3 k next. Computation of the Line Search Parameter, (3 k Any hybrid solution method that combines a globally convergent strategy with the fast local convergence of the Newton-Raphson iteration should be so designed that the full Newton step is attempted first. That is, f3 = 1 is first invoked. This is done to ensure k 2 8 See Dennis and Schnabel (1983), p. 117 See Dennis and Schnabel (1983), p. 120-125 A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 158 y Permissible Values of {5 k FIGURE D.l The Two Line Search Conditions. that rapid convergence rates of the Newton-Raphson iteration are fully exploited should the starting guess be close to the solution. If (J = 1 produces an unacceptable next iterate, k then a systematic backtracking of fa will be called upon. The backtracking strategy presented here avoids excessively small steps and the line search condition in Eq. (D.S) is thus, generally not required. Hence only the line search condition in Eq. (D.l) remains to be satisfied. Dennis and Schnabel (1983) has recommended that a = 10~ in Eq. (D.l). 4 This leaves only the f3 to be determined. Defining, k f(P) = f(x k + {3d ) k A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 159 we thus have the following known information: /(0) = /(**) and /'(0) = Vf(x )d k (D.10) k Since p = 1 is used at the start of an iteration, the following information is also known: k f(l) = f(x + d ) k If P = 1 is not acceptable, namely f(x +d ) k k k (0.11) k does not satisfy Eq. (D.7)* then a quadratic curve is drawn through the three known points given by Eqs. (D.10) and (D.ll). The P that minimizes this curve is then selected as P . This is displayed in Figure 2. Thus we k have for a general quadratic curve, y{P) = ap + bp + c 2 where constants a, 6, c are evaluated by fitting the curve through the three points resulting in, y(P) = [/(l) - /(0) - f'(0)]P + />)/? + /(0) (£U2) 2 Minimizing Eq. (D.12) yields Pmin = — ~> ^ 2[/(l)-/(0)-/'(0)] (0.13) Note that in Eq. (D.13), since /'(0) < 0 and that / ( l ) > /(0) + a/'(0) > /(0) + /'(0) we A thus have P min A A > 0. Also observe that since / ( l ) > /(0) + af (0) we get, 1 If / ( l ) 3> /(0), then p min < ^. This results in unnecessarily small P min and is not A desirable since it implies that f(P) is poorly modeled by a quadratic in this region. Dennis and Schnabel (1983) has recommended fixing a lower bound of /? ,„ > pj at the first m backtrack during an iteration. Thus the line search parameter at the first backtrack is given by, 1 > Pk = Pmin > TS * Le. /(I) > /(0) + a/'(0) ( D U ) A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 160 A (0,/(0)) slope = /(0) y = [/(l) - /(0) - /'(0)]/3 2 + / W + /(0) P P= l 'min Quadratic Backtrack. F I G U R E D.2 Suppose after the quadratic backtracking, the line search condition in Eq. {D.l) is still not satisfied. We now have four pieces of known information, the fourth being, f{P ) = f{x + P d ) k k k {D.15) k The more accurate cubic curve will now be drawn through these four points and is given by, y{p) = ap + bp + f'(0)(3 + /(0) z where fa\_ {D.16) 2 1 \l! -1 M ) - /(0) - firPW 1 in which /? , fi are two previous values of P . They need not necessarily be 1 and /3 ,„ r t k m obtained in the quadratic backtracking as the cubic backtracking will be applied repeatedly A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 161 to reduce fa until the line search condition in Eq. (0.7) is satisfied. Again, minimizing Eq. (0.16) yields, -6 + \/6 -3a/'(0) 2 Pmin = — (0-17) It can be shown that if a < | , Eq. (0.17) is never complex. Figure 3 shows the cubic backtrack model. As in the quadratic backtrack model, the following bounds recommended by Dennis and Schnabel (1983) applies: \Pr >Pk = Pmin > T5 A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 162 Application to Multi-Dimensional Nonlinear Problems By the following simple extension, the line search concept developed for one-dimensional problem will now be modified to handle problems with several variables. Consider the following system of nonlinear equations: F(x) = 0 (D.19) Solving Eq. (D.19) by Newton-Raphson iteration yields at the k x k+l where J(x ) k =x - t h iteration: J(x )~ F(x ) (D.20) l k k k is the Jacobian of F at x^. Again using the l norm of F(x) as a measure 2 of convergence we have, (D.21) /(x) = i F ( x ) F ( x ) r We have now converted the problem of solving a system of nonlinear equations to one of minimization. The descent direction from x * for the minimization problem depicted by Eq. (D.21) that is, dk satisfies, V/(x*) d r k < 0 We shall now show that the Newton direction at xk, —J(xk)~ F(xk) is a descent direction 1 for the minimization problem: Vf(x ) d = -F (x )J(x )J(x )- F(x ) T k T k 1 k k k k = -F (x )F(x ) T k k <0 (0.22) provided that F(x ) ^ 0. An alternative positive definite quadratic model can also be k developed for the system of nonlinear equations but this will not be pursued here. Hence 5 we see that the global technique for the one-dimensional problem can be applied to multidimensional problem if we define the objective function using the l norm and the descent 2 direction by the Newton's direction. 5 See Dennis and Schnabel (1983), p. 148-149 A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 163 The global algorithm presented so far, do suffer from failures and two notable examples are, i) that J(x ) k is singular. We have to ensure that J(x ) k conditioned resulting in J{x ) J(x ) T k k is sufficiently well- being positive definite. This is not a problem if the structure considered is properly restrained and stable. ii) that a local minimum of f[x) = \F(x) F(x) is not a root of F(x) = 0. While it is true that every root of Eq. (0.19) is a root of Eq. (0.21), the converse is not necessarily true. If the starting guess is close to one such local minimum, then the global strategy might converge to this point. There is nothing much one can do other than to restart from a new guess and hope for better luck this time. APPENDIX E Undeformed Roof of BC Place Stadium The data-set corresponding to the undeformed roof at B C Place Stadium is: Node N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 X 91 . 440 8 2 . 110 8 7 . 780 9 0 . 530 74. 520 6 8 . S80 6 8 . 580 6 8 . 580 6 8 . 580 6 8 . 580 54. 680 45. 720 45. 720 45. 720 45. 720 45. ,720 28. ,850 22. .860 22. .860 22. .860 22. .860 22. .860 0..0 0,.0 0 .0 0 .0 0 .0 0 .0 -28 .850 -22 .860 -22 .860 -22 .860 -22 .860 -22 .860 -54 .680 -45 .720 -45 .720 -45 .720 -45 .720 -45 .720 -74 .520 -68 .580 -68 .580 -68 .580 -68 .580 -68 .580 -82 . 110 -87 .780 - 9 0 .530 -91 .440 Coordinates Y 0.0 67.300 44.870 22.430 82.860 80.760 67.300 44.870 22.430 0.0 100.580 82.260 67.300 44.870 22.430 0.0 110.450 89.740 67.300 44.870 22.430 O.O 112.170 89.740 67.300 44.870 22.430 0.0 110.450 89.740 67.300 44.870 22.430 0.0 1OO.580 82.260 67.300 44.870 22.430 0.0 82.860 80.760 67.300 44.870 22.430 0.0 67.300 44.870 22.430 0.0 2 0.0 0.0 0.0 0.0 0.0 4 .410 7.010 12.250 13.430 13.850 0.0 10.550 16.810 19.070 20.900 21 .550 0.0 11.980 18.720 22.560 24.730 25.490 0.0 11.140 18.540 23.750 26.030 26.830 0.0 11.980 18.720 22.560 24.730 25.490 0.0 10.550 16.810 19.070 20.900 21 .550 0.0 4.410 7.010 12.250 13.430 13.850 0.0 0.0 0.0 O.O Note that the Z-coordinates are measured from the base of the roof. 164 PUBLICATIONS M.D. Olson and P.S. Han, "Interactive Analysis of Wind Loaded Membranes", 2nd Int. Conf. On Numerical Methods for Nonlinear Problems, Barcelona, Spain, April 1984. R.L. Johnston, G. Fairweather and P.S. Han, "The Method of Fundamental Solutions, An Adaptive Boundary Element Method for Problems in Potential Flow and Solid Mechanics," 5th ASCE Specialty Conf., Eng. Mech. Div., Laramie, Wyoming, August 1984. P.H. Han and P. Lukkunaprasit, " F i n i t e Strip Analysis of Framed Tube Structures," 3rd Int. Conf. on Tall Buildings, Hong Kong and Guangzhou, China, December 1984. P.S. Han, M.D. Olson and R.L. Johnston, "A Galerkin Boundary Element Formulation With Moving S i n g u l a r i t i e s , " J . Eng. Comput., Vol. 1, 1984. P.S. Han and M.D. Olson, "Static Analysis of Pneumatic Structures Loaded by Wind," IASS Int. Symp. on Membrane Structures and Space Frames, Osaka, Japan, September 1986.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Interactive nonlinear analysis of wind-loaded pneumatic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Interactive nonlinear analysis of wind-loaded pneumatic membrane structures Han, Peng Siew 1986
pdf
Page Metadata
Item Metadata
Title | Interactive nonlinear analysis of wind-loaded pneumatic membrane structures |
Creator |
Han, Peng Siew |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | A nonlinear analysis of pneumatic membrane structures loaded by wind, taking into account the fluid-structure interactions, is presented here. The analysis comprises two parts: the flow part which is solved using a boundary element procedure and the structures part which is analysed with a finite element approach. These two numerical techniques are then coupled to capture the effects of the fluid-structure interaction in the analysis. The wind is represented by potential flow and thus Laplace's equation governs. The boundary element method employed in the analysis, features moving singularities of the fundamental solutions. The singularities are placed on an auxiliary boundary located outside the domain of interest. Results from solved examples show the technique to be highly accurate for the computation of wind loads. The method is however, nonlinear with the nonlinearity occuring in the determination of the auxiliary boundary. The membrane, assumed to be perfectly flexible and to exhibit only exten-sional rigidity, is modeled by a small strain-finite deformation theory. A Galerkin finite element approach is formulated for the analysis of two-dimensional pneumatic membrane structures. The analysis is based on an infinitely long cylindrical membrane with an initially circular cross-section. For three-dimensional pneumatic membrane structures, the analysis is accomplished with a continuum mechanics-based incremental finite element formulation using a convected coordinate description. The metric tensor for any deformed configuration is derived and from this, a consistent measure of deformation is obtained. In both analyses, the governing finite element equations are solved via Newton-Raphson iteration with an 'automatic' line-search scheme. A pertubation technique is developed to check the accuracy of the proposed two-dimensional finite element analysis for low wind speeds. At higher wind speeds, comparison with available solutions are made. To assess the performance of the proposed three-dimensional finite element analysis, a square membrane and a cylindrical membrane are analysed. Good agreement with published results is observed. Having established the accuracies of the boundary element and the finite element methods on independent flow and structures problems respectively, interactive analysis is then carried out. Two pneumatic membrane structures are analysed: an infinitely long cylindrical membrane as a two-dimensional example, and the air-supported roof of the BC Place Stadium as a three-dimensional example. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-06 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062907 |
URI | http://hdl.handle.net/2429/27105 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1986_A1 H36.pdf [ 9.54MB ]
- Metadata
- JSON: 831-1.0062907.json
- JSON-LD: 831-1.0062907-ld.json
- RDF/XML (Pretty): 831-1.0062907-rdf.xml
- RDF/JSON: 831-1.0062907-rdf.json
- Turtle: 831-1.0062907-turtle.txt
- N-Triples: 831-1.0062907-rdf-ntriples.txt
- Original Record: 831-1.0062907-source.json
- Full Text
- 831-1.0062907-fulltext.txt
- Citation
- 831-1.0062907.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062907/manifest