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 ac-count 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 gov-erns. The boundary element method employed in the analysis, features moving singularities of the fundamental solutions. The singularities are placed on an auxiliary boundary lo-cated 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 ini-tially circular cross-section. For three-dimensional pneumatic membrane structures, the analysis is accomplished with a continuum mechanics-based incremental finite element for-mulation 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, com-parison 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 in-finitely 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 Prologue 1 1.2 Literature Review 3 1.3 Scope of Investigation and Thesis Layout 5 1.4 Nomenclature 8 2 Two and Three Dimensional Boundary Element Analyses 10 2.1 Introduction 10 2.2 Concept of Moving Singularities 11 2.3 Theoretical Formulation 12 2.4 Computer Implementation 15 2.5 Numerical Results 16 2.6 Summary 30 3 Two Dimensional Finite Element Analysis 33 3.1 Introduction 33 3.2 Assumptions and Theoretical Formulation 36 3.3 Galerkin Finite Element Formulation 42 3.4 Evaluation of Finite Element Matrices 45 3.5 Computer Implementation 48 3.6 Numerical Results 51 3.7 Summary . 67 Contents iv 4 Three Dimensional Finite Element Analysis . . . . 68 4.1 Introduction 68 4.2 Mathematical Preliminaries 69 4.3 Theoretical Formulation 74 4.4 Isoparametric Finite Element Solutions 84 4.5 Computer Implementation 94 4.6 Numerical Results 100 4.7 Summary 128 5 Summary and Conclusions 130 5.1 The Flow Analysis 131 5.2 The Structural Analysis 132 5.3 The Interactive Analysis 133 5.4 Recommendations for Further Study 134 References 135 Appendix A 140 Appendix B 144 Appendix C 146 Appendix D 155 Appendix E 164 Tables V 2.1 Location of singularities and the associated coefficients for rectan-gular membrane 19 2.2 Comparison of membrane deflection: Exact vs. B.E.M 20 2.3 Location of singularities and the associated coefficients for flow past rigid circular cylinder 23 2.4 Location of singularities and the associated coefficients for heat conduction in a rectangular prism 25 2.5 Comparison of temperature and its gradient for heat conduction in a rectangular prism 25 2.6 Location of singularities and the associated coefficients for flow past 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 . 1 1 2 4.7 Stress resultants at the apex of the roof, for the inflation problem . 1 1 4 4.8 Stress resultants at the apex of the rigid roof (V = 25 m/s) . 1 1 9 4.9 Maximum roof elevations (initial height=28.81m) 120 4.10 Stress resultants at the apex of the flexible roof 125 4.11 Maximum shear stress resultants and their locations 128 vi Figures 1.1 Wind loads on a pneumatic membrane 1 2.1 Problem definition and auxiliary boundary 13 2.2 Rectangular membrane 17 2.3 Discretization scheme and initial singularity locations for rectan-gular membrane, • Boundary Node, x Singularity 18 2.4 Loci of singularity movement in the direction of arrows, for rect-angular membrane (domain of problem shown shaded) 19 2.5 Discretization scheme, initial singularity locations, flow domain and boundary conditions for cylinder, • Boundary Node, x Singularity 21 2.6 Loci of singularity movement in the direction of arrows, for flow past circular cylinder (domain of problem shown shaded) 22 2.7 Pressure coefficient, Cp on a rigid circular cylinder 23 2.8 Discretization scheme and initial singularity locations for rectan-gular prism, • Boundary Node, x Singularity 24 2.9 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 for flow past rigid sphere. 27 2.11 Boundary elements for sphere and the flow domain 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 Figures vii 3.7 Comparison of membrane deflection 56 3.8 Membrane deformations for Example 1 (£„ = 1.91, C = \ x 10 - 4) . . . 59 3.9 Membrane deformations for Example 2 (<fa = 3.40, C = ± x IO - 4) . . . 59 3.10 Flow domain and boundary conditions for cylindrical membrane . . . 61 3.11 Membrane deformations at each global iteration (£ a = 0.276, C = £ x l 0 " 4 ) 63 3.12 Membrane deformations at each global iteration (£ 0 = 1.875, C = \ x IO"4) 63 3.13 Pressure loads acting on membrane in equilibrium (C = \ x 10 - 4) . . . 64 3.14 Influence of the loading ratio, <f0 on membrane deflections (C = J x 10~4) 65 3.15 Variation of membrane stress resultants Na, with internal pressure, pa (C = i x IO"4) . . 66 4.1 Coordinate system for a surface in three dimensional space . . . . 70 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 incre-ments: square membrane 103 4.10 Variation of central deflection with internal pressures for the square membrane 105 4.11 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 (<fa = 1.91, C = ^xlO - 4 ) . 107 Figures viii 4.14 Membrane deformations for Example 2 (£„ = 3.40, C - | x l O - 4 ) . . . 108 4.15 Stadium and flow domain . 110 4.16 BC Place Stadium, Vancouver, Canada I l l 4.17 Undeformed and initial roof profiles due to internal pressure . . . . 113 4.18 Stress resultants along the minor axis for the inflation analysis . . 114 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 the flow model investigations 119 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) 121 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) 122 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) 123 4.26 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 . 161 E. l Finite element grid of BC Place Stadium Roof 165 ix Acknowledgement The author wishes to express his immense gratitude and thanks to his ad-visor, 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 boun-dary 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 graph-ics and Mr. Fon Chow for providing free drafting of figures, too difficult to produce on a computer, are also much appreciated. The financial support 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. C H A P T E R 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 tension1 with the structural efficiency of a shell. Hence, pneumatic membrane structures 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 determina-tion 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 pneu-matic 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 avail-able 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 3 1.2 L I T E R A T U R E R E V I E W The British engineer, Frederick Lanchester (1917) appears to have been the first to intro-duce 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 struc-tures, Lanchester (1938) mentioned the potential of pneumatic membrane structures for enclosing large areas like those in aircraft hangers and sports arenas. Unfortunately, Lanch-ester 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 pro-posals (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 Sym-posium 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). Introduct ion / 1.2 4 With the advent of high-speed computers, numerical methods based on fi-nite element modelling became the main tools of analysis. Becker (1966) solved the prob-lem 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 mem-brane 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 Gaus-sian 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 mem-brane 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 struc-tures, very few published results are available. In an attempt to overcome this deficiency, Kunieda, Yokoyama and Arakawa (1981) presented some solutions obtained from an in-teractive 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.1 5 1.3 S C O P E O F I N V E S T I G A T I O N A N D T H E S I S 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 the fluid-structure inter-actions. 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 contin-ued 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 — Three-Dimensional Finite Element Analysis e) Chapter 5 — Summary and Conclusions 1.3.1 T H E F L O W A N A L Y S I S 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 auxil-iary 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 conven-tional 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 sig-nificantly 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 mem-brane, 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 A N A L Y S I S 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 three-dimensional pneumatic membrane structures. In Chapter 3, the solution of two-dimensional pneumatic membrane struc-tures 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 con-sidered. 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 mem-brane structures. Solution to the structural deformations is obtained through a contiuum mechanics-based incremental finite element formulation using a convected coordinate de-scription. 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 stiff-ness 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 non-orthogonal 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 A N A L Y S I S By coupling the two numerical methods in the manner discussed earlier, we generated some interactive results. All results pertaining to the two-dimensional pneumatic mem-brane 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 mem-brane is considered in the former while the air-supported roof of the BC 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 ap-pear. In the two dimensional formulation where a rectangular cartesian coordinate system is adopted, the notations will closely follow those given in Flugge (1973). Thus Nx, N^, Nx<f, denote membrane stress resultants; cxi s^ , 7X(^ for membrane strains and so on. In the three 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 in-dices, 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 differ-ent 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: 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 This quantity occurs in the undeformed configuration 0. The indices i,j j 1 as shown denote the quantity is a mixed second order tensor. dr r aad6a = aldOx + a2 d92 & = * ! Si + J Si + S3 M 7 £ / =A f 1 e 1 + M 2 6 + M 3 6 + - + A f 1 2 £12 NN y\ = N1 y\ + N2y^ + N3 yj + NA y\ 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. C H A P T E R 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 accu-rate 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 C O N C E P T OF MOV ING S INGULARIT IES 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 singu-larities 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 / 2.3 12 Mathon and Johnston (1977) were the first to introduce this concept of mov-ing 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 approx-imate 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) (2.3) a(x) + /?(x)u(x) + 7(x)c5u(x) fdn (Mixed) 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 / 2.3 13 ^ A u x i l i a r y Boundory F I G U R E 2.1 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 = £Cyf(« y,x) t^D,T (2.4) j'=l 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 / 2.3 14 boundary which is denned by, SRB = Bu ± 0 (2.5) For convenience, let Zk = (Ck, tk)T be the vector of unkowns. Then discretizing the boundary with M boundary elements and minimizing the boundary residue fkB yields, A/ ]T J[wkxB]idr = o or M JT f[wkBu}idT = 0 l<k<N (2.6) ,=i ff where T,- is the boundary for the Ith element, and the weighting functions wk are given by, (d(Bu)/dZk for least square minimization k \ du/dZk for Galerkin minimization where d dZk \dCk' dtk/ 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.4 15 2.3.1 S P E C I A L I Z A T I O N T O L A P L A C E ' S 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, in 2-dimensions: c(t-,x) = —In (\=^—r— U j U + W (2.9) in /-dimensions: ${tj,x) = |ty - x|2-' - (R + |t y|) 2"', / > 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 boun-dary 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 M Net— (for 2-D) and 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 results1 is likely, this practice is not 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 R E S U L T S 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 mem-brane 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 / 2.5.1 17 y F I G U R E 2.2 RecfcanguJar Membrane. lateral pressure, as depicted in Figure 2.2. The differential equation and boundary condi-tions for the deflection z are, d2z <Pz = _ dx2 + dy2~ q ( 2 n ) z = 0 on x = ±a; 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, « = * + Uix2 + y2) The differential equation and boundary conditions in terms of u are then given by, 3 2u d2u (2.12) + dx2 dy2 = 0 (2.13) « = \q(x2 + y2) on x = ±a; y = ±b Two and Three Dimensional Boundary Element Analyses / 2.5.1 18 F I G U R E 2.3 Discretization Scheme and Initial Singularity Locations for Rectangu-lar Membrane, • Boundary Node, x Singularity. The solution to Eq. (2.13) is, 16a2g ^ (_i)! n=l,3,5., cosh(n7r«/2a)l /nirx\ , . ? ,x . 1 r?—T7TTT cos ( — - J + \q(x2 + y2) (2.14) cosh(nrr6/2a)J V 2a / 4 ^ 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. All three showed little motion as the solution becomes nearly converged. Two and Three Dimensional Boundary Element Analyses / 2.5.1 19 o X-Axis F I G U R E 2.4 Loci of Singularity Movement in the Direction of Arrows, for Rectan-gular Membrane (Domain of Problem shown Shaded). I n i t i a l V a lues F i n a l V a l u e s X y Coefficient X y Coefficient 1.10 0.10 -0.0045 0.08 3.26 0.1552 1.15 0.65 -0.1145 1.21 0.69 -0.1760 0.40 0.60 0.0473 0.27 0.89 0.0380 Table 2.1 Location of Singularities and the Associated Coefficients for Rectangular Membrane Two and Three Dimensional Boundary Element Analyses / 2.5.2 20 Location Exact 2-D BEM Solution 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 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 Table 2.2 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 F L O W P A S T A R I G I D C I R C U L A R C Y L I N D E R 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 approxi-mated by a finite domain as depicted in Figure 2.5. Symmetry about the x-axis is invoked. The discretization scheme, initial locations of the singularities and the boundary condi-tions 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 X J X k y X X X \ du_ V 2 ^ dn » / x f X dn , u = 0 du X dn V ^ - R = I.O ' A0 x\ » IC X .0 -« 10.0 * X F I G U R E 2.5 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 Fig-ure 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 solution2 is shown in Figure 2.7. The exact pressure coefficient is given by, Cp{9) = l - 4 s i n 2 0 (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 func-tion, we see that the derivatives on the boundary are very accurately computed. The exact pressure coefficient is computed assuming an infinite flow domain. Since the diameter of our cylinder is very small compared to the dimension of our finite domain, we can justifiably make the comparison. Two and Three Dimensional Boundary Element Analyses / 2.5.2 22 FIGURE 2.6 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 23 Initial Values Final Values X y Coefficient X y Coefficient 11.00 3.00 -0.7080 14.43 3.06 -0.7507 11.00 11.00 -1.3351 14.28 9.89 -1.1937 8.00 11.00 -0.0624 8.73 14.40 -0.6406 -8.00 11.00 0.0624 -8.73 14.40 0.6406 -11.00 11.00 1.3351 -14.28 9.89 1.1937 -11.00 3.00 0.7080 14.43 3.06 0.7507 -0.70 0.20 0.0840 -0.27 0.49 0.0298 -0.40 0.40 0.1677 -0.32 0.14 0.3598 0.40 0.40 -0.1677 0.32 0.14 -0.3598 0.70 0.20 -0.0840 0.27 0.49 -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, Cp on a Rigid Circular Cylinder. Two and Three Dimensional Boundary Element Analyses / 2.5.3 A y 24 B3 B2 -T = 3 0 0 B4 dT V 2 T = 0 d n = 0 B l T = 0-B 6 - -12 13 - i — 3.0 -3.0 F I G U R E 2.8 Discretization Scheme and Initial Singularity Locations for Rectangu-lar Prism, • Boundary Node, x Singularity. 2.5.3 H E A T C O N D U C T I O N I N A R E C T A N G U L A R P R I S M Heat conduction in a rectangular prism is used as an example of a prescribed mixed boun-dary 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 - 3 . From all three plots of singularity motion, we see 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 as the solution approaches convergence. 25 Initial Values Final Values X y Coefficient X y Coefficient 3.20 1.00 0.0198 10.85 3.34 0.0560 3.20 3.20 0.0093 2.92 x 108 5.32 x 107 1.99 x 106 1.00 3.20 -0.0120 1.44 10.30 -0.0702 -1.00 3.20 -0.0236 -2.99 8.71 -0.0676 -3.20 3.20 -0.1030 -7.77 6.93 -0.2107 -3.20 1.00 -0.0599 -8.90 2.07 -0.1296 Table 2.4 Location of Singularities and the Associated Coefficients for Heat Con duction in a Rectangular Prism Location Temperature Temperature Gradient Exact 2-D BEM Exact 2-D BEM Boundary Points Bl(-3.0,0.75) 300.0 300.0001 -50.0 -50.0003 B2(-3.0,2.25) 300.0 299.9995 -50.0 -49.9986 B3(-1.0,3.00) 200.0 200.0002 -50.0 -50.0004 B4( 1.0,3.00) 100.0 100.0003 -50.0 -49.9995 B5( 3.0,2.25) 0.0 0.0003 -50.0 -49.9999 B6( 3.0,0.75) 0.0 0.0001 -50.0 -50.0001 Interior Points Il(-2.0,0.0) 250.0 250.0000 -50.0 -50.0003 I2( 0.0,0.0) 150.0 150.0000 -50.0 -49.9998 I3( 2.0,0.0) 50.0 50.0001 -50.0 -50.0002 Table 2.5 Comparison of Temperature and its Gradient for Heat Conduction in a Rectangular Prism T h e results for 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 gradient c o m p u t e d at vario u s boundary a n d interna] p o i n t s are s u m m a r i z e d i n T a b i e 2.5 along w i t h the exact solutions. 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 27 2.5.4 F L O W P A S T A R I G I D 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 solution3 derived for an infinite domain, is given by, Cp{9) = l-\sm20 (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 solu-tion 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 30 Initial Values Final Values X y z Coefficient X y z Coefficient 0.50 0.00 0.00 -0.0019 0.38 0.00 0.00 - 0.0003 0.42 0.27 0.07 0.0030 0.26 0.20 0.09 0.0056 0.36 0.30 0.18 -0.0098 0.12 0.07 0.06 - 0.5420 0.28 0.33 0.26 0.0046 0.26 0.24 0.25 0.0004 0.19 0.20 0.42 -0.0031 0.13 0.06 0.29 - 0.0037 -0.19 0.20 0.42 0.0031 -0.13 0.06 0.29 0.1760 -0.28 0.32 0.26 -0.0046 -0.26 0.24 0.25 -0.0004 -0.36 0.30 0.18 0.0098 -0.12 0.07 0.06 0.5420 -0.42 0.27 0.07 -0.0030 -0.26 0.20 0.09 - 0.0056 -0.50 0.00 0.00 0.0019 -0.38 0.00 0.00 0.0003 8.50 3.20 1.60 -0.3339 24.35 6.03 4.43 -29.8991 8.50 3.20 6.40 -0.3163 19.12 6.19 15.37 -31.7628 4.80 4.80 8.50 0.0284 10.83 16.13 15.23 -17.2196 3.20 1.60 8.50 0.0224 7.55 4.36 18.64 - 4.2419 -3.20 1.60 8.50 -0.0224 -7.55 4.36 18.64 4.2419 -4.80 4.80 8.50 -0.0284 -10.83 16.13 15.23 17.2196 -8.50 3.20 6.40 0.3163 -19.12 6.19 15.37 31.7628 -8.50 3.20 1.60 0.3339 -24.35 6.03 4.43 29.8991 -4.80 8.50 1.60 -0.0380 -16.85 16.32 5.45 21.4340 -3.20 8.50 4.80 -0.0196 -6.00 17.19 3.52 1.7581 3.20 8.50 4.80 0.0196 6.00 17.19 3.52 - 1.7581 4.80 8.50 1.60 0.0380 16.85 16.32 5.45 -21.4340 Table 2.6 Location of Singularities and the Associated Coefficients for Flow Past Rigid Sphere 2.6 S U M M A R Y The distinctive feature of the boundary element method presented here is that the singu-larities 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 / 2.6 31 F I G U R E 2.12 Pressure Coefficient on a Rigid Sphere. Two and Three Dimensional Boundary Element Analyses / 2.6 32 is not seen as a m a j o r drawback. U n l i k e c o n v e n t i o n a l b o u n d a r y element methods w h i c h require the same number of s i n g u l a r i t i e s a n d b o u n d a r y elements, the number of singu-l a r i t i e s i n the proposed f o r m u l a t i o n can be s i g n i f i c a n t l y reduced w i t h o u t any apparent loss of accuracy. T h i s is an i m p o r t a n t advantage i n higher d i m e n s i o n a l problems. F r o m the plots of the l o c i of s i n g u l a r i t y m o t i o n , we see that as the s o l u t i o n progresses towards convergence, the s i n g u l a r i t i e s 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 the d o m a i n of the problem. W h e n the s o l u t i o n is v e r y near convergence, we generally observe m i n i m u m movement of the s i n g u l a r i t i e s . C o m p u t e d values o n the boundary, as w e l l as those i n the int e r i o r , are accurately obtained. I n the c o m p u t a t i o n of de r i v a t i v e s o n t h e boundary, con-v e n t i o n a l b o u n d a r y element techniques can y i e l d very p o o r results w h i c h is not the case for the present method. It s h o u l d be p o i n t e d out t h a t the proposed m e t h o d w o u l d be more expensive 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 element approach because the m e t h o d is not o n l y n o n l i n e a r but also requires e x t r a i t e r a t i o n s i n c o m p u t i n g the o p t i m u m s i n g u l a r i t y locations. One w o u l d therefore have t o balance t h i s 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 curate de r i v a t i v e s on the boundary. T h e proposed m e t h o d w i l l be used i n the flow analysis f o r c o m p u t i n g the fluid l o a d i n g o n p n e u m a t i c membrane str u c t u r e s . C H A P T E R 3 Two Dimensional Finite Element Analysis 3.1 I N T R O D U C T I O N T h e 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 membrane subjected to w i n d loads is presented here. Since p n e u m a t i c membranes e x h i b i t such large deformations i n w i n d , i t is essential that one i n c o r p o r a t e i n t o the analysis, the fluid-structure interactions. I n t h i s thesis, this w i l l be done by c o u p l i n g the finite element p r o g r a m w h i c h pre d i c t s the struc-t u r a l deformations 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 computes the w i n d loads. T o derive the p e r t i n e n t d i f f e r e n t i a l equations f o r the s t r u c t u r a l analysis, a c y l i n d r i c a l pneu-m a t i c membrane s t r u c t u r e whose cross-section is c i r c u l a r i n the undeformed c o n f i g u r a t i o n is considered. B y t a k i n g the c y l i n d e r to be s u f f i c i e n t l y l o n g compared to i t s cross-sectional dimensions, the p r o b l e m becomes two dimensional, thereby p e r m i t t i n g some s i m p l i f i c a -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 deformed i n t o a n o n - c i r c u l a r c y l i n d e r by means of a plane motion. T h e c y l i n d r i c a l membrane held down at the edges, is s u p p o r t e d by i n t e r n a l pressure pa and l o a d e d by e x t e r n a l pressure pw, due t o w i n d b l o w i n g i n the broadside direc-t i o n as shown i n Figure 3.1. G r a v i t y forces are neglected i n comparison w i t h the pressure loadings. A s m a l l s t r a i n - f i n i t e d e f ormation membrane the o r y is derived f o r the s t r u c t u r a l analysis. T h e r e s u l t i n g second-order i n t e g r o d i f f e r e n t i a l equation is solved v i a a G a l e r k i n finite element scheme us i n g Newton-Raphson 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 made 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 mem-branes 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 supris-ingly, 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 ex-pressed 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 the fluid-structure interactions of pneumatic mem-branes 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 defor-mation computation from the solution of the integrodifferential equation. The iterative process is stopped when the structure attains an equilibrium configuration that is consis-tent 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 integrodiffer-ential 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 po-tential 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.1 Also, in the study involving cylindrical pneumatic membranes, Kunieda et. al. presented results that showed at most, a 4% difference be-tween 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 real flow model. But as discussed in Chapter 4, flow separation still occurs in these structures, with re-attachment then taking place at the leeward side Two Dimensional Finite Element Analysis / 3.2.1 36 3.2 A S S U M P T I O N S A N D T H E O R E T I C A L F O R M U L A T I O N Figure 3.1 shows a cylindrical membrane inflated by internal pressure, pa and acted upon by external pressure, pw. Despite the large structural deformations, the internal pressure 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. 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, 3.2.1 E Q U I L I B R I U M E Q U A T I O N S dNx 1 dNx<l> 0 (3.1) dx R d<f> dNx+ 1 dN+ 0 (3.2) dx R d<j> (3.3) where Nx, Nx<^ are the usual membrane stress resultants. Two Dimensional Finite Element Analysis / 3.2.2 37 N I R<M> jc ,U F I G U R E 3.2 Positive Stress Resultants. 3.2.2 S T R A I N - D I S P L A C E M E N T R E L A T I O N S For a t h i n c y l i n d r i c a l s h e l l u n dergoing large deformations, the s t r a i n components exi e^, 7 I (^ are given by, dx dx H[6)2+©2+0] (3-4) £* - 5 (Ij + ") + [(£)2 + - »)2 + S + ")'] + £» <3'5> ^ = Ua? + di) + RTxa* + RTX\M + w) + Rdi {at - v) (3-6> where u, w, tw are displacements i n the x ,y, z d i r e c t i o n s a n d s0 is a p r e s t r a i n t h a t w i l l be determin e d later. Two Dimensional finite Element Analysis / 3.2.4 38 3.2.3 C O N S T I T U T I V E L A W S A s s u m i n g s m a l l s t r a i n s despite the large deformations, t h e c o n s t i t u t i v e r e l a t i o n s f o r a membrane 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 ) where E, G, / i are the elas t i c a nd shear m o d u l i , a n d Poisson's r a t i o s r e s p e c t i v e l y and t is the membrane thickness. 3.2.4 S I M P L I F I C A T I O N S If the c y l i n d e r is su f f i c i e n t l y l o n g compared to i t s cross-sectional dimension, the deforma-tions can be considered to be two dimensional. T h e c y l i n d e r t h u s degenerates i n t o a s t r i n g i n the y — z plane. Hence ex and a l l d erivatives w i t h respect t o x can be set to zero. A l s o , we c a n assume t h a t the r a d i a l displacement component, w is l a r g e r t h a n the circumferen-t i a l displacement, v. T h u s we w i l l ignore the 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 terms i n comparison w i t h the to a n d i t s derivatives. T h e v a l i d i t y of t h i s a p p r o x i m a t i o n w i l l be checked when the membrane is resolved u s i n g the three-dimensional f o r m u l a t i o n where such a p p r o x i m a t i o n s are not intr o d u c e d . T h e results of t h i s c hecking w i l l be 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 equations therefore s i m p l i f y to, 1 dNx<t> _ = 0 (3.10) R d<f> L = 0 (3.11) ldN4>_ R d<f> M 1 + R «" w + F n a)= RiPa ~Pw) (3-12) Two Dimensional Finite Element Analysis / 3.2.4 39 T h e strain-displacements relations reduce to, ex = 0 (3.13) = 5 S + w) + + O'] + e° (3ll4) = 0 (3.15) F r o m Eqs. (3.10) and (3.11) we have, N<t» Nx4> = constants (3.16) and f r o m Eqs. (3.9) and (3.15) we conclude t h a t Nx+ = 0 (3.17) In view of Eq. (3.13) we have 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^ (3.19) where _ /1 - MXMA V EA ) E l i m i n a t i n g f r o m Eqs. (3.14) a n d (3.19) we have, T h e p r e s t r a i n , e0 is evaluated at no w i n d c o n d i t i o n (i.e. v = tw = 0). Thus, we have f r o m Eqs. (3.12) and (3.20), e0 = CRpa (3.21) Eq. (3.20) then becomes, Two Dimensional Finite Element Analysis / 3.2.5 40 E l i m i n a t i n g (dv/d<P) f r o m Eqs. (3.12) a n d (3.22), y i e l d s the d i f f e r e n t i a l e q u a t i o n f o r the problem, ,d?w (dw d<f>2 ' \ \d<f> where V _ 2 — — ) d<f> d<j> 2 i 2 — r r + — - 2 — - — + w ,.2 + 2Rw = f(<j>) (3.23) /(*) = 2R> [ l + C{N+ - RPa} -T h e c y l i n d r i c a l membrane is fixed at the edges, r e s u l t i n g i n the f o l l o w i n g b o u n d a r y con-di t i o n s , n amely »(*!) = v{92) = «;(*!) = u»(*2) = 0 (3.24) T h e s t r u c t u r a l analysis of a pressure-loaded s t r i n g is hence defined by the n o n l i n e a r equa-t i o n d e p i c t e d i n Eq. (3.23), used i n c o n j u n c t i o n w i t h the b o u n d a r y c o n d i t i o n s given i n E q . (3.24). S o l u t i o n to the non l i n e a r e q u a t i o n w i l l be ob t a i n e d v i a a finite element scheme. There are however two unknown q u a n t i t i e s i n E q. (3.23) and these are the hoop stress , and the r a d i u s of cu r v a t u r e R((f>). T h e y w i l l be det e r m i n e d next. 3.2.5 H O O P S T R E S S , T h e hoop stress, is o b t a i n e d by i n t e g r a t i n g E q. (3.22) and i n v o k i n g t h e b o u n d a r y c o n d i t i o n s l i s t e d i n Eq. (3.24) to e l i m i n a t e the v. T h i s yi e l d s , J {27? [ ( $ ) + w + d<f> C*?R{4>)d<f> *1 (3.25) Observe f r o m Eq. (3.11) t h a t N$ is a constant a n d can thus be taken out of the i n t e g r a t i o n sign as shown i n Eq. (3.25). Two Dimensional Finite Element Analysis / 3.2.6 41 i i Deformed Configuration F I G U R E 3.3 Deformation of an Initially Circular String. 3.2.6 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 = (R0 + rv) cos <j> + v sin <j> (3.26) z — (RQ + to) sin <j> — v cos <f> (3.27) 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 and %=($)/ (10 - {%) (0) / (10 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>) = L 3/2 -(RQ + 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 strain-displacement relations, the v motions in the above quadratic terms are neglected in com-parison with the tw motions. Thus we have, 13/2 RW = , „ ,.,„. . ,.,a (3-29) {RQ + tw) tw" + (2tw'2) + (RQ + w) 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, i.e. w(<f>) = NjWj 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 F I G U R E 3.4 Typical Membrane Element. Minimizing the residue over the domain of a typical finite element results in the equa-tion 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> (3.33) K\NL) = Pi = <t>i j fNid* <t>i (3.34) (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 lw 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, dTi d awj uWj i.e. —± = - '* wk - 3.39 Defining the following quantities we have, Gii = ^ h J = 1, 2 (3.40) ^ d K ^ Hij = d * wk i, j, k = 1, 2 (3.41) 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 Ne finite elements, we have, Two Dimensional Finite Element Analysis / 3.4 45 E {it! ( 2 * [WW + NW)»i"j\ + + °R2Pa) #} Similarly, substituting Eq. (3.30) into Eq. (3.29) yields the radius of curvature, 3.4 EVALUATION OF FINITE ELEMENT MATRICES Invorder to be able to integrate some of the terms derived previously, it is necessary to 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, £{*[ f i ^W + ^ ) + % a « i ^ ' n=l ^ L + A^(tw1 + w2) + 2CR2paA<j>\ > n (3.46) 2C(*2 {R}n n=l Similarly, substituting Eq. (3.30) into Eq. (3.45) and computing the average radius of curvature we have, 3/2 (, x 2 I 3 / 2 " { + ( * o + " 1 ) 2 } { ( ^ ) 2 + ( * > + ^ 2 ) 2 } . — — - j - • — - = 2 ( ^ ) + (*o + "i) 2 2 ( = ^ ) + {Ro+ 102)1 (3.47) The various stiffness matrices, 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 iK(NL)] = T2 (3u,1 + u;2 + to^) («, 1 + u , 2 - 5 I ^ i ) Since we have, 5.-v = 1^2 .^ 21 ^22. where R + {3wl + tu2) ( u»1 + u>2) [5] = <S12 2^1 From Eq. (3.40), the stiffness matrix G,y is given by, -u>2) 2 -<o 2) 2 - W2) 2 (wl - « > * ) + 2R + 2R and this yields1 <7wf TJtcJJ L dwi dw-t where 7\ = R2 A# <C + • -R (P« ~ P J } (3.48) (3.49) (3.50) (3.51) 72 = RA<f> | 2 + C (2JV, - 3* P a) - 3 * f a ^ - P « ) j and (dN^/dwi), (dR/dwi) are defined in Appendix B. The stiffness matrix, is given 1 see Appendix B. Two Dimensional Finite Element Analysis / 3.4 4 7 from Eq. (3.41) as, Hi>' ~ -dw-(NL) WK But as shown in Eq. (3.49), K\^L^ is a linear function of wk. Thus we have, dK) (NL) ^ - W k = K^NL) Hence The external load vector P,, is obtained by from Eqs. (3.30) and (3.35) as, [H\ = [jr<™>] = 12 (3.52) Pl = P2 = R2&<f> (3.53) 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 displace-ments are then solved and added to the previous displacements to get the current displace-ments. 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.1 48 3.5 C O M P U T E R I M P L E M E N T A T I O N The governing finite element equations derived for the structural analysis were implemented in program N F E A M 2 3 . To capture the fluid-structure interaction effects, the finite element program is coupled to the boundary element program presented in Chapter 2. This coupled program is named I N F E A M 2 4 . Figure 3.5 depicts the structure of I N F E A M 2 . Starting with an assumed profile of the cylindrical membrane, the flow anal-ysis 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 C R I T E R I O N By specifying that the two profiles agree within a set tolerance, we are using a displacement-type 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 s Nonlinear finite Element Analysis of Membranes: 2-D * Interactive Nonlinear Finite Element Analysis of Membranes: 2-D Two Dimensional Finite Element Analysis / 3.5.1 49 F I G U R E 3.5 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 50 3.5.2 I T E R A T I O N 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 solu-tions. 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,6 fa is used in the solution algorithm. This idea of a line search 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).6 Basically, the concept of a line search algorithm involves seeking a new solution twj.+1 in a given descent direction dk as follows, at the fcth iteration: ™k+i = v>k + fad-k (3.56) where fa > 0 is the line search parameter that makes tuj.+ 1 an acceptable next iterate. The term line search refers to a technique for selecting the . 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 pk < 1, a retardation effect is produced in the iterative process. On the other hand, * 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. 5 Also known as an under-relaxation parameter when 0 < < 1, see Cook (1981) cit. loc, p. 359 6 cit. loc, p. 116 7 See Appendix D for more details 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 R E S U L T S This section is presented in two parts: non-interactive analysis from NFEAM2 and in-teractive 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 depicting fluid-structure interactions. The results are given in the part on interactive analysis section. 3.6.1 N O N - I N T E R A C T I V E A N A L Y S I S As the title implies, the fluid-structure interaction 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 wt is an unknown characteristic displacement introduced to non-dimensionalize the 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, Pw(4>)=P. cos(o^) (3.59) where p» is a characteristic wind pressure and a is an arbitrary parameter. This quantity pt has to be prescribed for a given wind. At low wind speeds, p. must be necessarily small while the internal pressure in the membrane pa, can be set arbitrarily large. This 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, wt=2R[^) (3.61) 8 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, 9" + e[9,2 + 92]+9 = Y£ 1 + C(JV^ - RPa) -Rpa(l — ecosoxp) 4> (3.62) where $2 N+ = c{*2- »l} / ^ ^ + ^ + 2eg}d* + Rp< (3.63) 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>) = gM) + egx{4>) + e2&{<!>) + 0{ez) N<j> = N0 + eN1 + e2N2 + O{e3) (3.64) (3.65) Substituting Eqs. (3.64) and (3.65) into Eqs. (3.62) and (3.63) result in the following sets of linear recussive relations, 9Z + g0=l2-[cN1 + $] + l1cos(a<fi) g'l + gx = I [CN2 + g - (&)2] - 1 ^ cos(a^) - g'02 - g02 g^ + g2 = i[cN3 + ^ + ( ^ y - 2 ^ + 2 [(ft) ~ ft cos(o^) - 2^1 - 2gQgx (3.66) and AT0 = RPa 1 C ( * s 2-^7)S%d4> $2 (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(e1) solutions, the number of terms requiring integration jumps to 20. For this reason, the solutions of the perturbation analysis will be obtained to at most 0(e2) 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 equa-tions': 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° * 2 = 150° C = | x IO"4 m/kg F I G U R E 3.6 Membrane Roof Subjected to Low Wind Pressure. 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 Finite Element Analysis o{*1) 0(S>) 300.00 307.74 303.62 304.34 Table 3.1 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 / 3.6.1 56 F I G U R E 3.7 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 — $i)- 9 This term vanishes at some combinations of $'s (eg. $ x = 0°, $ 2 = 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 mem-branes 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 cur-vature, 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 pa, in addition to the wind pressure pw. 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 U2 Two Dimensional Finite Element Analysis / 3.6.1 58 where pa is the density of air. Another quantity that will determine the deformation response of the mem-brane 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 defor-mation at a fixed f0. 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 FEM 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~ ooooo 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 A N A L Y S I S 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 mathemat-ical 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 mem-brane 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 61 dn = 0 i i dn y i \ 1 Rigid Cylinder / / / / 1 / 1 / \y ^-Flexible Cylinder >\ I*-\ \ \ \ \ ) If x / \ f — • L F I G U R E 3.10 Flow Domain and Boundary Conditions for Cylindrical membrane. 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 con-ditions 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, pa = 1.25 k g / m 3 , RQ = 10 m (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 No. of Boundary No. of No. of Finite Section Elements Singularities Elements Vertical 30 x 2 3 x 2 -Horizontal 48 4 -Membrane 50 6 49 Total 158 16 49 Table 3.3 Discretization Schemes The 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 £ a = 0.276 and £„ = 1.875 respectively. In contrast with a rigid 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 a , on membrane deflec-tions. The results are similar to Kunieda et al. (1981) who have adopted a real flow model 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 deflections10 in the present The inward deflection at the windward side and the inward deflection at the leeward side Two Dimensional Finite Element Analysis / 3.6.2 63 F I G U R E 3.11 Membrane Deformations at Each Global Iteration (£a = 0.276,C = | x 1CT4J. F I G U R E 3.12 Membrane Deformations at Each Global Iteration (£a = 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~4J. 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 a at different wind 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 / 3.6.2 65 F I G U R E 3.14 Influence of the Loading Ratio £a, on Membrane Deflections (C = between the total stress and the stress due only to the internal pressure, changes rapidly at small pa. It then decreases with increasing internal pressure. This is attributed to the fact that the membrane becomes more and more 'rigid' at higher pa. Also, recalling that 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 pa. A difficulty that quickly arises in the attempt 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 perturba-tion analysis, is obtained from the final pressure loads acting on the membrane. From Two Dimensional Finite Element Analysis / 3.6.2 66 to A ! < ft O o F I G U R E 3.15 Py, = 102 kfi/m! p w = 57 kg/m 2 P„ = 25kg/m8 p w = 6 kg/m* P w = 0 kg/m 200 250 Internal Pressure pa , kg/m 2 Variation of Membrane Stress Resultant , with Internal Pressure P a (c=\x i o ~ 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(e1) term for the hoop stress, Nl is derived from Eq. (3.67) to be, AT, = aZi {1 - cos (A$)} - Z2 sin (A$) + (a2 - l) aa,Z 3 1 (a 2 - I) a {(7A$sin(A$) + (C + ^)[2 - A$sin(A$) - 2cos(A$)]} (3.71) where Zx = cos(a$!) + cos(a$2) Two Dimensional Finite Element Analysis / 3.7 67 Z 2 = sin (a$2) — s m ( a $i) Z3 = A$sin(A$) - 2{1 -cos(A$)} The results of the perturbation analysis, accurate to 0(el), are shown by broken lines in 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 S U M M A R Y A two-dimensional interactive analysis of cylindrical pneumatic membranes subject to wind loads is presented. The flow, modeled with an inviscid fluid, is analysed using the boundary element technique. The structural analysis is based on a large deformation but small-strain 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 fa. Also the additional hoop stress due to wind changes rapidly at smaller pa, and then decreases with increasing pa. C H A P T E R 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 struc-ture 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 prin-ciple 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 el-ement procedure. The incremental continuum mechanics equations are used to develop the governing finite element equations. To evaluate the integral associated with the ex-ternal 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 dimensional finite element program has been developed for the struc-tural analysis. The program, called NFEAM3 1 uses the warped isoparametric quadrilat-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 au-tomatic 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 cou-pled program named INFEAM3, 2 is then used to analyse the air-supported roof of the BC Place Stadium under the action of wind loads. 4.2 M A T H E M A T I C A L P R E L I M I N A R I E S Only the relevant mathematical relationships are given here. For a more complete presen-tation, 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 quan-tities 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 a2 F I G U R E 4.1 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, 62) with respect to a fixed rectangular cartesian coordinates, x% 3 as follows, * < = x,'(0a) (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 = xigi (4.2) 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° ~ d9a 8 i = (4.3) where a*a = (dz'/<90o) yields the cartesian components of the base vectors. The contra-variant surface base vectors, a a are given by, aa.ap = 8$ (A A) where the Kronecker delta, Sjj is defined as, 5c = i 0 forct^P; ( 4 . P \ 1 for a = /?, P not summed. \ ' ) The covariant metric tensor on the surface is given by, °a/? = aa» a/? (4.6) This yields, fl«3 = a « a J 9ij = a/?a (4.7) Similarly the contravariant metric tensor on the surface can be computed from, a o A a A / ? = $| (4.8) The normal vector to the surface at P, a 3 is given by, a 3 = ai x a 2 (4.9) Substituting Eq. (4.3) into Eq. (4.9) yields, « 3 = A Si = <»i<*2& x 8j Three Dimensional Finite Element Analysis / 4.2 72 where a 3 denotes the cartesian components of the normal vector. 4® = <*1 «2 eijk8k = o'i «2 eijkV99ktgi 4 = « i °2 ei}k\/99H i.e. (4.10) where e.-yj. is the permutation tensor defined as, eijk = >/J7eiyjfc 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 0>\0>2 ~ a">a 2"1 a\a\ - a\a\ ,1„2 ,1„2 axa2 - a2a1 Dennoting the magnitude of the normal vector by a 3 we get from Eq. (4.9), a 3 = l«3l = >J4 4 9ij (4.12) (4.13) If e 3 represents the unit normal vector to the surface at P, we have, a 3 ax x a 2 a 3 |aj x a 2 | (4.14) where e'3 are the cartesian components of the unit normal vector and are given by, V. c 3 / > = — 4 t „2„3 ,2 -° 3 1 a 3 "3 <»3 ,2„3 -.1„3 a l ° 2 ~~ a l a 2 11 11 (4.15) The magnitude of the surface base vector for the covariant components is given by, = ( a a . a a ) 1 / 2 = v ^ (no sum) (4.16) Three Dimensional Finite Element Analysis / 4.2 and the contravariant components are given by, \aa\ = v/a™ (no sum) The direction cosine between the surface base vectors is given by, a 1 . a 2 = 182)003(8!, a 2) or a12 cos(au a 2) \ / 0 i i a 22 From trigonometry we have, sin(a l f a 2) = yjl - cos 2 ^, a 2) • / ^ / an°22" 8111(8!, a 2) = W — V <*n a12°12 La22 or sm(au a2) = J V an a 22 where a = |°a/3| = °lla22 - a12a12 The line element of the surface at P, ds is given by, (ds)2 = dr *dr (ds)2 = aad9a.aj3d9^ i.e. (is) 2 = aQp dOa dBp The area of a surface element is found from, dA = \drx x dr2\ dA = \d9l a, x dB2 a 2| i.e. <M = dfl1 <f02 \ax\ \a2\sin(a1, a 2) In view of Eqs. (4.16) and (4.19), dA is thus given by, dA = ^d9ld62 Three Dimensional Finite Element Analysis / 4.3 74 4.3 T H E O R E T I C A L 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, 1C and a neighboring configuration, 2C. This is shown in Figure 4.2. It is 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. Three Dimensional Finite Element Analysis / 4.3.1 75 4.3.1 C O N V E C T E D S T R A I N T E N S O R A symmetric tensor in a deformed configuration, VC can be defined as, "lap = i - V ) (4-23) where "a a /j, °a a ^ are the metric tensors in the deformed configuration, VC and the unde-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)2-(d0S)2 = 2''<1af}d9ad9P (4.24) A symmetric, incremental convected strain tensor, 7a/3 can then be defined from Eq. (4.23) as, la? = 2lafi ~ l1a(i (4.25) i.e. 7 o j 9 = ± ( 2 a a / ) - laafi) (4.26) Expanding the covariant metric tensor, 2aap in Taylor series about the deformed configura-tion, 1C we have, - v * + - '*•) + \VM (v - '*•') (V - V)+... Neglecting cubic and higher order terms we have, oQ / J - a a / l « - ^ p - x + 2 a l i , a l x i x x> (4.27) where x* = (2x* — 1xt) represents an incremental displacement term. In view of Eq. (4.27), Eq. (4.26) can be recast in the following form, 7o/3 = Cap + Vafi (4.28) where c a ( 3 is the linear incremental strain tensor given by, (<•»> and r)a/3 is the nonlinear incremental strain tensor given by, _ ! d2 (laafi) • • ^ ^ W ? 1 ^ ( 4 - 3 0 ) Three Dimensional Finite Element Analysis / 4.3.2 76 4.3.2 C O N V E C T E D S T R E S S 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 ten-sors, this convected stress tensor will be defined with respect to the undeformed configura-tion °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 configura-tion VC as shown in Figure 4.3. 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 dvSa along the curvilinear coordinates, 6a and simply duS along the third side of the triangle. The orientation of the sides are defined by unit outward normal vectors "v a and "v, as shown in Figure 4.3. When the triangle 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, duS = vvdvS or in terms of the contravariant base vectors, "a a we have, 2 vaa vvdvS = Y -7=dvSQ a = l (4.31) (4.32) where d vSa represents the length of the sides along the 6a coordinate lines given by, dvSa = y/^d9a (no sum) (4.33) If uva denotes the covariant components of "v with respect to the base vectors aa, we have, vvaV^dvS = dvSa For convenience, the following covariant quantity is defined, dvS„ Eq. (4.34) can thus be written as, Also from Eq. (4.8) we have, vvad"S = dvS*a v„a\ v„ ca which upon expansion and inverting yields, "a 1 1 "a 1 2' 1 "°22 "o21 "a 2 2 ~ "a -^21 (4.34) (4.35) (4.36) (4.37) (4.38) where a — an a22 — O12 a 2 i Three Dimensional Finite Element Analysis / 4.3.2 78 Let vt and "ta 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, vtdvS = %dvSa (4.39) 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 ; a v / V ' ' t a (4.40) Note in Eq. (4.40), vt is an invariant to coordinate transformation for any arbitrarily chosen "v and "va is a covariant vector. Thus it follows that \/vaaa vta is a vector with a contravariant character, a. We may therefore define the following contravariant quantity, "ta* = v^a 5""ta (no sum) (4.41) The above contravariant quantity, vta* can be decomposed with respect to either a contra-variant or covariant basis at P as follows, or "ta* = "E% V (4.43) vEa^ and VE^ are tensors by Quotient rule and are called convected-contravariant Euler 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, dvT = vtdvS (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, dvT = vtdvS = und°S (4.45) With the aid of Eqs. (4.36), (4.39) and (4.41), we can rewrite Eq. (4.45) as, d VT = vta* d "S* = vna* d °S* (4.46) 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, vta* = vEafi "ap = VE$ "a* (4.47) or "t a* = vFa? \ = VF$ V (4.48) and n = n K a^ = (4.49J or n = 9 a^ = a M (4.50) Of the many possible definitions of the stress tensors, two will be of special interest, namely the convected P-K stress tensor, una^ and the convected Euler stress tensor, vEa& . The 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 ba= n ^ apd ba from which "EaP d "S* = "n a / ? d°S* (4.52) But from Eq. (4.51) we have, 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, vt d VS = "ta* d "5* = vEa& "ap d "S* (4.54) Expressing the base vectors "a^ in terms of the unit vectors vep we have VH % = —J= (no sum) (4.55) v HP Manipulating with Eqs. (4.35), (4.38) and (4.55) we have, vapdvSl^^a%dvSa In view of above, Eq. (4.54) becomes, vtd"S = "Eal}\fia d "Sa ve0 (4.56) utdvS = vEapdl,Save0 (4.57) in which vEaf> = \fovEafi (4.58) Since VE°P is associated with the physical quantities dvSa and "e^ as shown in Eq. (4.57), vEAP 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 vtduS = vnd°S = YJ V*HP"EPd°Sa (4-59) P=i Substituting Eqs. (4.35) and (4.51) into Eq. (4.59) we have, "nd°S = vnafi^dvSave? "n d °5 = ^ d "Sa ve? (4.60) where vnafi = V ^ V * (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 R E L A T I O N S For any configuration UC, the convected P-K stress tensor and the convected strain tensor are related as follows, where vCa^pX is the elastic tensor at configuration "C. For isotropic material, the elastic tensor, as given in Green and Zerna (1968)4 is, 2(1 + /*) where E, \i are the elastic modulus and Poisson's ratio respectively, t is the membrane thickness and vaa$ is the metric tensor in configuration VC. The incremental convected (4.63) P-K stress tensor can be defined as, naP = 2Ga0pX 2^ _ l^pX 1 ^ i.e = 2C°^px l p X + ( 2 C a ? p x - l C a f i p x) l l p X (4.65) The elastic tensor 2Ca/3pX is an unknown apriori. However for small strain problems, we can assume negligible change in the elastic tensor in the increment. i.e. 2CapP\ w \Capp\ ^ 4 6 6 j Hence Eq. (4.65) becomes, Further simplification of Eq. (4.67) is necessary as ipX has a nonlinear component. This is done as follows: substituting from Eq. (4.28) we have, 4 cit. i oc , p. 389 Three Dimensional Finite Element Analysis / 4.3.4 82 The nonlinear incremental strain tensor, r\p\ is negligible within an increment and hence 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 E Q U A T I O N S O F E Q U I L I B R I U M A consistent and rigorous formulation of the incremental equilibrium equations based on virtual work approach is presented as follows. Let VA denote the area of the deformed surface "S loaded by a convected Euler stress field, vEaP. The principle of virtual dis-placement for a virtual field of convected strains, 6v"tap requires that, jj "E"? 6 (-7^) dA = »Wext (4.70) "A where vWexl is the external virtual work in configuration VC. Since the deformed surface 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 = ^Wext (4.71) °A The incremental virtual work equations are then derived by applying the above virtual work expression to two neighboring configurations, 1C and 2 C . Linearization of these equations results when the configurations are assumed to be infinitesimally close. In configuration 2 C , Eq. (4.71) becomes, / / V ^ ( 2 7 ^ ) d ° A = V e x t (4.72) °A Three Dimensional Finite Element Analysis / 4.3.4 83 Substituting Eqs. (4.25) and (4.64) into Eq. (4.72) yields, °A That is, Ij ( V * + 6(lafi) d°A = *Wexl (4.73) °A In view of Eq. (4.28), Eq. (4.73) becomes, II (ln^ + n^)6(eQ, + r}a,)dQA = 2Wext °A 11 n"? S^d°A + 11 V » 6Vap d°A = *Wext -Jj V 6eap (4.74) °A °A °A Using Eq. (4.67) to eliminate naP in Eq. (4.74) we have, /j lCa^x 7„A 6lafid°A + JJ ln* 6r,ap dQA = - JJ V 6ea/, dQA (4.75) °A °A °A which represent the nonlinear equations for the incremental analysis. The nonlinear term in the above equation is linearized as follows, 7PA hap = (epA + VP\) 6{eap + Vap) « epA Seafi (4.76) In view of the above, Eq. (4.75) becomes, JJlC^xepX6eQ0doA + JJln^6r,afld°A = *Wext - JJ^6eapd0A (4.77) °A °A °A elastic stiffness geometric stiffness unbalanced load From Eq. (4.77), the linearized incremental equations of equilibrium are obtained by ob-serving that the equation holds for arbitrary variations of compatible strains. These lin-earized 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 configura-tion. Solution to these equations is achieved via a finite element modelling. This will be presented next. Three Dimensional Finite Element Analysis / 4.4 84 4.4 I S O P A R A M E T R I C F I N I T E E L E M E N T S O L U T I O N A displacement-based finite element 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 VC is chosen for 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, 9a . Global coordinates, vxx and global displacements components, "J7* of P are given by, V = NM "fy (4.78) TP = NM Tty (4.79) where vx*N are the nodal coordinates of the ^ node in the I t h direction at configura-tion "Ify is nodal displacement defined in a similar manner and is the shape function of the A/ t h node. The shape functions for WIQS are, Nl = J(l + 0i)(l + 02) J V ^ J (1-0*) (1 + 02) (4.80) N* = \(i-el)(i-e 2) JV4= i ( l + 0i)(l-02) Substituting Eq. (4.78) into Eq. (4.2) yields the position vector of the point P, VT = N" vx\ g i (4.81) Differentiating Eq. (4.81) yields the covariant surface base vectors, a ~ d6a Substituting from Eq. (4.81) we have, dNM • "*a = -gga- v*\ gi = X g i (4.82) Three Dimensional Finite Element Analysis / 4.4 85 K2 «3 Si F I G U R E 4.4 Typical Finite Element. where dNN 1/ » "a* = (4.83) 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, °A V gafipX dlxi + i V « °A d ^Zd9ld92 (4.85) where from Eq. (4.22), d°A = V®ad9l d92. For external loads acting on a finite element surface, the external virtual work expression is given by, = JJ2V 8uiV2a~d91d92 (4.86) 2 A ext where V is the » t h component of the externally applied traction vector, 5u* is the * t h com-ponent of the virtual displacement vector and 2A is the deformed area in configuration 2C. The displacement components can be calculated as, u' = V - V = xi (4.87) Since the external load vector comprises only normal pressure, the traction vector can be written as, (4.88) 2*» _ 2 „ 2 .« t = p e3 where 2p is the magnitude of the normal pressure and 2e 3 is the i t n component of the unit normal vector, 2 e 3 in the deformed configuration. Substituting Eqs. (4.87) and (4.88) into Eq. (4.86) yields, 2Wext = JJ2P 24 6xis/*a'd9l d92 (4.89) Substituting Eq. (4.89) into Eq. (4.85) we have, JJ 14 [ dlx% dlxJ = JJ [«*•' 2 p 2 e 3 ] v ^ d * 1 d92-±JJ 2 [ dxx*dxxi J j °A ,i lap d(laap) 6x% lri dlx* 0ad91d92 (4.90) Three Dimensional Finite Element Analysis / 4.4 87 To simplify the incremental virtual work equations and thus facilitate computer program-ming, the following finite element matrices are introduced, A"1 0 0 N2 0 0 N3 •• [<?] = \Pa} = where NNa = dNM/d9Q 0 Nl 0 0 N2 0 0 0 0 0 0 JVl 0 0 N2 0 • • • N* n 0 0 N2, 0 0 N*a 0 N'a 0 0 N% 0 0 0 N}a 0 0 N2, 0 m — J I — / "V 1 "~ 2 i/_3 t / _ l i/_2 i/_3 i / 2 i/„ 3 \ r - \ . f - \ xl xl xl x2 x2 x2 • • • XA x\ x\ I uc\2 {%.} = { "a2 "a3 , In view of Eqs. (4.92), (4.93) and (4.94), Eq. (4.83) becomes, {%,} = [Pa] TO Also Eq. (4.84) becomes, "*afi = TO' [^] TO where [Pafi] = [P.] T fa] (4.91) (4.92) (4.93) (4.94) (4.95) (4.96) (4.97) Substituting Eqs. (4.91) and (4.93) into Eq. (4.78) and the result into Eq. (4.90) we have,5 V °A //(* 2 A "^a/?) lra^pA g(S*) ! pQ,/ 2 e 3 ] ^adOUQ2- II i <fiad9ld92\ (4.98) Three Dimensional Finite Element Analysis / 4.4 8 8 where 7 t h term otherwise 6t'=l l> a t 1 0, l Since the equation holds for arbitrary variations of the nodal coordinates 6 £ j , we have d2C d(l(1af)) lra0p\ di\\) a1? d^J //(' °A V de2} zJ °A de2 (4.99) ext Eq. (4.99) represents the incremental equations of equilibrium for an element. However, it is impossible to evaluate the first term on the l.h.s. of Eq. (4.99) because the quantities associated with the deformed configuration 2C are not known apriori. The situation is further complicated by the fact that the external loading is deformation dependent. The problem of converting the integrals over 2A to integrals over lA have been discussed by Oden and Key (1970), Oden (1972)6, Larsen and Popov (1974), and Bathe, Ramm and 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 lC we have, 2Wext — 1Wext + d(lWext),j + (higher order terms) (4.100) Examining the second term, we have i.e. d(lWext) e J _ \ [ t o K p ) ' ^ 4 . ^ l c ^ , , , £ f V 5 ) (4.101) 6 cit. he, p. 242 Three Dimensional Finite Element Analysis / 4.4 89 Observe that this second term represents a stiffness-like quantity, arising from the follower-load character of the pressure loading and hence contributes to the system tangent stiffness matrix. Unfortunately, this is an unsymmetric contribution. This unsymmetric load stiff-ness 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. How-ever, 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, 2Wext « / / [lp <?„ V 3 ] v W d92 + I jj 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, ira/3p\ d( l f lpA) + -JJ E*V A + P-Q^ry/la+ P e^-dT^r d6ld92 = ff[lpQ»1*] ^d^d92-Jf i °A d1? 0ad9ld92 (4.103) Eq. (4.103) depicts the linearized incremental equilibrium equations derived using a vir-tual 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.7 The unsymmetric stiffness terms were 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 - F I (4.104) where K.ffi, Kjp , Kjp are the elastic, the geometric and the unsymmetric stiffness matrices respectively and are defined as, K (E) IJ iff °A u d(laaf}) lrafiPX dC«p\) dH1 0ad9ld92 K (G) IJ °A 0ad91d92 (4.105) (4.106) K IJ -ffo \d ^ P ) V V ^ , t n ^ ' ^ ^ I ^ n V d ^ -JJ Q i I W e * V a + p~dqrVa+ P e^-d^J-d9ld92 (4.107) and Pj, Fj are the load vectors due to external and internal actions respectively. They are defined as, = ff[lpQul<&] V^Zd9ld9' lA Fi °A l„a0 d( l a «/?) \fia~d91d92 (4.108) (4.109) 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). Three Dimensional Finite Element Analysis / 4.4 91 and QI^T^J = (P°PKL + P<XPLK) 6I 6J where Papu is the (J, J ) t h term of [Pap] • Thus we have, (4.110) (4.111) * ffl^Kl^M 1Ca^(>X (PafJKL + Pa/UK) (PpXMN + PpXNM)S}S^V^Ld9ld92 Ku = \jj (P*f*KL + PaPLK) ti SfV^dd1 d92 °A Pi = i /j ^ (PaPKL + PafiLK) ti ^ (4.112) (4.113) (4.114) 0 ,1 in which lCALIPX is given by Eq. (4.63) as, IQOPPX _ Et 2(1 + n) I and lna/3 is given by Eqs. (4.23) and (4.62) as, \„ap 1 f)X , l„aX \„fip , (-5L) V> (4.115) i.e. (4.116) 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 - lA -"3 yfa which upon differentiation yields, 3 l*s9u dl£J l a 3 e5 3$ a 1 ^ e3 e3 9kl From Eqs. (4.10) and (4.95) we have, a 1 ^ (4.117) Three Dimensional Finite Element Analysis / 4.4.1 92 where % = ejkl Vggli PiJ + l4 P{j} (4.H8) The third term in Eq. (4.107) is evaluated from Eq. (4.20) as, djyfia') = 1 dla = , Qp d{laafi) 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, KIJ = JJ QiJ P l l ^ \ e3 e3 S j - e 3 e 3 SjJ + ^3 1 f l«/»l^ ( P a p K J + P a , J K ) V^d9ld92 (4.120) Having developed the governing finite element equations for the incremental analysis, we shall discuss next, the solution strategy. 4.4.1 E Q U I L I B R I U M I T E R A T I O N 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 fcth iteration is, Error*** = 2Wext - jj 2n°^ 6 (27$) 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 rep-resents 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 fcth iteration is, Jj i<7Q^ A A e J J 6eafi d°A + JJ lna? 6*$ d°A °A = %ezt - JJ V**"1) 6 d°A (4.122) Substituting the various finite element arrays and after some manipulations similar to that in deriving Eq. (4.103), we have, {// [i (^l<^^) + l f r $ $ ) ] J * " 1 * = JJ[lpQulei] Vhidelde* d9lde2 i ff Uap I d(2«ap) . X 2]] \ d*£' + d 2i Id^^ K I °A ' J 0ad9id92 (4.123) (4.124) where the incremental displacements are updated as follows, £j(fc) = £J{k-i) + A£J(k) 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, {Kff + K<ff - jrg>} A(W = P, _ rf-i) (4.125) where the arrays PI} K{jf , fljj* and A J ^ are given by Eqns. (4.108), (4.112), (4.113) and (4.120) respectively. The array Fjk~^ is obtained by substituting Eqs. (4.110) and (4.111) into Eq. (4.123) to yield, i f _ 1 ) = J / / ( V [ 2 ? ( P A F I K L + P A F I L K ) 6} °A + {PapKL + PapLK)ti^jiJ\fk ^ ^~ad9ld92 (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 un-balanced 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, in-voke 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 com-mand 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 fash-ion, 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. 8 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 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 / 4.5 96 START > begin clock • check RESX\RT code > set up all arrays w.r.t. a l-D master array I READ INPUT DATA • job title & control information • nodal coordinates • element connectivities • boundary conditions • material information • loading parameters I DATA ANALYSIS • check memory allocation for data input • diagnosis by DOCTOR • data organisation • printing Fatal Errors START LOAD INCREMENT ITERATION I CONSTRUCT GLOBAL QUANTITIES • compute element stiffness matrix • compute element unbalanced load vector • assemble elemental quantities into global I CONSTRUCTION TECHNIQUE > 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 I CORRECTIONS • improvement algorithm • line search parameter SOLUTION OF EQUATIONS No < NITER • choice of solution technique via ISOLN > use of a compacted unsymmetric solver, SOL No INVOKE > NITER RESTART FACILITY UPDATE LOAD INCREMENT F I G U R E 4.5 STRAINS & STRESSES ' compute and print strains i compute and print stresses Structure of NFEAM3. OUTPUT 1 print results & stop clock > print program statistics > print errors (if any) 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 SK in A is required. This is achieved by establishing the column height array L M and the connectivity array K L D . 9 Thus we see that the efficiency of this skyline 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, ISOLN .EQ. 1 - Triangularization of Stiffness Matrix .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 tech-nique. 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 cri-terion. 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 1 Q. J ^26 0 I K$2 K33 ^35 0 i 1 ^ 1^41 0 -^43 -^44 -^45 0 < j >. ^ 0 0 J -^53 -^54 -^55 -^56 (a) Actual Stiffness Matrix (Unsymmetric). SKD = (#11, #221 ^33> #44> ^ 5 5 » ^ 6 6 ) S K U = (#12, #23) ^14> 0> ^34 ) ^35> ^45> ^26> Of °i ^Se) SKL = (#21, #32» -Ktl 1 Of -^43J ^531 ^ 5 4 > ^ 621 0> 0> ^65) (b) Arrays for Storing Elements of S K . F I G U R E 4.6 Storage Scheme used for Global Stiffness Matrix. Three Dimensional Finite Element Analysis / 4.5 99 energy10 during each iteration is compared to the initial internal energy increment. Bathe 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 ZJW = £ J ( * _ 1 ) + BETA<*) A f ' W (4.128) Also, incorporated in NFEAM3 is a RESTART facility. If for some rea-sons, 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 execu-tion 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.1 100 4.6 N U M E R I C A L R E S U L T S As in the two-dimensional analysis, this section is presented in two parts: the non-interactive 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 N O N - I N T E R A C T I V E A N A L Y S I S 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 exam-ple 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 - 3 in E = 15 x 106 psi (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 / 4.6.1 101 6 C 4 5 B 1 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 ° a d increment for convergence is found to be 0.459 while the average number of iterations required to achieve convergence at each load step is 5.3. Three Dimensional Finite Element Analysis / 4.6.1 102 Deflection in inches Points on Membrane 1 2 3 4 5 6 Linear Theory 0.164 0.323 0.366 0.691 0.798 0.928 Irvine (1976) 0.135 0.264 0.299 0.565 0.653 0.759 3-D FEM 0.142 0.268 0.296 0.572 0.656 0.763 Table 4.2 Comparison of Deflection Response for the Square Membrane A comparison of the deflection response is presented in Table 4.2. For con-venience, 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 ana-lytical 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 incremen-tal 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 / 4.6.1 103 F I G U R E 4.8 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 Incre-ments: Square Membrane. Three Dimensional Finite Element Analysis / 4.6.1 104 Incremental Tension in lb/in Locations in Membrane A B C Linear Theory Irvine (1976) 3-D FEM 0.0 0.210 0.242 0.0 1.878 1.955 0.0 3.667 3.770 Table 4.3 Comparison of Tension Increments for the Square Membrane 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 conver-gence. 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 / 4.6.1 105 F I G U R E 4.11 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, pa = 1.25 kg/m 3, 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 F E M F I G U R E 4.13 Membrane Deformations for Example 1 (£, = 1.91, C = | x IO -4,). Hoop Stress Resultants, in Kg/m Uemura (1972) 2-D F E M 3-D F E M 818.7 864.1 886.3 Table 4.4 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 two finite element solutions. This suggests that the approximations introduced in the strain-displacement relations for the two-dimensional analysis are valid, as such approx-imations are not used in the three-dimensional analysis. Also, for reasons mentioned in Chapter 3, the agreement with Uemura's solutions is poor. U E M U R A (1972) ° ° 0 0 2 - D F E M FIGURE 4.14 Membrane Deformations for Example 2 (£a = 3.40, C = l x 10_4j. Hoop Stress Resultants, in Kg/m Uemura (1972) 2-D FEM 3-D FEM 1227.3 1271.9 1299.3 Table 4.5 Comparison of Hoop Stress Resultants for Example 2 Three Dimensional Finite Element Analysis / 4.6.2 109 4.6.2 I N T E R A C T I V E A N A L Y S I S Having studied the performance of the three-dimensional finite element program, we con-sider 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, pa = 1.20 kg/m 3 (4.132) and for the roof, Li = 182.9 m, Bx = 224.3 m, Hx = 26.8 m, H2 = 33.5 m (4.133) E = 0.690 GPa, fi = 0.30, t = 7.62 mm, pa = 345 Pa 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 B o u n d a r y No. of B o u n d a r y No. of No. of F i n i t e S e c t i o n Nodes S i n g u l a r i t i e s E l e m e n t s R o o f 50 6 38 W a l l 60 4 -F l o w D o m a i n 171 14 -T o t a l 281 24 38 Table 4.6 Discretization Schemes F i n a l l y , t h i s r e p o r t is d i v i d e d i n t o three stages: (a) I n f l a t i o n analysis (b) F l o w m o d e l i n v e s t i g a t i o n (c) I n t e r a c t i v e a n a l y s i s (a) Inflation Analysis To i n i t i a t e the s o l u t i o n process f o r membrane deformations due to w i n d , the i n f l a t i o n p r o b l e m is first solved. T h i s involves c o m p u t i n g the s t r u c t u r a l response of the roof due only t o i n t e r n a l pressure, pa. T h e undeformed roof geometry w h i c h is stress-free, is shown i n Figure 4.17.11 T h e d a t a set p e r t a i n i n g to t h i s undeformed roof is t a b u l a t e d i n Appendix E. It s h o u l d be emphasized here t h a t the f o l l o w i n g t e r m i n o l o g y is adopted: undeformed profile refers to the state p r i o r t o i n f l a t i n g the roof, w h i l e initial p r o f i l e refers to the state after the roof 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 he p r e s c r i b e d f u l l l o a d is reached. A l s o at each l o a d increment, a 'lo c a l ' i t e r a t i o n is set u p t o compute t h e corresponding converged pr o f i l e . T h i s converged p r o f i l e t h e n becomes the s t a r t i n g p r o f i l e for the next l o a d increment i t e r a t i o n . T h e average number of 'local' i t e r a t i o n s i s 3.1 w h i l e the average l i n e search p a r a m e t e r encountered f o r the 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 / 4.6.2 113 F I G U R E 4.17 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 roof12 is shown in Figure 4.18. As expected, Nx appeared to be fairly constant while Ny increases towards the apex of the roof. At the apex, the shear stress resultant, Nxy 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 Nx and Ny roles are interchanged. Hence, along this axis, Ny is fairly constant while Nx 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 Elliptic-Paraboloid Finite Element Resultants Shell Theory Solution kN/m kN/m Nx 44.10 47.01 Ny 33.22 36.60 0.0 0.0 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), flow separation over a large portion of the roof was encountered. Three Dimensional Finite Element Analysis / 4.6.2 116 F I G U R E 4.19 Model Test of BC Place Stadium (Parkinson (1980)). Three Dimensional Finite Element Analysis / 4.6.2 117 POTENTIAL FLOW REAL FLOW (PARKINSON (1980)) (A) S E H K o I—I w a o. o INITIAL PROFILE POTENTIAL FLOW REAL FLOW (PARKINSON (1980)) m — •« « * i ft n r; n r i Ti .t .-; »• - •« „ (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 mi-nor axis are shown in Figure 4.22. Again, we see for the case of the potential flow model, that Nx is fairly constant while Ny increases towards the apex. Similarly, Nxy 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, Nxy Three Dimensional Finite Element Analysis / 4.6.2 118 CM P O T E N T I A L F L O W R E A L F L O W ( P A R K I N S O N (1980)) CO I" « ID W Ui w Pi «r> cu i i (A) O . OO o CO o I N I T I A L P R O F I L E P O T E N T I A L F L O W R E A L F L O W ( P A R K I N S O N (1980)) • .*» ."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 Elliptic-Paraboloid Shell Theory kN/m Potential Flow Model kN/m Real Flow Model kN/m 74.17 91.53 78.91 *y 50.52 58.87 69.49 Nxy 0.0 0.0 4.37 Table 4.8 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 w OT RIGID MODEL FLEXIBLE MODEL OT o. o INITIAL PROFILE DEFORMED PROFILE (B) o o. OT OT W RIGID MODEL FLEXIBLE MODEL (C) S .S o. 00 to CM INITIAL PROFILE DEFORMED PROFILE F I G U R E 4.23 (D) 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) F I G U R E 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 / 4.6.2 124 F I G U R E 4.26 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 two-dimensional interactive analysis in which the flow past an infinitely long cylinder was studied. Recall in those results14 the symmetry in the deformed profile was destroyed. This 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 three-dimensional 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 coefficient15 for a circular cylinder is -3.0 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 Nx is generally constant while Ny increases towards the apex. At that point, Nxy 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 Wind Speed V = 15 m/s Wind Speed V = 25 m/s Wind Speed V = 35 m/s Nx 54.82 81.25 110.35 Ny 41.35 53.78 70.84 0.10 0.13 0.28 Table 4.10 Stress Resultants at the Apex of the Flexible Roof defined as the ratio of the maximum pressure to the dynamic pressure of the flow, jp a U2 Three Dimensional Finite Element Analysis / 4.6.2 127 F I G U R E 4.28 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 / 4.7 128 Wind Speed, m/s Maximum Nxy, kN/m Locations in Roof X Y Z 15 11.40 -41.0 71.1 17.0 25 20.48 -65.1 82.3 4.62 35 31.65 -65.1 82.3 4.66 Table 4.11 Maximum Shear Stress Resultants and their Locations 4.7 S U M M A R Y A three-dimensional interactive analysis of pneumatic membrane structures is presented. The flow, 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 mem-brane. 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 blow-ing at the three selected wind speeds are also presented. It is shown that if the effects of fluid-structure interactions are ignored, the wind loads computed can err on the non-conservative 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 account fluid-structure interaction. 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 anal-ysed 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 de-scription 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 struc-tures: 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 A N A L Y S I S 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 formu-lation 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 equa-tion in two and three-dimensions. They are, deflection of a stretched membrane, flow past a rigid circular cylinder, heat conduction in a rectangular prism and flow past a rigid sphere. To accurately gauge the performance of this boundary-based1 numerical method, as opposed to the finite element 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 A N A L Y S I S 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 T W O - D I M E N S I O N A L 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.3 133 5.2.2 T H R E E - D I M E N S I O N A L 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 governing finite element equations. These finite element equations are solved using Newton-Raphson 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. Com-parison with published results is made, and good agreement is attained in both the exam-ples. 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 A N A L Y S I S Having discussed the performance of the boundary element and the finite element pro-grams, interactive results for wind-loaded pneumatic membrane structures are then gener-ated 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, on membrane deformations and the variation of membrane stress resultants, with the Summary and Conclusions / 5.4 134 internal pressure, p0 are presented. It is seen that the additional stress due to wind blowing at a specific velocity, changes rapidly at smaller pa, and then decreases at increasing p0. This is because the membrane structure has becomes so 'rigid' that there is little 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 the fluid-structure interaction 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 for fluttering and 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 Cylin-drical Membranes J. Struct. Div., ASCE, Vol. 109, p. 297-313 AHMADI, G. and GLOCKNER, P.G. (1984) Effect of Imperfections on Ponding Instability J . Eng. Mech. Div., ASCE, Vol. 110, p. 1167-1173 ALLGOWER, E. and GEORG, K . (1980) Simplicial and Continuation Methods for Approxi-mating 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 Pacific J. Math. Vol. 16, p. 1-3 ASCE (1979) 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. BATHE, K . J . (1982) Finite Element Procedures in Engineering Analysis Prentice Hall, New Jersey, 735 p. BATHE, K.J . , and CIMENTO, A.P. (1980) Some Practical Procedures for the Solution of Non-linear 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 Defor-mation 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 Defor-mations 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 AIAA Journal, Vol. 10, No. 8, p. 1107-1108 BIRD, R.B., STEWART, W.E. and LIGHTFOOT, E.N. (1960) Transport Phenomena 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 Com-puting Science Technical Report #67 Bell Lab., Murray Hill, N.J., 10 p. BREBBIA, C A . and DOMINGUEZ, J . (1977) Boundary Element Methods for Potential Prob-lems Appl. Math. Modeling Vol. 1, p. 372-378 136 CHRISTIANSEN, S. (1967) On Kupradze's Functional Equations for Plane Harmonic Prob-lems Function Theoretic Methods in Differential Equations (Eds. R.P. Gilbert and R.J. Weinacht), Pitman, London, p. 205-243 COOK, R.D. (1981) Concepts and Applications of Finite Element Analysis 2nd Ed., John Wiley k Sons, N.Y., 537 p. CRISFIELD, M.A. (1979) A Faster Modified Newton-Raphson Iteration Computer Meth-ods in Applied Mechanics and Engineering, Vol. 20, No. 3, p. 267-278 CRISFIELD, M A . (1983) An Arc-Length Method including Line Searches and Accelerations Int. J. Num. Meth. Engrg., Vol. 19, p. 1269-1289 DENNIS, J.E. Jr. and SCHNABEL , R.B. (1983) .Numerical Methods for Unconstrained Opti-mization and Nonlinear Equations Prentice Hall, New Jersey, 378 p. DENT, R.N. (1971) Principles of Pneumatic Architecture The Architectural Press, Lon-don 236 p. DHATT, G. and TOUZOT, G. (1984) The Finite Element Method Displayed (transl. CANTIN, G.) John Wiley & Sons, 509 p. DONNELL , L.H. (1976) Beams, Plates and Shells McGraw Hill, New York 453 p. F L U G G E , W. (1973) Stresses in Shells 2nd Ed., Springer Verlag 525 p. GEIGER, D.H. (1970) U.S. Pavillion at Expo 70 features Air Supported Cable Roof Civil Eng., ASCE, p. 48-50 GOLDSTE IN, A.A. (1967) Constructive Real Analysis Harper ic Row, N.Y., 178 p. GREEN , A.E. and ADKINS, J.E. (1970) Large Elastic Deformations 2nd Ed., Oxford Univ. Press, London 324 p. G R E E N , A.E. and ZERNA, W. (1968) Theoretical Elasticity 2nd Ed., Oxford Univ. Press, London 457 p. H A N , P.S. and OLSON, M.D. (1986) Static Analysis of Pneumatic Structures Loaded by Wind submitted to IASS Int. Symp. Membrane Structures and Space Frames Osaka, Japan HAN, P.S., OLSON, M.D. and JOHNSTON, R.L. (1984) A Galerkin Boundary Element Formu-lation with Moving Singularities J. Eng. Comput. Vol. 1, p. 232-236 HART-SMITH, L.J. and CRISP, J.D.C. (1967) Large Elastic Deformations of Thin Rubber Membranes Int. J. Eng. Sc., Vol. 5, p. 1-27 H A U G , E. and POWELL , G.H. (1972a) Finite Element Analysis of Nonlinear Membrane Struc-tures Rep. No. UC-SESM 72-7, Univ. Calif., Berkeley H A U G , E. and P O W E L L , G.H. (1972b) Finite Element Analysis of Nonlinear Membrane Struc-tures Proc. IASS, Pacific Symp. Tension Struct. Space Frames, Tokyo and Kyoto, p. 165-175 HEISE, U. (1978) 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 137 HOPPER, M.J . (ed.) (1973) Harwell Subroutine Library Catalogue Theoretical Physics Division, A.E.R.E. Harwell, England 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 AHMAD, S. (1980) Techniques of Finite Elements Ellis Horwood Ltd., Chichester, 529 p. IRVINE, H.M. (1976) Analytical Solutions for Pretensioned Cable Nets J. Eng. Mech. Div., ASCE, Vol. 102, p. 43-57 IRVINE, H .M. (1981) Cable Structures 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 HAN, P.S. (1984) The Method of Fundamental So-lutions, an Adaptive Boundary Element Method, for Problems in Potential Flow and Solid Mechanics Proc. Fifth ASCE Specialty Conf., Eng. Mech. Div., 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 Mem-branes Z. angew. Math. Phys., Vol. 15, No. 608, p. 608-628 KUNIEDA, H., YOKOYAMA, Y . and ARAKAWA, M . (1981) Cylindrical Pneumatic Membrane Structures Subject to Wind J. Eng. Mech. Div., ASCE, Vol. 107, p. 851-867 KUPRADZE, 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. Engi-neers, 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 Spher-ical 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 Mem-branes 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, (eds. W. and B. Pilkey), The Shock and Vibra-tion Monograph Series MORRIS, N.F.(1981) Wind-Structure Interaction of Flexible Roofs Proc. Symp. Long Span Roof Structures, Struct. Div., ASCE, p. 297-312 NAGARAJAN, 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 J. Struct. Div., ASCE, Vol. 93, p. 235-255 ODEN, J.T. (1969) Discussion of "Finite Element Analysis of Non-Linear Structures" J. Struct. Div., ASCE, Vol. 95, p. 1379-1381 ODEN, J.T. (1972) Finite Elements of Nonlinear Continua 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 Struc-tures 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 J. Eng. Mech. Div., ASCE Vol. 94, p. 79-101 OLSON, M.D. (1982) CIVL 538: Advanced Topics in the Finite Element Method, Class Notes University of British Columbia OLSON, M.D. and HAN, P.S. (1984) Interactive Analysis of Wind-Loaded Membranes Proc. 2nd Int. Conf. Num. Methods Nonlinear Problems, Barcelona, p. 501-512 ORTEGA, J.M. and RHEINBOLT, W.C. (1970) Iterative Solution of Nonlinear Equations in Several Variables Academic Press, N.Y., 572 p. OTTO, F. (1962) Zugbeanspruchte Konstruktionen Ullstein Fachverlag, Berlin. Translated and Published as Tensile Structures MIT Press, Cambridge, Mass. 320 p. 139 PARK INSON, G.V. (1980) Wind Pressure Distributions on a Model ofB. C. Place Amphithe-atre for Phillips Barratt Engineers <fe Architects, 23 p. PATTERSON, C. and SHEIK, M.A. (1981) Regular Boundary Integral Equations for Stress Analysis Proc. Third Int. Sem. Boundary Element Methods Irvine, Ca. p. 85-104 PATTERSON, C. and SHEIK, M.A. (1982) 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 O W E L L , M.J.D. (1970) A Hybrid Method for Nonlinear Equations Numerical Methods for Nonlinear Equations, (eds. Rabinowitz, Gordon and Breach), p. 87-114 ROSS, E.W. Jr. (1970) Large Deflection of an Inflated Cylindrical Tent J. Applied Mech., ASME, No. 69-WA/APM-9, p. 845-851 SIMONS, J.W. and P O W E L L , G.H. (1982) Solution Strategies for Statically Loaded Nonlinear Structures UCB/EERC-82/22, Univ. Calif. Berkeley STEVENS, H.H. (1942) Air Supported Roofs for Factories Architectural Record, U E M U R A , M . (1972) Membrane Tension and Deformation in Cylindrical Pneumatic Struc-tures Subject to Wind Loading Proc. IASS, Pacific Symp. Tension Struct. Space Frames, Tokyo and Kyoto, p. 199-210 ZIENKIEWICZ.O.C. (1977) The Finite Element Method 3rd Edit., McGraw Hill, London, 787 p. 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 71 in these results, the equations for Dirichlet and Neumann boundary conditions are realized. The Laplace's equation in /-dimensional space is given by, d2u d2u d2u d2u 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 or + i w £ { c ^ ( t " x ) } ) ] , . < ' r = o S/L r ifil('w +*ww r, 1 See Eq. (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 r . Ckd(/dtlk (a(x) + tf(ty,x)) d r = 0 \<k<N (A.Z) where x — (xi, x2) are the coordinates of the boundary node tk = (tlk, t2k) are the coordinates of the singularity ^ , x ) = - log V ^ - ^ + ^ - x ^ 2 R + >/'?* + <2* 2(<ijfc - xx) + 2* life d$ 2(r2Jt - x2) 2k + 2t 2k dt2k [tlk - xxf + (t2k - x2)2 R^Tqk+t\k+qk 21(x)Ck H(tk,x) = /3(x)Ck$(tk,x) + (AA) . . 3 x i , . (hk-Xl)-g^ + (t2k-X2) ^ dx2 (<lJk-Xl) 2 + (*2Jfc-^2)2 The expressions (dxl/dn) and (dx2/dn) are the direction cosines with the xx and x 2 axes 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, M N EE/ f f(t*,x) ' Ckds/dtlk Ckd$/dt2k \{a(x) + n{tj)x)) KCkd$/dtlk, dT = 0 l<k<N, l>3 (A.5) Boundary Element Equations / A 142 where x = ( s l f x 2 , 1 3 , . . . , X/) denotes the b o u n d a r y nodes = (*iJti hk> hki • • • > '/it) denotes the s ingular i t ies ?(t f c,x) = [(t l f c - X t ) 2 + (<2* " * 2 ) 2 + • • • + ft* " I/)2] 3/2Jk (2-0 2 2-J - ( / - 2 ) ( t 1 J t - z 1 ) [(<lJfc - l l ) 2 + (*2Jb " *2? + ' •' + ftfc - I/)2] //2 -(* - 2)(*2jk - x 2 ) [ ( * i * - x , ) 2 + («2ib - *2) 2 + • • • + ft* ~ x / ) 2 ] ^l-l 1/2 + (' ~ 2)*2Jfc +4 + -- + 'f*[* + >/'?*+4 + - + 4] /-1 (A.6) - ( / - 2 ) ( t / J t - x , ) 1/2 + (hk ~ * i ) 2 + ihk ~ * 2 ) 2 + • • • + ft* - * / ) 2 ] l-l (/ - 2)7(x)CJfc ^(ftAlx) = /?(x)C i b f(t j f c >x) + ft* " * i ) 2 + (hk ~ *2)2 + ••• + ft* - x , ) 2 -.f/2 X ihk - x . ) ^ - + ft* - x2)-^ + • • • + ft* - xt) Q n dxi As before, the expressions (dxt/dn), (3x 2/9n),..., (dxi/dri) are the direction cosines with the xl, i 2 , . . i i 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 ? ( t y , x ) 3=1 (A.7) 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 l f x2, xz, . . . directions, namely, vXl, vx%, vX3, . . . are given by, du dXi That is, 3 = 1 2f ty - x . ) (l-2)(ti]-xi) The resultant velocity, V(x) is given by, • = 1,2 772 * = 1> 2,3, / > 3 dxi dx2 dxi V(x) = The Pressure, P(x) is given by Bernoulli's equation which is, P(x) = PQ(x) + ip{(Vb(x))2 - (K(x))2} (A.8) (A.9) where PQ(X) and K0(x) are the datum pressure and free stream velocity respectively. APPENDIX B Stiffness Matrix G t J r dP> The stiffness matrix G,y is defined in Eq. (3.40) as, OWj where P, is given by Eq. (3.53) as, P l = P 2 = R2A<p j l + C{N+ - RPa) - R[Va~^ P w ) | Since Pj = P2 , the subscript can be dropped. That is, From Eq. (-B.2) we have, dPi dP dwj du)j + = RA<j>^R^C + 2 + C{2N<j>-3RPa)-R(Pa-PwY N} J dtVj 3R(Pa " >g-p„)1 dR\ N+ J a»j J where from Eq. (3.46) we have, dw1 2C($ 2 -dR L, \KJn n=l dwi [9x(w\ + tof) + &«>i«>2 " 4CRzpa&<}>] + R2A<t>} (B.l) (B.2) (5.3) dR N t { 1 \ - E \ i [^i^i + ^ + &«>i«>2] + Atfun + u;2 + 2CJ22pa) J J (B.4) 144 Stiffness Matrix / B 145 dN. dw2 2C7($2 dR dw2 _ dR d w * t in which Q± and Q2 are defined as, 3 + A<f>2 f» = l 3A<(> 92 = A<f>2-6 3&<p From Eq. (3.47) we have, dR dwi (2fl2 + je2)2 dR dw2 I \(2Q2 + y ^ 2 + k\ ' (2& + z2y S o c 2 ) 2 (B.6) (2fi2 + £2)2 ) +2QZ + 4Q2R2- QZJ + Zl< in which R\ and £ 2 are defined as, Q = * l = (*b + «l) *2 = (*0 + «>2) (B.7) The stiffness matrix G,-y is obtained by substituting Eq. (B.Z) into Eq. (5.1) which yields, (B.8) ON* dN* "I 3ici 5 tea J a « dR n ai? dR where ?t = R2 &<p{C + } 72 = RA+ |2. + C - 3RPa) - 3RiPaN~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, R(Pa~Pu>Y 2RwH> + (wl + w2- Iw+Rt) + 2Rw = 2R2 where 1 + C{N+-Rpa)- (C.l) f{^(rvl + w2)+w + CR2Pa}d<p = ~ 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, + + ™2) + » = R [ l + C[N+ - RPa) - 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(il) = w($2) = 0 (C.6) 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 wt, given by, w(<f>) = tw, g(<*>) (C.7) Also, the wind pressure pw(<f>) can be non-dimensionalized as $(<f>), by defining a charac-teristic wind load p, as follows, PM) = P.$(<f>) (C.8) When a wind is prescribed, it fixes the value of p, and is thus a known quantity. Substi-tuting Eqs. (C.7) and (C.8) into Eqs. (C.4) and (C.5) we get, 1 + C(N<j> - RPa) R(pa-P.t(<l>)) N. (C.9) and *1 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, 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, s = Pa (C.12) (C.13) Perturbation Solutions In view of Eqs. (C.ll) and (C.12), the integrodifferential equation becomes, 9" + e[9,2 + 92]+9=lr 1 + C{N+ - RPa) -Rpa(l-V$(<t>)) (C.H) where C($ 2 - $o 1 (C.15) Perturbing g(<f>) and as follows, 9{<t>) = 9o(<f>) + Witt) + e292{4>) + 0{e3) N<j> = NQ + e N l + e2N2 + 0(e3) (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? + <7o2) + e2{92 + 92 + HA + 2<7O0i) + • • • = ^ [ { i + W - J R p J - * } + e{cjv1 + ^ ( ^ + ^ ) ) } +*> {CN2 - ((ft) 2 " ft - fe^w) } + *>{CN3 [ 2 f t ft - ( f t ) 3 - ft + ( ( f t ) 2 - ft) } + •••] and N0 + + e2N2 + • • • = Rpa + e < 9o2 + 92)d<f>}+--- (C19) We can set the first term on the r.h.s. of Eq. (C.18) to zero and this yields,1 {l + C(N0-RPa)-!fc}=0 (C + ^)(N0 - RPa) = 0 Since (C + ^ ) ^ 0 then (NQ - Rpa) = 0 N0 = Rpa i.e. (C.18) (C.20) Observe that Eq. (C.20) corresponds to the linear solution for the hoop stress of a pressure-loaded 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? + go2) + e2{g$ + g2 + + 2<70ffi) + • • • = 3 { C » i + ft + J f ( « } + | [CN9 + ft - (ft)' - ft^(*)} J { ^ 3 + ft + (ft)3 - 2ft ft + [(ft)2 - ft] J jW) } + • • • (C.21) + 1 Similar conclusion can also be arrived from Eq. ((7.19) 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? + 9o2) +s2{92 + 92 + ^ gWi + 2<7o0i) + • • • = h {CN, + ft + + l{cN2 + %- (ft)' - ft^)} + T H + ft + ( * ) 3 " 2 f t ft+ [ ( f t ) ' " ft] rf'>} + • • • 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) + e2?2($i) + ° ( £ 3 ) = 0 (CM) g0($2) + egi(92) + e2g2($2) + 0(e3) = 0 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='2 (cNi + ft) + | cos(c^) (C.26) where N0 = Rpa and Sb(*i) = 0o(*2) = 0 (C.27) Perturbation Analysis / C 151 Solution to Eq. (C.26) is given by, gQ{<p) = A0 cos <f> + BQ sin <f> + ^ {C + ^ ) - J_ cos(c^) (C.28) where \BQS 2s in($ 2 -* 1 ) l / l ( S 1 ) c o s $ 2 - / l ( $ 2 ) c o s $ l J l " J with Observe that Ni is unknown at this stage and will be determined in the 0[el) solution. 0(e1) Solution The governing equations from Eqs. (C.19), (C.23) and (C.24) are as follows, + 9i = J (CW 2 + ft " (ft ) 2) ~ ^ft C ° S ( ^ } " ^ " ^° 2 ( C - 3 ° ) where and <7i(*i) = <7i(*2)=0 (C.32) Substituting Eq. (C.28) into Eq. (C.31) yields the first order approximation to the hoop stress, namely, _ 2 Integrating yields afcosa^ + cosa$2}(l — cos(A$)) — {sina$2 — sina$x}sin(A$) "i — p v- (C33) (a 2 - l)a{CA$sin(A$) + (C+ ^)[2 - A$sin(A$) - 2cos(A$)]J where Perturbation Analysis / C 152 Solution to Eq. (C.30) is given by, gi(<t>) = Ai cos<f> + Bi sin <f> + a x cos(2a<^ ) + a 2 cos(a + 1)^ + a 3 cos (a — 1)<£ + a4sin(a + l)<f> + a$ sin(a — 1)^ + a^sin^ + ai<f>cos<f> + a% cosa<£ + a 9 (C.34) where - l °1 _ 8 ( a ' ^ - l ) ( 4 o 2 - l ) a3 = 2(a-2)to-l)a B o °2 - 2a ( a+ l ) ( a+2) °4 - 2a ( o+ l ) ( a+2) a 6 = - i ^ 0 ( C + 3 r } - ) (C.35) °8 - 2A r 0(a 2 -l ) 2(a2-l)2 ° S _ 2 ( a - 2 ) [ a - l ) a a 7 = i J V 1 B 0 ( C + ^ ) The integration constants AuBi in Eq. (C.34) are determined from the boundary condi-tions shown in Eq. (C.32) and this gives, {Bi J s in (* 2 -* i ) t / ^ ^ / 2 ( * 2 ) a n * i ~/2(^i)sin$ 2 (C.36) where L ) C O S $ 2 - / 2 ( $ 2 ) c o s $ l /2(^) = ox cos(2or0) + a 2 cos(a + l)<f> + o 3 cos(a — l)<f> + o 4 sin(a + l)<f> + a 5 sin(a — 1)^ + a 6 4> sin ^ + 07 <f> cos ^ + ag cos(ar<£) + a 9 (C37) Again N2 is unknown and will be determined in the 0(e2) solution. 0(s2) Solution Since it is intended that the results of the perturbation analysis be evaluated up to 0(el) 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 N2. This information is necessary for the complete solution of the O^e1) problem. Hence, we have from Eq. (C.19), Perturbation Analysis / C 153 N* = C{*22-*l)Ii9l+9'o2 + 92°)d4, ( a 3 8 ) Substituting Eqs. (C.28) and (C.34) into Eq. (C.38) we have, 2 V f •^ 2 = "FHTk *~T / S - ^ l c o s ^ + B\sm $ + a i cos(2a^) + a 2 cos(a + 1)# + a 3 cos(a - 1)<£ + a 4 sin(a + l)<f> + o5 sin(a - 1)<^ + a 6 ^ sin <f> + a 7 #cos<f> + a 8 cos(a^) + a 9 + 42 , B 2 , 0L(r. i\2 , <*2 + 1 cog(2afl + 2 ( ^ 1 ) Ka - 0 c o s ( a + W ~(a+l) C 0 8 ( a - *M + 2 ( ^ 1 ) Ka " J) s i n ( a + l)+ + (a + !) sin(« " + iV^o (C + jfe) cos^ + iV^o ( c + jfe) sin* - ^ t ^ f cos(a^) J <ty Integrating yields, 2/tt(sin $ 2 - sin $x) + 2fe2(cos $ 2 ~ c o s $ 1 ) + 2fe3 sin(A$) (C39) 2 ~ C A $ sin( A$) + (C + ^){2 - A $ sin( A$) - 2 sin $ x sin $ 2 - 2 cos $j cos $2} ' in which A $ = $ 2 - $ 1 ^1 = /3(^2) sin - / 3($i) sin $ 2 2^ = /3(^2) cos $! - / 3($i) cos $ 2 / i 3 = 61(sin 2a$ 2 — sin 2a$x) + 6 2 ( s m a $ 2 ~ s m + 63{sin(a + 1)$2 _ sin(a + l)$i} + fc4{sin(a — 1)$2 - sin(a - l)$i} + 65{cos(a + 1)$2 - cos(a + + 6e{cos(a - 1)$2 - cos(a - l ) $ i } + b7 (C.40) where /3(</>) and the constants bi... bj are defined as follows, h{4>) = ai cos(2a^) + a 2 cos(a + 1)* + a 3 cos(a — 1)# + 04 sin(ar + 1)* + 05 sin(a — l)<f> + a6<f> sin <f> + a7<f> cos<j> + a 8 cos(or^ ) - | (ft) ^ + B ? + £ ( c + A ) J + i ^ ± L ] Perturbation Analysis / C 154 and *1 = —a 4(a a-l)(4o; 63 = Ao [, 2(a+l)» V h = Bo f. 2(a+l)» V b7 = ' . = ^ { ' + 1 ^ } <C-42> - 5 0(COS$ 2 - C03$! - $ 2 8 m $ 2 + ^ 1 Sin - g ( ^ ) 2 A $ 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 algorithm1 incorporated 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, dk 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 Mhnost any starting point 155 A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 156 The line search parameter, fik in Eq. (D.l) exhibits these features in the solution iteration: Pk ' < 1 decelerating N-R iteration = 1 standard N-R iteration (£-2) k > 1 accelerating N-R iteration In well-behaved problems, (3k varies with the iteration cycle in this manner: usually, at the starting guess, {3k < 1. As the solution approaches convergence, f3k —• 1 and when the solution becomes very close to convergence, fik > 1. The implementation of the acceleration 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 prob-lems. 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 for-mulated, we shall convert the problem of solving a nonlinear equation to one of minimizing an associated function. This is achieved by forming the l2 norm of Eq. (D.3) which is, /(x) = i | F | 2 {DA) A new solution xk+it is obtained from the current solution xkt by requiring that /(xjt+i) < f(xk). This is achieved by computing xk+l in a descent direction dki from xk. Mathe-matically, dk is a descent direction from xk if the directional derivative of / at xk in the direction dk is negative, namely, V/(x f c )4 < 0 (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 fik > 0: f{xk+/3kdk)<f{xk) {D.6) 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 /3k > 0 namely, (3k should satisfy these two relationships: f{xk + (3k dk) < f{xk) + a(3k Vf{xk)dk {D.l) and V/ (x ) t + i K « V/Cxjfe + Pk dk)dk > XVf{xk)dk {D.S) where a G (0,1) and A G (a, 1) are constants. Eqs. {D.l) and {D.S) are satisfied simul-taneously 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 fol-lowing powerful result: that there exist j3k > 0 satisfying Eqs. {D.l) and {D.S) for any given descent direction dk and that any algorithm that generates solution xk obeying Vf(xk)(xk - xk+i) < 0. a Q d Eqs. {D.l) and {D.S) at each iteration is essentially globally convergent.3 A systematic procedure for computing the line search parameter fik is given next. Computation of the Line Search Parameter, (3k 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, f3k = 1 is first invoked. This is done to ensure 2 See Dennis and Schnabel (1983), p. 117 8 See Dennis and Schnabel (1983), p. 120-125 A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 158 y Permissible Values of {5k F I G U R E 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 (Jk = 1 produces an unacceptable next iterate, 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~4 in Eq. (D.l). This leaves only the f3k to be determined. Defining, f(P) = f(xk + {3dk) A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 159 we thus have the following known information: /(0) = /(**) and /'(0) = Vf(xk)dk (D.10) Since pk = 1 is used at the start of an iteration, the following information is also known: f(l) = f(xk + dk) (0.11) If Pk = 1 is not acceptable, namely f(xk+dk) 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 Pk. This is displayed in Figure 2. Thus we have for a general quadratic curve, y{P) = ap2 + bp + c where constants a, 6, c are evaluated by fitting the curve through the three points resulting in, y(P) = [/(l) - /(0) - f'(0)]P2 + />)/? + /(0) ( £ U 2 ) Minimizing Eq. (D.12) yields Pmin = — ~> ^ (0.13) 2[/(l)-/(0)-/'(0)] Note that in Eq. (D.13), since /'(0) < 0 and that /(l) > /(0) + a/'(0) > /(0) + /'(0) we A A A thus have Pmin > 0. Also observe that since /(l) > /(0) + af (0) we get, 1 If /(l) 3> /(0), then pmin < ^. This results in unnecessarily small Pmin 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 /?m,„ > pj at the first backtrack during an iteration. Thus the line search parameter at the first backtrack is given by, 1 > Pk = Pmin > TS ( D U ) * Le. /(I) > /(0) + a/'(0) A Globally Convergent Solution Algorithm for Newton-Raphson Iterations / D 160 A (0,/(0)) slope = /(0) y = [/(l) - /(0) - /'(0)]/32 + / W + /(0) 'min P = l P F I G U R E D.2 Quadratic Backtrack. 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{Pk) = f{xk + Pkdk) {D.15) The more accurate cubic curve will now be drawn through these four points and is given by, y{p) = apz + bp2 + f'(0)(3 + /(0) {D.16) where fa\_ 1 \l! -1 M) - /(0) - firPW 1 in which /?r, fit are two previous values of Pk. They need not necessarily be 1 and /3m,„ 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 2-3a/'(0) 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 t h iteration: xk+l =xk- J(xk)~l F(xk) (D.20) where J(xk) is the Jacobian of F at x^. Again using the l2 norm of F(x) as a measure of convergence we have, /(x) = i F ( x ) r F ( x ) (D.21) 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 * ) r dk < 0 We shall now show that the Newton direction at xk, —J(xk)~1 F(xk) is a descent direction for the minimization problem: Vf(xk)T dk = -FT(xk)J(xk)J(xk)-1F(xk) = -FT(xk)F(xk) < 0 (0.22) provided that F(xk) ^ 0. An alternative positive definite quadratic model can also be developed for the system of nonlinear equations but this will not be pursued here.5 Hence we see that the global technique for the one-dimensional problem can be applied to multi-dimensional problem if we define the objective function using the l2 norm and the descent 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(xk) is singular. We have to ensure that J(xk) is sufficiently well-conditioned resulting in J{xk)T J(xk) 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 BC Place Stadium is: Node Coord ina tes N X Y 2 1 91 . 440 0 . 0 0 . 0 2 82. 110 67.300 0 . 0 3 87. 780 44.870 0 . 0 4 90 . 530 22.430 0 . 0 5 74. 520 82.860 0 . 0 6 68 . S80 80.760 4 .410 7 68. 580 67.300 7.010 8 68 . 580 44 .870 12.250 9 68 . 580 22.430 13.430 10 68 . 580 0 . 0 13.850 11 54. 680 100.580 0 . 0 12 45. 720 82.260 10.550 13 45. 720 67.300 16.810 14 45. 720 44.870 19.070 15 45. 720 22.430 20.900 16 45. ,720 0 . 0 21 .550 17 28. ,850 110.450 0 . 0 18 22. .860 89.740 11.980 19 22. .860 67.300 18.720 20 22. .860 44.870 22.560 21 22. .860 22.430 24.730 22 22. .860 O.O 25.490 23 0. .0 112.170 0 . 0 24 0, .0 89.740 11.140 25 0 .0 67 .300 18.540 26 0 .0 44.870 23.750 27 0 .0 22.430 26.030 28 0 .0 0 . 0 26.830 29 -28 .850 110.450 0 . 0 30 -22 .860 89.740 11.980 31 -22 .860 67.300 18.720 32 -22 .860 44.870 22.560 33 -22 .860 22.430 24.730 34 -22 .860 0 .0 25.490 35 -54 .680 1OO.580 0 . 0 36 -45 .720 82.260 10.550 37 -45 .720 67.300 16.810 38 -45 .720 44.870 19.070 39 -45 .720 22.430 20.900 40 -45 .720 0 .0 21 .550 41 -74 .520 82.860 0 . 0 42 -68 .580 80.760 4 .410 43 -68 .580 67.300 7.010 44 -68 .580 44.870 12.250 45 -68 .580 22.430 13.430 46 -68 .580 0 . 0 13.850 47 -82 . 110 67.300 0 . 0 48 -87 .780 44.870 0 . 0 49 -90 .530 22.430 0 . 0 50 -91 .440 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, "Finite 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 Singularities," 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062907/manifest