osmvmou OF AN ARBITRARY TRIANGULAR ' pun: BENDING STIFFNESS MATRIX m ns APPLICATION To LARGE omncnou snail. PROBLEMS Thesis for the Degraefic! Ph D ‘ Ml.CHiG..£s.N STATE UNIVERSITY _ George Lasker 19.66 mass W "L LIBR/f?" I. Michigan Sta-“.6 University ’U E “I". i \xw: ‘. V: RL’I’A g‘éi-Y ABSTRACT DERIVATION OF AN ARBITRARY TRIANGULAR PLATE BENDING STIFFNESS MATRIX AND ITS APPLICATION TO LARGE DEFLEC TION SHELL PROBLEMS by George Lasker This investigation is concerned with a discrete model formulation and solution of certain shell structure problems. The discretization requires a modeling of the shell by flat triangular plate elements which can have any side lengths and thickness. Each element is assigned an independent set of deformed configurations and its elastic properties are described by two stiffness matrices associate; with membrane and bending stresses. Displacement and slape compatibility conditions are satisfied at the corner points of elements and in a limited sense these conditions are then satisfied along the common edge of adjacent members. The main objective of the investigation is to obtain an explicit representation of a bending stiffness matrix for an arbitrary triangular plate element and to examine its applicability to small and large deflection plate and shell problems. Two bending stiffness matrices are obtained. One is associated with a set of deformed configurations that permit compatibility between elements to be completely satisfied for some problems and to be satisfied to a high degree for problems in general. The other is associated with a set of deformed configurations that relaxes 810pe compatibility but appears to give better numerical results for some problems. GE OR GE LASKER The solution to the geometrically non-linear problem is obtained by a formulation consisting of a sequence of linear solutions which enable equilibrium conditions to be approximately satisfied with respect to the deformed configuration. A computer program is given for a class of axially loaded shell of revolution problems having a symmetrical or an asymmetrical deformed configuration describable by a half period strip. Most interpretive and computational operations are performed internally from a small amount of input .data describing the undeformed geometry, material properties, and boundary conditions. Numerical results are obtained for several shallow conical shells exhibiting a snap-through type of instability. These results compare favorably with both experimental and numerical results given by several other investigators. DERIVATION OF AN AR BITRARY TRIANGULAR PLATE BENDING STIFFNESS MATRIX AND ITS APPLICATION TO LARGE DEFLECTION SHELL PROBLEMS BY George Lasker A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOC TOR OF PHILOSOPHY Department of Metallurgy, Mechanics and Materials Science 1966 ACKNOWLEDGEMENTS I wish to express my gratitude to Professor William A. Bradley for the guidance and encouragement that he extended to me throughout this investigation. ii TABLE OF C ONTENTS ACKNOWLEDGEMENT ................. LIST OF TABLES ..................... LIST OF FIGURES ..................... NO TA TIONS ........................ I. INTR ODUC TION .................... l. ‘1 1 II. GENERAL FORMULATION ............. 2.1. Some general properties ........... 2. 2. Discrete method ............... 2. 3. Simplified model ............... 2. 4. Matrix relating nodes and members ..... 2. 5. Interpretation of non-linear problems . . . . 2. 6. Shell stiffness matrix . . . . . ....... 2. 7. Member-shell parameter transformation matrix... ........... III. TRIANGULAR MEMBER MATRICES ........ WWW WNH o. A symmetric form of the strain energy expression . . . .............. The member displacement configurations . . . Member-shell transformation matrix . . . . Member stiffness matrix . . . . . . . . . . wwww «1001A Continuity conditions between adjacent members . . ................ 1. Preliminary remarks ............ 2. Previous developments ........... 3. Present investigation ............ Member coordinate systems ......... Geometric parameters. . ......... . . Pa ge ii vi vii ix D—fi \Oyh" ll 11 16 24 25 26 29 31 38 3 8 4o 43 48 58 '72 76 Page IV. FORMULATION FOR SHELLS OF REVOLUTION . . . 89 H ' e 4. Interpretation of problem in terms of discrete elements .............. 89 4.2. Coordinate systems . . . . . . . . ..... 92 4. 3. Reduced set of generalized displacement parameters................. 93 4. 4. Node coordinate transformation matrix . . . 97 4. 5. Formulation of shell stiffness matrix ..... 98 4. 6. Shell imperfection . . ............. 102 4. 7. Boundary conditions and loading ....... 104 4. 8. Computer program .............. 107 V. COMPUTATIONAL RESULTS ..... . ........ 108 5.1. Some general remarks ............ 108 5. 2. Linear results ................. 109 5. 3. 'Non-linear results ............... 114 5. 4. Conclusions .................. 122 BIBLIOGRAPHY .................. . . . . . 123 iv Appendix A . Appendix B . Appendix C . Appendix D . OBLIQUE C OOR DINAT ES ........... Some properties of oblique cartesian coordinates systems ...... . . . . . . The plane stress strain energy expression referred to oblique cartesian coordinates . The thin plate bending strain energy expression referred to oblique cartesian coordinates ........... . ..... CONS TR UC TION OF TRANSF OR MATION MATRIX 13(0) ................. (i) DERIVATION DETAILS FOR MEMBER STIFFNESS MATRICES . . . . ........ Membrane stiffness matrix for triangular member. . . ................ Bending stiffness matrix for triangular member ......... . ......... SOME ADDITIONAL RESULTS ON THE TRIANGULAR PLATE BENDING STIFFNESS MATRIX ............. An alternate form of member bending stiffness matrix . . . . . . . . ...... Member transformation matrix for bending . . . . . .............. COMPUTER PROGRAM ............ Sequence of computer operations ....... List of Fortran symbols . . .......... Fortran program ............... Binary deck structure ......... . . . . Page 126 126 130 132 134 135 135 139 142 142 144 147 147 149 152 171 LIST OF TABLES Table Page 3.1 Functions used to construct displacement modes. . 54 3. 2 Membrane stiffness matrix ............ . 74 3. 3 Bending stiffness matrix . . ............ 75 4.1 S matrix ........ . ..... . ....... 99 5.1 Comparison for a uniform stress field ....... 110 5. 2 Displacements along a radial line for a concentric plate (Fig. 5.1) . . . ......... 112 C1 Some definite integrals ............... 137 D1 Alternate bending stiffness matrix ......... 146 vi LIST OF FIGURES Figure Page 2.1 Hyperplane formed by load vectors ..... . . . . 23 2. 2 Triangular member i ................ 33 3.1 Membercoordinates................. 39 3. 2 Covariant and contravariant node base vectors . . . 39 3. 3 Member parameters and oblique coordinates . . . . 41 3. 4 Unit normal and tangent vectors ........... 41 3. 5 Covariant components of node displacement relative to oblique coordinates . . . . . . . . . . 50 3. 6 Displacement modes associated with node rotations. 52 3. 7 Lines of intersection of planes given in Table 3. 1. . 54 3. 8 Symbols used to designate subdomains of triangle ....................... 54 3. 9 Graphic description of functions Z12 and Z13 . . . . 56 3.10 Force vector subjected to a'virtual displacement . . 56 3.11 Relationship of generalized forces to node variables . ..... ' .......... . ..... 63 3. 12 Relationship between displacement components normal to an edge and node rotation components . 79 3.13 Variation of normal slope for secondary bending modes.. ...... ................81 3.14 Notation for two adjacent triangles . . . . . ,. . . . 83 3. 15 Type of approximation implicit in bending modes . . 87 4.1 Shell geometry . . . . . ...... . ........ 90 4.2 Triangulation . . ......... . ......... 90 4.3 Node coordinates ........ . ....... . . . 91 vii Figure Page 4.4 Ordering of triangular members . . . . . . . . . . 95 4. 5 Distribution associated with generalized coordinatesv3.......‘........... 105 5.1 Comparison for auniform stress field . . . . . . 110 5. 2 Vertical displacements along a radial line for a concentric flat plate with fixed boundary conditions..................... 113 5. 3 Load-deflection curve for a concentric flat Plate 0 O O 0 O O O O O O O O O O O O 0 O O O 0 ll 5 5. 4 Axial load-deflection curve for a shallow cone. . . 116 5. 5 Axial load-deflection curve for a shallow cone. . . 117 5. 6 Axial load-deflection curve of the shallow cone shown in Fig. 5. 5 for various boundary conditions..................... 118 5. 7 Displacements along radial line for the cone of Fig. 5.5with'free boundary conditions . . . . 119 5. 8 Displacements along radial line for the cone of Fig. 5. 5 with hinged boundary conditions . . . . 120 5. 9 Displacements of node points along radial line for the cone of Fig. 5. 5 with the various boundary conditions of Fig. 5.6 . . . . . . . . . 121 A1. Coordinates and base vectors . .......... 127 A2. Base vectors .................... 127 viii NO TA TIONS A scalar quantity is represented by a lower case Greek or Latin letter except the letter A. A bar above a lower case letter implies the quantity is a vector (not a column matrix). A capital Greek or Latin letter is used to designate a matrix with the exception of A . The elements of matrices are scalars with the exception of J, H, Hi which have elements that are vectors and 9 which has elements that are vector operators. The transpose and inverse matrices of B are respectively designated by BT and B-1. Two groups of rectangular coordinate systems, associated with triangular members and node points, and one general coordinate system are used. A bracketed subscript refers a quantity to the appropriate coordinate system and/ or node point. Greek, Latin, and the 0 subscripts respectively refer quantities to node, triangular member, and general coordinate systems. In Chapter II we refer to triangular member 1 with node points 0., 5, 'Y3 in Chapter III we refer to a triangular member with node points 1, 2, 3; and in Chapter IV we refer to a triangular member i, j with node points ((11, BI), (a 2’ [32), (C13, 03). ix 1'2’3 E —1-2-3 e,e,e El, 82, e3 '51 .52 33 (l)’ (1)’ (1) 31(1)’ 32(1)’ 23(1 area of triangular member subsets of points of member middle surface constants matrix relating E to AV constants matrix relating A to AU matrix relating AU to AV submatrices of CR matrix relating AA to AV matrix relating J(0) to J(G) magnitude of slope discontinuity between adjacent members flexural rigidity matrix of deformation parameters contravariant, covariant base vectors of 61.. 62, Q3 oblique coordinates base vectors defined in Fig. 3. 2 (Y) ga ' gbk k 3(3) '3‘! Eif Young' 8 modulus __ T T T ‘[F(1,1>' F(Y,Y)' Fow] generalized forces associated with Ui -[f1 f2 f3 m m m ] - H.811, (Yd—'1), (Y,e1)’ 1(Y.81)' 2(y,e1)’ 3(Yael) force components are referred to base vectors 81(1)’ e2“), e3“) and momlent :gmpogients are referred to base vectors e“). ea): e(1) =[f1.,f2.,f3.,m1.,m2.,m3.] (Y1) (Y1) (Y1) (Y1) (Y1) (Y1) components referred to i-member coordinate of node point Y force vector of node point Y matrix relating stress resultants to deformation parameters k = l, 2, 3; functions used to construct Bil imperfection function column matrix with elements 171! column matrix with elements Bil specified displacement configurations of shell specified displacement configurations of member i xi (0) ‘ [ Jwr :(or J(0)] (1) [Jm' :(1)'J(1)] 1'3 (11) [J(<1rJ(a)' Jun] K m(Y)’ 1 m(Yi)' n1, n2, Am _- m F1 (Y 2 (Y1. 3 ) )“m 3 (Yi)i matrices of base vectors for O-general, i i-member, a -node coordinate systems shell stiffne s 8 matrix member stiffness matrix (includes rigid body displacements and rotations) subsets of points along side of member side lengths of member matrix relating Ai to components of node displacements and rotations matrix relating node displacements and rotations referred to different coordinate systems matrix relating A i to AUi submatrix of Mi matrix relating A i to AVi covariant, contravariant stress resultant (moment) components moment vectors of node point Y i-coordinate components of r71”) unit vectors xii :PR 1X7).1AP iflnl.n2.7) Q Ei(n1.n2.'r) ql1q21q3 r,¢,z ~ 1 2 r(n .n .7) 1.(Y) S column matrices of generalized loads generalized implied load matrix generalized residual load matrix load intensity, load intensity increment parameters surface load vector generalized load distribution matrix surface load distribution vector force components defined in Fig. 3.11b cylindrical coordinates position vector of shell middle surface position vector of node point Y matrix relating triangulation node points to member node points; also used to designate middle surface of shell same as S but for member i matrix of stress resultants thickness of member unit vectors T T T AU U (p.p)] z: [AIN1.1)’ (2,2y ... , A xiii T T T . . AU. =AU.,AU.,AU. ofmmbr1W1th 1 [ (a1) (111) mm] 8 8 nodes 0. [3 Y T 1 2 3 l 2 3 AU = A , A , A , A9 , A9 , A9 (YY) [ u(YY) u(YY) u(YY) (YY) (YY) (YY)] ua relative average axial displacement of the two ends of a shell of revolution Aul . ,Au2 . ,Au3 . i-member coordinate components of Au (Y 1) (Y 1) (Y 1) (Y) Au1(YY)'AuZ(YY)’Au3(YY) Y -node coordinate components of A11”) .. 1 2 - . . u (n , n ,7), Au shell middle surface displacement vectors Bi, Afii member i middle surface diaplacement vectors EL. displacements relative to S; 1 GR rigid body displacements of S; i 1.1”), A11”) displacement vectors of node point Y V column matrix of generalized coordinates Vi a subset of V VI generalized coordinates associated with generalized displacements h! We’ We. work of surface and boundary loads on shell; 1 on member i ' WI’ WI strain energy of shell, of member i i WI’ ,WI membrane, bending strain energy 1' 2 xiv wl, w2, w3 w1(1)' W2(1)' “’30) x1 x2 x3 r Z (Y)’ (Y)’ (Y)’ (Y)’ (Y) XXX 1’ 2’ 3’ 1 2 3 Y 1 Y 1 Y 1 2 3 YR! YR! YR 6(...) 61, ..., 66 R,Z covariant components of T1 referred to 12 3 oblique coordinates L1, Q 1 Q covariant components of A11( referred to base 1) vectors 31 32 E3 (1)’ (1), (1) O—general rectangular coodinate system, fixed in space components of T”) column matrices of the above components i-member rectangular coordinates, redefined after each linear increment rectangular coordinates of member i, member rigid body rotations (Sec. 3. 5) are zero relative to them a ~node rectangular coordinates 1 2 = 12,1. 11 £2 diagonal matrix with elements Pi member stiffness matrix a small finite change in the internal A'rk a virtual change displacement increment constant middle surface deformation parameters XV 611’ 622' £12 1511’ 322' $12 1 (YY)’ 2 A9 A9 1: (YY)’A (YY) covariant strain components of L1, L2, £3 coordinates covariant strain components in rectangular coordinates curvilinear coordinates of a surface Y -node coordinate components of AH”) rotation vector of node point Y covariant curvature components T T =[A1, .Am] =[11. .112] member 1 generalized coordinates = 2(1 - v) Poisson's ratio oblique coordinates fixed to member 1 matrix relating E to A1 matrix with elements 0'1 generalized forces of member 1 xvi covariant, contravariant stress components parameter associated with various deformed shell configurations rotation vector of node point Y measured relative to coordinates yllR’ vi. 3'; matrix operator residual load relaxation constant xvii I. INTR ODUC TION l. 1. Preliminary Remarks In engineering, particularly in the aerospace field where weight minimization is of the utmost importance, thin walled structures, consisting of thin bars, plates and shells, are widely used. To facilitate the formulation and solution of broad classes of these structural analysis problems and related dynamic and aeroelastic problems, discrete element methods using matrix notation are now widely used. Their increasing use has very closely followed the improvements and increasing availability of digital computers. These methods are part of the so called ”matrix methods of structural analysis" put forward by Langeforsl and Argyris. The components of these structures are characterized by their flexibility, i. e. , their relatively small resistance to bending and torsion. When they are loaded the resulting displacements are frequently comparable in magnitude to their linear dimensions. The classical (linear) theory of elasticity and in particular the theory of shells is based on the assumption that the displacements of points in the body are infinitesimal which in turn permits a formulation of equilibrium and compatibility conditions with respect to the undeformed geometry. If a theory requires the formulation of the conditions for equilibrium and compatibility to be with respect to the deformed geometry, as in reality it is, then it is said to be "geometrically non-linear. " If the constitutive equations are non-linear, then the theory is said to be "physically non-linear. " In the following, the Kirchhoff-Love hypothesis is assumed. Some books on shell theory used as source material are given by references 3, 4, and 5. The discrete model method used requires that the shell geometry be modeled by flat triangular plate members. Each triangle is assigned 12 independent deformed configurations specified so that adjacent member displacement continuity and in a limited sense (see Sec. 3. 7) normal sloPe continuity is maintained between triangular members by maintaining these conditions at corner points. These are used to obtain triangular plate stiffness matrices in the generalized displacement sense which in turn are used to construct the stiffness matrix of shell structures. The modeling is discussed in detail in Chapters 11 and HI. This investigation is devoted to the formulation of arbitrary plate flat triangular‘stiffness matrices associated with the membrane and bending stress states, and their applications to geometrically non-linear thin elastic shell problems. The investigation is restricted to shell materials which are isotrOpic,homogeneous, and linear elastic, and numerical results 'are presented for axially symmetric geometries. The method is applicable to somewhat general configurations and boundary conditions and can be directly extended to include physically non-linear materials. As formulated here, this method may be interpreted as consisting of two parts, i. e. , the linear problem and the linear incremented extension of the linear problem into the non-linear range. The linear problem is closely related to the Rayleigh-Ritz method; however, the concept of a minimizing sequence cannot be interpreted directly, at least not for the more general problem considered. Each element of such a sequence would be associated with a different triangulation. The inconsistency of the require- ment is discussed later. The method is interpreted as a direct method of the calculus of variations which gives an approximation to the problem, i. e. , a set of arbitrary constants (generalized coordinates) associated with a set of admissible shell displacement modes (generalized diSplace- ments) are determined so that the integral of potential energy is minimized. An admissible shell displacement mode is interpreted here as one which forms a compatible field, 8 i. e. , satisfies compatibility conditions and displacement boundary conditions. In this case the potential energy is bounded from above. 8’ 9 A compatible field is not necessarily an equilibrium field, i. e. , it does not necessarily satisfy the equilibrium equations and stress boundary conditions. If a solution yields a compatible field and also an equilibrium field then it is said to be exact. In its relationship to the matrix methods of structural analysis the method us ed is part of the matrix displacement method. 2 This is also called the direct stiffness method. The non-linear problem is essentially a step by step procedure based on the linear formulation. The iteration can be interpreted as consisting of two parts; i. e. , the load is advanced in increments (since the solution may encounter unstable regions the load can increase or decrease in increments) and it searches for a deformed configuration in equilibrium with the specified load. Equilibrium is here interpreted to mean with reSpect to the deformed configuration. 1. 2. Previous Deve10pments Much of the impetus for the development of discrete model formulations to shell problems, at least during the past two decades, has come from the aircraft industry and in particular from applications to dynamic and aeroelastic problems. The various discretization procedures used to approximate the behavior of a structure can be classified according to their properties of either satisfying compatibility but not equilibrium, or satifying equilibrium but not compatibility, or violating both equilibrium and compatibility. It is desirable that a procedure admit to a refinement which in the limit converges to the exact solution and/ or converges monotonically. Hrennikoff10 developed and McHenry11 improved on the ”frame work analogy method" in which an analogy, consisting of a beam element lattice, is made to the plane stress problem. It was later gener- alized to include bending by Parikh and Norris. 12 This method implicitly relaxes both compatibility and equilibrium conditions so that solutions based on it do not in general form either compatible or equilibrium fields. Many procedures similar to this have been and are now being used in both static and dynamic applications. We note, in this connection, that some finite difference formulations to some differential equations in elasticity implicitly relax both equilibrium conditions and compatibility conditions, and satisfy these conditions in general only in the limit as the mesh size is made smaller. A plane stress triangular plate element stiffness matrix, in the generalized displacement sense, was put forward in a paper by Turner, Clough, Martin, and TOpp. 13 The stiffness matrix is associated with three independent deformed configurations which have displacements that vary linearly in all directions and strains that are consequently uniform over the entire element. This matrix is now widely used. It was given a different form by Argyris. one which he calls the natural or invariant form. In the absence of bending, solutions based on this matrix give displacements which form a compatible field for the triangulated model of plane or curved surfaces and for linear or non-linear problems. If bending is present, then compatible fields in general are obtained only for plane surfaces and geometrically linear problems. In general this matrix does not yield equilibrium fields; how- ever, de Veubeke8 gives an alternate approach which yields equilibrium but not compatible fields and uses it, in what he calls a dual treatment, to obtain upper and lower bounds to static influence coefficients. In order to satisfy compatibility between a triangular plate and a beam segment de Veubeke15 generalizes the plane stress stiffness matrix to include parabolic variations in displacements along its sides. The formulation requires nine independent deformed configurations for each triangle and that interelement displacements be satisfied at corner points and at the midpoint of sides. In their paper Turner et. all3 also present a membrane stress rectangular plate stiffness matrix. As pointed out by Melosh16 this matrix does not in general yield a compatible field. Compatibility conditions are satisfied in the interior of elementsand at node points; however, gapping may result, 1. e. , displacement compatibility between elements is not necessarily maintained. In the same paper Melosh presents a rectangular stiffness matrix which yields compatible fields. This matrix is associated with five independent deformed configurations. Argyrisl4 gives a presentation of the so-called natural forms of these nodes. All the deformed configurations have displacements that vary linearly along the sides of the member, however, the diSplacements associated with them have terms that are quadratic. Melosh gives a sufficiency condition for increasing the number of rectangles so that with each refinement the potential energy monotonically approaches a minimum. In essence he requires a refined subdivision, in a sequence of subdivision, to be so constructed that its displacement field can contain any displacement field of a coarser subdivision. This sufficiency condition can be interpreted with respect to the plane stress triangular plate elements. Displacement -modes for flat surface plane stress problems obtained by a sequence of subdivisions satisfying Melosh's sufficiency condition can be used to obtain a minimizing sequence similar to the Rayleigh-Ritz type; however, the concept of completeness necessary for convergence to the exact solution has not been demonstrated and would, undoubtedly, require additional sufficiency conditions, as noted by de Veubeke. 15 ArgyrisI4 presents a paralleIOgram plane stiffness matrix and indicates a method for constructing one for a plane quadrilateral panel which he has obtained. A11 yield compatible fields. In a report by Bogner, Mallett, Minich, and Schmitl7 the authors give the displacement modes for constructing the stiffness matrix for a curvilinear rectangle associated with any orthOgonal curvilinear coordinate system. For rectangular coordinates it reduces to those of Melosh. Melosh16 gives a rectangular bending stiffness matrix. As pointed out by Pian31 the displacement modes do not in general maintain 810pe compatibility between elements and consequently solutions based on this matrix do not yield compatible fields. Melosh18 previously had given a bending stiffness matrix. Bogner et. a1. 17 present displacement modes which can be used to construct bending stiffness matrices for curvilinear rectangles, which yield compatible fields and can be used with Melosh's sufficiency condition to obtain monotonically converging sequences. These rectangular member displacement modes and those of Melosh do not form an independent set of deformed configurations in as much as they include rigid body displacements. 19 Clough, 15 Adini, and Zienkiewicz15 use a polynomial to numerically calculate the coefficients of a rectangular plate bending stiffness matrix, i. e. , the polynomial forms a set of displacement modes which in this case include rigid body displacements. Zienkiewicz15 extends this to quadrilateral plates. Clough15 and TocherzO use similar polynomials for triangular plates. Several objections can be raised with this type of procedures as pointed out by the authors. With the exception of the rectangle the use of these polynomials in general results in displacement discontinuities between members and in all cases s10pe continuity is not maintained between members. A member, and in particular a triangular member, will in general have different stiffness coefficents depending on its orientation with respect to the coordinate system of the polynomial even after they have been properly transformed so that components of node rotations and displacements are with respect to the same coordinate system, i. e. , the stiffness coefficients are not uniquely defined by the polynomial. This points out the desirability in selecting deformed configurations that reflect the geometric prOperties and symmetries of the member and which are independent of pure rigid body displace- ments as Argyris did in his plane stress formulations. The matrix methods of structural analysis were originally developed to facilitate linear formulations of complex problems. These methods of analysis have been given three classifications: the displacement formulation, the force formulation and the combined formulation. They have been shown to be equivalent to stationary energy principles, i. e. , the displacement formulation is equivalent to the principle of stationary potential energy, the force formulation is equivalent to the principle of stationary complementary energy, and the combined formulation is equivalent to Reissner's principle 21’ 22 A useful survey on linear structural of stationary energy. . . . . 23 analyS1s 13 given by Argyris.. During the past few years various investigators have extended these linear procedures to include geometric non-linearities, and . . . . . . , 7,24-2 both conservat1ve and non-conservat1ve material non--11near1t1es.l4 l 8 1. 3. Present Investigation This investigation was primarily directed at obtaining an plate explicit representation for an arbitrary flat triangular‘bending stiffness matrix and to studying its applicability to geometrically linear and non-linear shell problems. Two bending matrices were obtained. The form in which these matrices are used requires a representation of node variables with respect to rectangular coordinates. A direct representation in this form is, however, so awkward that it has little value. The matrix is represented in the form T M2i F2 M2 (1.1) i i 10 where F21 (Table 3. 2 or Table D1) is the bending stiffness matrix with reSpect to a set of generalized variables and M21 (D. 6) is a transformation matrix that relates the generalized variables to a set of node variables referred to rectangular coordinates. The generalized variables were selected so as to take advantage of the symmetric and geometric properties of an arbitrary triangle and consequently the coefficients of both I72. and M2 have a 1 1 particularly simple form. A derivation obtained by representing all quantities with resPect to a rectangular coordinate system was found to be for all practical purposes prohibitive due to the very large volume of algebra required. By using oblique coordinates and representing the bending strain energy expression in terms of a set of deformation parameters that reflect the geometric pr0perties of an arbitrary triangle, the derivation was performed with relatively little algebra. Since membrane behavior is in general present, the plane stress matrix of Turner et. a1.13 is used. It is, however, generalized by including three additional inplane diSplacement modes associated with components of node rotation normal to the plane of the triangle. These were included in order to remove the possibility of obtaining singular shell stiffness matrices. 11. GENERAL FORMULATION 2.1. Some General Pr0perties Let 1:071, 172, ’T) be the position vector of points on the middle surface of a thin elastic shell possessing a positive-definite strain- energy function, quadratic in the components of strain. The parameters 771, n2 are the curvilinear coordinates of the surface and the parameter 'r is associated with various deformed config- urations. When the shell undergoes a continuous deformation from configuration '1:(nl,n2,'rl) to configuration ?(n1,n2,72), the parameter 'r varies continuously over an internal (7'1, 7'2). In particular the configuration of the undeformed middle surface is designated by ;(n1, n2, 0). The displacement vector of a point 171, n2 with respect to the interval (71, T2) is defined by — 1 2 _ — l 2 — 1 2 u(n.n.'rZ-'rl)—r(n.n.'rZ)-r(n.n.71) (2.1) . . — l 2 . . The displacement vector relat1ve to r (n , n , 0) 15 de31gnated by _ 1 2 _ 1 2 — 1 2 u(n.n.'r)=r(n.n.T)-r(n.n.0) (2.2) The external force per unit area acting on the shell surface is designated by Emlmzur) = p('r)?1(n1.n2.'r) (2.3) where p('r) is called the load intensity parameter and aml, 172,7) is called the load distribution vector. The load distribution vector ha 8 the pr Ope rty 11 12 ff 5 ° 13 dA = constant (2.4) S where the integration is carried out over the entire middle surface area S. If the interval (0, '1') is partitioned into n subintervals (0.3% 01.72% '°°’(Tkal’TkL ""(Tn11’7) (2. 5) = A71,ATZ, ..., A'rk, ..., A'rn then the displacement vector can be expressed in the form _ 1 Z — l 2 — l 2 um1n.fl=um.n.An)+u.+um.n.Am)+ _ 1 2 (2.6) + ... +u(n 2n 1A7) 11 In general it is assumed that < < < < < < < 0 7'1 7'2 7k “In-1 'r . (2.7) This investigation is concerned with shell deflections of sufficient magnitude so as to require a formulation of the equilibrium and compatibility conditions in the deformed state. The resulting non-linearities are dealt with by a formulation consisting of a sequence of linear intervals. Each linear interval has a one to one correSpondence with an element of the sequence {Aqi‘} . The linear problem intrinsic in the non-linear problem is thus one of starting with a deformed configuration TM], 772, 'rk) and seeking a deformed configuration ?(nl , nzxr k+1). 13 The following simplified notation is used when convenient: HI I: HI :5 p—n :5 ”-4 WV Cil ll El :5 H :1 .1 WV 5 = Emlmznrk) (2.8) ;' =r+A1-1 1_1' =T1+Au 15" = 13+A1-3 All symbols preceded by A are interpreted as finite functional changes in the interval A‘Tk . It is assumed that Au together with its first partial derivatives with respect to n1 and n2 are sufficiently small, in accordance with linear shell theory, for all partition intervals A'rk . The assumption of smallness in the interval ATk implies that geometry changes in the interval are small, and that superposition of displacement configurations and corre3ponding surface loads associated with the interval is admissible. The middle surface stress resultants and deformation parameters defined in linear shell theory are respectively designated by the six element column matrices 14 T . 1 2 6 1 2 T = [t1(n.n.'rk).....t(n.n.'rk)] (2.9) T _ 1 2 1 2 E — [61(1) .n .Tk). 66(1) .n .Tk)] (2.10) They are related by T = GE (2.11) where G is a 6x6 symmetric positive definite matrix with elements that are assumed constant throughout the range of '1' . We define AT, AE, T', and E' in accordance with (2. 8). Then AT = GAE (2.12) E' = E+AE (2.13) T' = T+AT (2.14) The strain energy, associated with the infinitesimal middle surface area dA, and strains E is given by _ 1 2 _ 1_ T dWI —dWI(r) ,n ,Tk) - 2T EdA (2.15) Substituting from (2.11) into (2.15) we obtain dWI = ~12—ET G E dA (2.16) The total strain energy is p—a T lezjéfE GEdA (2.17) The total strain energy after the next increment is similarly related to E' by .—n T W'Izzfst' GE'dA (2.18) 15 The change in strain energy WI due to a virtual change 6E is 6W = l—ff(ET+6ET)G(E+6E)dA-l—ffETGEdA 1 2 5 2S = %ff{ETG6E+6ETGE+6ETG 6E} dA (2.19) 5 On dropping higher order terms and noting that G is Symmetric (2.19) reduces to ff ETGéE dA s 6WI (2. 20) ffTTéEdA 5 Similarly the virtual change in the total strain energy after the next increment is related by ff E'TGéE dA s 1 6W1 ll ff {ETG6E+AETG6E} dA S (2.21) ff {TT 6E+AETG6E}dA s (SWI + ff AETG 6E dA S The external work due to the surface load 13 and the virtual displacement 61-1 is awe = fsf'fi- 65 dA (2.22) Similarly for 13' we obtained 6W' 2 f f '5' - 6T1 dA e s (2.23) 16 = 6W +ffAfJ- 51311.1. e S The principle of virtual work requires 0 = 6WI - 6We (2.24) and 0: 6W' - 6W' 1 e (2.25) = 6WI+ ffAETG 6EdA- 5w -fpr-6T1dA s e s From (2. 24) and (2. 25) we obtain {fpr-éudA-ffAETGéEdA} =0 (2.26) s 3 Let Q designate a six element matrix operator that relates the deformation parameters to the displacements such that AB = GAG (2.27) This Operator is obtainable from linear shell theory and in the form used is given by (3. 27). _ l 2 — 1 2 If a1, a2 are scalar constants and ul(n , r) ), u2(n ,n ) are vector functions then the linear property of this Operator requires {Mala +a 11 1 2 2) = :11 521-11 + 212 $2132 (2,23) 2. 2. Discrete Method The term discretization is used here to imply the reduction of the problem from a formulation in terms of unknown functional quantities 17 whose domain contains all points of the middle surface to a formulation in terms of a finite discrete set of unknown parameters. If the discrete set contains n unknown parameters then the discretization is said to have n degrees of freedom. The discretization is based on Rayleigh's method, in which the displacement vector is approximated by a linear combination of independent displacement configurations that satisfy compatibility conditions and displacement boundary conditions. Implicit in an 11 degree of freedom discretization is the existence of an n element set of independent displacement configurations designated in matrix formby T ‘— 1 2— .. 2— H =[h1(n.n.r). hn(n1.n.r)] (2.29) The construction of the vector functions hi depends on a knowledge of 1:011, 172, 7k) and as used here the Hi are functionally dependent on Tml, n2, 7' In accordance with the step-wise k)' linearization discussed above the functions hi are assumed to be independent of T, in the interval AT and corrected after each such interval. Each function hi is Specified so as to satisfy diSplacement boundary conditions and compatibility conditions. The displacements are then approximated by - T Au = H AV (2. 30) where AVT = [Av], Avn] (2.31) 18 are scalar parameters associated with the interval AT and called the generalized coordinates. Substituting (2. 30) into (2. 27) we obtain AE = QHT AV A 6xn matrix is defined by B = nHT Then (2. 32) has the form AE : BAV The virtual changes 61-1 and 6E are then related to 6V by 65: HT6V 6E: B6V Let pT — T =ffp H M s T szAp H dA s T — HT Q =1;qu dA The dot product is used since the elements of H are vectors. (2.32) (2.33) (2.34) (2.35) (2.36) (2.37) (2.38) (2.39) The P, AP, and Q are respectively called the generalized load matrix, the generalized load increment matrix, and the generalized load distribution matrix. On substituting (2. 34), (2. 35), and (2. 36) into (2. 26) we obtain 19 T T {fpr-H 6VdA-ffAVTB GB6VdA} = o (2.40) s s Noting that 6V and AV are not functions of 171,172 and substituting from (2. 38) into (2. 40) we obtain {APT - AVTffBTGBdA) 6V = o (2.41) 5 Eq. (2. 41) can be satisfied for all virtual changes 6V only if APT-AVTffBTGBdA = o (2.42) s The matrix K, called the shell stiffness matrix, is defined by K = ffBTGBdA (2143) S Since G is symmetric and positive definitive, K is symmetric and positive semi -definite and if in addition det K > 0 then K is symmetric and positive definite. In the physical sense the stiffness matrix is positive semi- definite if the . elastic System is not tied down, i. e. , if an adequate number of constraints have not been imposed so as to prevent rigid body displacements or rotations of the entire structure. Consequently, the displacements associated with a given loading are not unique and the relating coefficient matrix K is not invertable. Substituting (2. 43) into (2. 42) and taking the transpose of the entire expression we obtain AP = K AV (2. 44) 20 The generalized displacement parameters are chosen so that the displacement boundary conditions may be Specified by specifying a subset of these parameters AVZ. The matrices in (2. 44) are partitioned as follows: = (2. 45) where AVl are unspecified and AV are Specified generalized 2 coordinates, and AP are specified and AP2 are unspecified 1 generalized forces. The matrices API and AP2 are resPectively called the generalized load increment and generalized reactions increment. If in (2.45) the det Kll > 0 then K11 is invertable and we can solve for the unknown generalized displacements and reactions. We then obtain -1 .AVI = K11 {APl - KIZAVZ} (2.46) _ -1 APZ _ K21K11 {APl — Klev2)+ ansv2 (2.47) If a load vector ”15 = pq is in equilibrium with a distribution of stress resultants T , then from (2. 20), (2. 22) and (2. 24) it follows that ffp-éfidAszTTéEdA (2.48) s s 21 Substituting from (2. 35) and (2. 36) into (2. 48), noting that the resulting expression must hold for all 6V, and then substituting from (2. 37) we obtain P=ffTTBdA (2.49) S The above is a relationship between the externally applied generalized load P and internal distribution of stress resultant T for an elastic system in equilibrium, 1. e. , it is a form of the equilibrium equation. In the sequel we use the equation relating the load intensity parameter p(7'k) to the generalized load matrix P and generalized load intensity matrix Q. From (2. 3), (2. 37), and (2. 39) it follows that P : pQ (2.50) QTP = pQTQ (2.51) T p = Q—T—P (2.52) Q Q The matrix H, as already indicated, is functionally dependent on '1': , and consequently B and K are also. The distribution of stress resultants T is associated with configuration T = T071, 112, TR), but it is determined from quantities defined with resPect to configuration ;(nl, n2,7k_1). Consequently the generalized forces implicit in T do not necessarily conform to the specified load distribution of Q. The matrix PI’ associated with 'r = 7' and k called the implied generalized load matrix, is related to T by (2. 49) and is given by 22 pI = ffTTBdA (2.53) S The load intensity p = p(7'k) is defined by QT P1 p = T (Z. 54) Q Q The matrix PR , associated with 7k and called the residual generalized load is defined by P : P - pQ (2. 55) The above two definitions (2. 54) and (2. 55) have a conceptually useful interpretation. We multiply the terms of (2. 55) by pQT and obtain 11er PR = pQT P1 - pZQTQ (2. 56) Solving for QTPI in (2. 54), substituting into (2. 56) and dividing the resulting expression through by the scalar p we obtain Q P = o (2.57) It then follows that Q is orthogonal to PR' If we interpret Pi, PR’ and Q geometrically as vectors in a hypersPace,then these three vectors may be interpreted as lying on a hyperplane, pictorially shown in Fig. 2.1, and the vector PI may be interpreted as being equal to the sum of two orthogonal vectors PR and p0. 23 hyperplane Fig. 2.1. Hyperplane formed by load vectors. It is desirable to choose a generalized load increment APl so as to minimize the next residual load. There does not appear to be any convenient method of making a best choice for APl. It is taken in the form AP = Ale - wP 1 (2. 58) R1 where Ap is called the load intensity increment and w is called the residual load relaxation constant. The matrices Q1 and PR in l (2. 58) are submatrices of Q and PR obtained by partitioning in accordance with (2. 45). Substituting (2. 58) into (2. 46) we obtain AV = K;1{Ale-wP 1 - KlZAVZ } (2.59) R1 24 Note that 'i-i = I- + All (2.60) .. AV =r+[H1H2] 1 AV2 The matrices H, B, K are constructed after each linear increment from a knowledge of T . This is symbolically represented by HG), 13(5), KG). The algorithm employed in determining the elements of the sequence { T011, 112, 7k)} is obtained by substituting (2. 59) into (2. 60) and is given by $1 = ¥ +[ H163) 112(5)] Kl‘i(?){Ale - wPRl - K12(?)AV2} ...sz '...) (2.61) Associated with every deformed configuration T071, 172, 'rk) we have a residual generalized load matrix PR(7k). If PR(7k) = 0 then - 1 r(n , flank) is a configuration in equilibrium with the specified loading. In general, elements of { Tml, n2,7k)} do not satisfy this condition and consequently their acceptability is determined by interpreting the smallness of the magnitude of PR(7k). 2. 3. Simplified Model To simplify the analysis of this problem the shell is interpreted as a model consisting of flat triangular plate elements. This process can be interpreted as a p node In element triangulation of T011, 172, 0) into curvilinear triangles and an isomorphism which deforms the 25 curvilinear triangles into simple triangles without materially disturbing the relative position of node points. The assumptions implicit in the simplified model cannot always be justified. This restriction is not, however, very severe for many applications. In the following, quantities associated with node points have Greek subscripts, quantities associated with triangular elements have Latin subscripts and quantities associated with the shell in general either have no subscripts or the subscript 0. 2. 4. Matrix Relating Nodes and Members The formulation of the shell problem in terms of discrete parameters defined with respect to the simplified model described above deals with configurations of node points and triangular members formed by lines joining node points. This leads to the definition of the mathematical sets N and S . The node set N is defined as the set of all node points of the simplified model and the member set S is defined as the set of members Si or of ordered node triplets (o. , (3, Y) that in turn define all triangular members of the simplified model. A correspondence between elements Si and ordered triplets (a, (3, Y) is formed and represented in matrix form by S. = [a1 (3. Y] (2162) 1 The dependence of S on N is symbolically represented by S : S(N) (2.63) 26 Elements of S and N are ordered respectively by integer numbers (2.64) l: 2: 00-, Cl, ..., p The functional dependence implicit in (2. 63) is explicitly represented by an mx3 matrix S of integer elements and is F5 5 s S r— ‘- 11 12 13 ° ' s =[sij] = = (2.65) 811 312 513 a ‘3 Y The Symbol S will serve the dual but analogous role and will also refer to the set of all points on the middle surface of the simplified model, and similarly Si will refer to the subset of S associated with the ith triangular. member. The quantities T, 11, and A171 are defined over S and the quantities ;i’ iii, and Ai-ii are defined over Si . 2. 5. Interpretation of Non-Linear Problem The member Si is displaced and deformed from some initial undeformed position.. Part of the rigid body displacements are described by a triangular member 3' associated with Si (the node points of Si are designated by (1', (3', Y'). The triangle Si is identical to the triangle Si in its undeformed state and is fixed to the deformed member Si so that point 8', line B'Y ', and plane a 'fl'Y' lie respectively on point (‘3, line BY and plane QBY. 27 The displacements of points on Si can be represented in the form Au. 2 Au . + Au . (2.66) 111 = uR + uL. (2.67) where A13. , T1 are displacements of S! and A11 , T1 are R1 R L. i 1 1 Li displacements of Si relative to S; . Since AER- and ER- are rigid body displacements, the 1 1 magnitude of the strain deformation parameters correSponding to them are identically zero. From (2. 27) we obtain (2. 68) = QAfiLi The displacements uL_ are approximated by 1 — _ T uLi - HiA i (2.69 and consequently AB = HTAA (2 70) L. i i ° 1 where Hi is a column matrix of linearly independent vector valued functions and A i is a column matrix of scalar parameters called the member generalized coordinates. Eq's (2. 69) and (2. 70) imply that Ai(7k) = AAi(A7k) + + AAi(A'rk) (2.71) 28 In a manner similar to that used to obtain (2.33 ), (2. 34), (2. 38), (2.43), and (2. 44) we obtain II. = 52H: 1 1 AE: IIAA. 1 1 Azi = ijAp HldA 1 (2.72) (2.73) (2.74) (2.75) (2. 76) where 21 is called the member generalized force matrix and Pi is called the member stiffness matrix. The matrix Pi is symmetric and positive-definite. From (2. 71) it follows that 21: Pi‘Ai Let AT _[J\1T 23T : [2T, 1 II = [III, I" 1" z o (2. 77) (2.78) (2. 79) (2. 80) (2. 81) 29 Then AB = TIA (2.82) 23 = FA (2. 83) The total strain energy in the shell for 7' = 7k is wI =JZ—2TA =% ATI‘A (2.84) The displacement of a point measured relative to the undeformed shell is interpreted as the sum of two parts. One part is associated with certain rigid body displacements and rotations of the entire triangle and the second part is measured relative to these displace- ments. The linear relation (2. 83) is associated with displacements obtained after the first part has been subtracted off, and is assumed valid through the range of 7' . The non-linearities are, therefore, limited to those resulting from displacements associated with rigid body displacements and rotations of triangles Si . 2. 6. Shell Stiffness Matrix The shell generalized coordinates V and generalizedforces P are related to Z and A of the same elastic system by the strain energy expres sion 8wI = awe = PT6V = zTaA (2.85) The two sets of generalized coordinates, in accordance with the step- Vvise linearization used, are related by 30 A A = CVAV (2.86) and consequently 8A = cvav (2.87 In general the elements of CV change with 7' so that CV(7k) ,1 CV(Tk+1) (2.88) Substituting (2. 87) into (2. 85) we obtain PT6V = zchav (2.89) Eq. (2. 89) can hold for all virtual displacements 6V only if P = c 2: (2.90) T V From (2.34), (2. 82), and (2. 86) we obtain AEleAA = IICVAV = BAV (2.91) and c ons equently B = IIC (2.92) Substituting (2. 92) into (2. 43) and noting that the elements of CV are not functions of 111, n2 we obtain KszBTGBdA =ffCTV'IITG‘II'CVdA s s (2.93) _ T T _ T T _CV{fo1'I GTIdA}CV—CV[SfoIIiGTIj (1.1).]cV ij 31 The integration in the last of the integrals of (2. 93) is carried out over points common to both Si and Sj since “IIi is defined over Si and IIj is defined over Sj' Since the triangles do not overlap and from (2. 75) is follows that o if 1,1 j ffnfon. dA = (2.94) 85 J P. 111:1“ i j 1 Then from (2. 81), (2. 93), and (2. 94) we obtain _ T K — CV l"CV (2.95) 2. 7. Member - Shell Parameter Transformation Matrix If the generalized coordinates are interpreted as components of node displacement and rotation vectors then the associated components of generalized forces, being related through the strain energy as in (2. 85), can be interpreted as components of force and moment vectors. The shell generalized coordinates and generalized forces having this interpretation are respectively designated by the column matrices U and F. The formulation in terms of U and F is desirable since most boundary conditions and in a limited sense, load distributions can be irterpreted directly in terms of these parameters. The undeformed simplified model has a continuous middle surface consisting of flat triangular members. Two adjacent undeformed triangles have an intersection line and intersection angle. We have required the displacement vector A11 to satisfy compatibility conditions on the entire middle surface. This can 32 be satisfied only if continuity and intersection angles are preserved for all points on the intersection line. These requirements are not strictly satisfied nor essential. The requirements are useful, however, inasmuch as the approximating surfaces that satisfy them also satisfy compatibility conditions with respect to the simplified model. In this case the elements of the shell stiffness matrix can be interpreted as satisfying certain upperbound requirements. 9 It is desirable to choose the functions Hi so that whenever compatibility requirements between two adjacent triangles are satisfied at their common node points they are satisfied at all points on their intersection line. The degree to which this condition is satisfied is discussed in Sec. 3. 7. The member displacements obtained by substituting (2. 70) into (2. 66) are AB. = AB + H? AA. (2.96) 1 R1 1 1 These displacements can be directly related to AU. In the formulation used the member displacement vector (2. 96) is related only to components of AU associated with its nodes. Each node has six elements of AU associated with it and consequently each member has 18 degrees of freedom. The rigid body displacements AER- 1 have six degrees of freedom associated with them, i. e. , three translation components and three rotation components. In addition AAi has degrees of freedom associated with it equal in number to its dimension. Since the member has 18 degrees of freedom and six can be interpreted as rigid body displacements it follows that the 33 dimension of AA must be 12 or reducible to 12. That is, the formulation used requires that all possible displaced configurations Aui be linear combinations of six independent rigid body displace- ments AER- and 12 independent deformed configurations H?AAi . 1 Consider a triangular member 1 with node points a, B, Y and a reference point 0 as shown in Fig. 2. 2. We associate a rectangular cartesian coordinate system r(c1) Fig. 2. 2. Triangular member 1. with each of the quantities i, 0, a, (3, Y and refer to them as the i-member, O-general, 0. -node, 8-node, and Y -node coordinate systems. The base vectors of these coordinate systems have the £0 m . , . , . and th mat of th 8 base vectors 18 r 3(1) 3(1) 3(1) e r1x e e designated by T *2 “'3 ] (2.97) -:-1 J I = a ’ c , I (1) [Jm J(1) J<1> The bracketed subscripts on a quantity refer it to a node point and/ or a coordinate System. These subscripts are used only when necessary for clarity. 34 The matrices J and I“) are related bya 3x3 matrix (0) 131(0)) where the indices are interpreted by the relations .1 : D“) .1. (2.98) -1 T Dlm) = Dw) (2.99) m m 1m (i) 's definedb (0) 1 Y r . ' D1“) 0 m _ m D(0) — O D (i) (2.100) 1m The matrix D The position vector of node point a relative to point 0 is Tm) and its O-general coordinate components are xléa), X(Zo.)’ x(3a) . The matrices of their components are respectively designated by X1, X2, X3. The elements of these matrices are arranged in accordance with the ordering in (2. 64). The i-member coordinate system is redefined after each interval A 7k so that its origin is incident to point a , base vector 3%) is normal to plane a (3 Y and base vector It) lies along line a 8 and directed toward (3. The base vectors J“) can be represented in terms of the position vectors 1'0, rfi, ;Y as follows fl): flm-RQ_ “m'fiwl 3.3. _ {(1711) ' I'(11))X(rhtl ' rm)” (2.101) (1) |('r' _; )xG -? )l ((3) (a) (Y) (9) T2 : ':'3 X31 Jm (0 35 By representing the position vectors in (2. 101) in terms of the elements (0) mm of Jw), the elements of D can be obtained directly and are given in Appendix B. The base vectors J(a) of the 0. -node coordinate system as used here do not change with 7' and are specified so that the (0) transformation matrix 131(0) can be constructed directly. Then we obtain (a) _ (0) (11) DH” _ Dunfhm) an0m Node point a has a displacement vector A11”) and rotation vector 135(0). The O-general coordinate components of 2 2 - " 1 3 1 0 . Aum) and A9(a) are Au(a0)’ Au(a0)’ Au(0.0) and A9(c1 of),A ((10)’ A9?C1 0.). The a -node generalized diSplacement matrix is T _ 1 2 3 1 2 3 AU(a 0) ‘ [2% or Au(nor Au(<10)’ A9(111w AG(<1 or A9010)] (2.103) The generalized displacement matrix is T T T T AU = AU , AU , ..., AU 2.104 [ (1.1) (2, 2) (11.11)] ( ’ The column matrix AUi is defined by T T T T AU. = AU ., AU . , AU . 2.105 1 [ (<11) ((11) 0/1)] ‘ ’ We can construct a transformation matrix Mi defined by AA.: M AU. (zuom 1 1 1 36 We partition Mi in the form p—- —1 M13] AU(ai) (2.107) AU(131) AAi : [M11 M12 LAUWi) — The matrix Mi is explicitly given by (3. 74). The following relationship follows from (2.102) AU _ 11‘“) o o [ rAU — (111) (i) ( ) (ca) _ (3 AUmi) _ o D“) c: ) AUmB) (2.108) Y AU< Y1) O O 13(1) AU(YY) ... L. _ _ .3 Substituting (2.108) into (2.107) we obtain _ (a) ((1) (v) ‘ ‘ AAi ‘[M11D(1) , Mi2 11(1), Mi3 13(1) ] AU(aa) (2.109) ”(1111) AU _ (YY)_J We define the transformation matrix C, which is similar to CV, as follows AA = CAU We partition C into 6x6 submatrices in the form 011 C12 cipj c - = [C (2.110) jk] m1... mp 37 The algorithm for constructing C can be obtained by examining matrix S defined in (2.65) and (2.109) given above. Then for jzl, ..., m and k=l, ..., p Cjk = 0 exceptfor CJ11 - M31 13(3)) where n = ij III. TRIANGULAR MEMBER MATRICES 3. 1. Member Coordinate Systems Consider a triangular member with middle surface points referred to rectangular coordinates yl, yZ, y3 and y;, yé, y; and oblique coordinates £1, £2, £3. The member is displaced from an initial position as shown in Fig. 3.1a to a position as shown in Fig. 3. lb. The coordinates £1, £2, 23 are fixed to the member as symbolically shown in Fig. 3.1b. The coordinates Viv YR’ y; are fixed so that the components of rigid body rotation of the triangle (defined explicitly at the end of Sec. 3. 5) relative to them vanish. The triangular member rigid body rotations are directly related to node displacements so that if the node displacements relative to yfi, y; coordinates will lie respectively on them. The £1, Q2, 4.3 coordinates have covariant the yl,y2,y3 coordinates vanish, the yL, unit base vectors e1, e2, '53 , and contravariant base vectors .. -2 _ . 2 e1, e , e3 . The coord1nates yl, y ,y3 have base vectors .—1 +2 4-3 3(1), 3(1), 3(1) and are the i-member coordinates of Sec. 2. 7. In addition we introduce three oblique coordinate systems associated with the three node points. Their covariant and contra- variant base vectors are as shown in Fig. 3. 2. The inplane covariant components of each system are parallel to lines defined by node points and two of the contravariant components are perpendicular to these lines. The coordinates yl, yz, y3 are redefined after each interval A7 k so that the displacements measured relative to all the coordinates 38 (a) Fig. 3.1. Fig. 3. 2. 39 Member coordinates (b) Covariant and contravariant node base vectors e(2.) 40 defined in this section are small in accordance with linear theory. 3. 2. Geometric Parameters Consider a triangle of side lengths £1, £2, 13; node points 1, 2, 3; included angles 411, LIJZ, 413; and referred to oblique cartesian coordinates £1, £2 as shown in Fig. 3. 3. The coordinates Q1, L2 are components of a position vector 5 = 4131 + 1232 (3.1) where El’ :32 are the covariant base vectors (see Appendix A). The dimensionless parameters 0. and 8 are defined by 1 a 2 4_ I3 (3.2) (1 = %i 2 Wedefine _l 2 2 2 8.1 — :(-£1+£2+13) _l 2 2 2 _1 2 2 2 a3 _ Z(£1+£Z-f3) The area A of the triangle is related by _1 2 2 2 2 2 2 4 4 41/2 A _ Z(2£1£2+2£2£3+2£3£1 -11 -12 -13) (3.4) 41 Fig. 3. 3. Member parameters and oblique coordinates ,\/ \/"l 1' T '2 ”3 Fig. 3. 4. Unit normal and tangent vectors 42 The trigonometric functions are related by 2A 2A 2A sinq.) = , SinkiJ : , SinLIJ = 1 1223 2 £311 3 2122 (3. 5) a a a _ l _ 2 _ 3 COSLpl—TT, COS¢Z—W1 COSLIJ3—fl— 2 3 3 l l 2 The unit normal vectors 1-11, 1—12, 33 and tangent vectors T1, T2, T3 are defined in Fig. 3. 4, and are related to the base vectors by - _ . —1 . —2 nl — SlnLlJze +51n¢3e E = -Sinlll El (3 6) 2 1 ° .. , —2 n3 — - Sln LlJl e _ _ —1 _2 sin 1113 _ sin 1112 _ tl — - cos 1112 e + cos L113 e — - sianl el + sin “1’1 82 _ _ -1 -2 _ .. tZ—-coqu1e -e — -e2 (3.7) .. -1 -2 - t3—e+cos¢1e — el The displacement vector of points on the member in terms of its covariant components is - 1 2 —1 l 2 —2 l 2 —3 u=W1(€.{-.)e +W2(§.§)e +W3(€.§)e (3. 8) —l —2 -3 2 W1 ((1,8) e + w2(o,(3) e + w3(o,(3) e The rotation vector is related to 1—1 by 3w 8w - l 2 1 3 — _ 3 _. 9(§,§)= - ( -—_—e)+e e (309) S1n1p1 81,2 e1 811 2 3 43 Substituting from (3. 2) and (3. 5) into (3. 9) we obtain 3w 3w 3_ e 6(ap)=1(1 3 z - 1 EZ)+0 2K 3 as 1 2 8o (3°10) 3 The rotation component 9 3 is not defined in general. It is, however, given an explicit definition for the three node points at the end of Sec. 3. 5. The normal distances as shown in Fig. 3. 3 are designated by b b b and are related by 1’ 2’ 3 _ 2A _ 2A _ 2A 3. 3. A Symmetric Form of the Strain Energy Expression The plane stress strain energy expression associated with the oblique coordinates C1, £2 (see Appendix A) is W1 _ 6t ff{'—‘12‘—‘(€11 '2C°S Ll’1612‘L'522) _ 2 o . . 1 2(1-1/ )s1n1111 Si 8111 1.111 2 1 Z -#(€11€ZZ-612)}d§ d§ (3112) where 6 is the Young's modulus, V is the Poisson's ratio, , u = 2(1-V), and where the strain components are related to the components of displacement by awl 11 341 6 : __ (3.13) 44 Introducing the change of variables (3. 2) into (3. 13) we obtain Bwl 611 :11 aa 8w2 €22 31—5—55— (3.14) 8w 3w 6 = L( 1 +--——---—-) 12 2 13 88 £2 80. For convenience we define the deformation parameters 61, 62, E3 in the form 2 2 61‘ £3 €11+12 e22'2’2’23'512 e = 1 26 (315) 2 2 22 ' 2 53 ‘ ’3 611 The parameters do not have dimensions of strain, but they are, however, related to the elongations in directions defined by the three sides of the triangle. The parameters are related to diSplace- ments by awl 8w awz awz - ___. ___1 ___. ___. 51‘ ’3‘aa ' 38) +’2‘88 ‘ 6o) 8W2 62 : [2 __‘8—8- (3.16) 8w _ 1 63 ‘ ’3 86 Solving for the strain components in terms of the deformation parameters in (3.15) we obtain 45 1 6 = —-€ 11 2 3 13 e - J—e (317) 22 '— 12 2 ° 2 6 = 1 {E +6 '6} 12 21213 3 2 1 Substituting from (3. 2), (3. 3), (3. 5), and (3.17) into (3.12) we obtain after some rearranging of terms 2 2 a a . 2 1 2 Z W = ___a— ff{_l__(_l.+u)€ +_.( +u)€ 1l 404,2) Si 4A A2 1 4A A7 2 3.2 aa aa 1 3 2 2 12 2 2 3 + To.(—+“)€3+IA'(—T'“)€1€2+ZK(—_2-"”€2€3 A A A Z a3‘11 5t T +-( ~u)€€}dad(3=——-——-ffE G dad8 4A "TA " 3 1 40-1/2) Si 1 1E1 (3.18) where T _ and _. r 32 aa aa- 1 1 l 12 l 13 —4A(—2 ”J 4—A'(—“2"“) inky”) aa a2 33. 1 21 1 2 1 2 3 c = —(—~u) —(—+u) ——+——-u) 1 4A A2 4A A2 4A A2 a'2 ——1 (——a3al -u) —1 (Lie12 -u) —‘ (—3 +11) 4A A2 4A A2 4A A2 (3. 20) 46 The plate bending strain energy expression (see Appendix A) is 2 2 { 1 8 W3 8 W3 W = . f f - 2cos L): -—— I2 2 811le1 Si 31112111 81,1811 latlaéz 82w 2 82w 82w 82w 2 3 3 3 3 1 2 + “'2'"'2 ' 1 1 2 2 '(' ‘1 2) Nd” d4 8!, 8L 81’; 8Q 8§ 8Q. 8Q 8Q (3.21) where ,d is the flexural rigidty. Note that 2 2 8 W3 : 8 W3 841811 1; a 2 82w 82w 7‘35 = T‘T3‘ (3.22) 81; 8C, 12 88 2 2 8 w3 _ 8 W3 81.1812 — 1 1 80. as 2 3 We define the deformation parameters €4,65,€6 by 82w3 82w3 2 82'w3 e z + - 4 8&2 a‘32 '8'688 82w3 65 = —2— (3.23) 315 E -_- 82w3 6 8112 From (3. 22) and (3. 23) we obtain 82w ___.}. z -1_ E 841811 1: 6 82w3 1 : -- € (3.24) 31.2312 1: 5 82w3 1 —— = (-€ +6 +6) 8§18§2 21213 4 5 6 Substituting from (3. 2), (3. 3), (3. 2 a2 2 a a _ D 1 1 2 1 a2 2 1 3 2 W12 ‘ 4ff{4A(—2+”’€4+ “41""? “’65 4A’AZ+”’€6 si A. A a a 32 a3 3a 2 2 1 44,}1 2-u)e4€5+ -4-A(—23-- 10656 6 +4-A(-a—3-z—— -IJ)€6€4 }d<1d A2 A2 =9ffET GEdad8 (325) 4 1 2 . where G1 is given by (3.20) and ET 2 [6 e E ] (3 26) 2 4 5 6 ° If the middle surface deformation parameters (2. 10) are taken in the form (3.16) and (3. 23) then the matrix operator 9 in Chapter II is given by 48 8 8 8 8 - ‘3‘3‘3‘+5'§)61 ”ass-saw a .. ’2 87362 a _ 1 —- e 52 = 3 8a 1 (3.27) 32 -2 32 + 82 E 80.2 80.88 882 3 2 m e 882 3 a2 E 80.2 3 Also the matrix G in (2.11) then has the form ___.12. 91 o 4(1-1/ ) G: 60 (3.28) 0 4 G1 where G1 is given by (3. 20). 3. 4. The Member Displacement Configurations 2 The vector of displacements measured relative to the L1, Q , {.3 coordinates in terms of its covariant components is EL(§’. 4") = “(11, 12) 31 + w2(§’. 1.2) 3" + w3(§’. 1.") 3.3 (3.29) In terms of the dimensionless parameters (3. 2), (3. 29) is EL(6,8) = w1(c1,8) :1 + w2(o.,8) 32 + w3((1,8)33 (3.31) 49 Six of the deformed configurations Hi are associated with membrane stresses and consequently have the component W3 = 0. The other six are associated with bending stresses and have the components wl = w‘2 = 0. Three of the deformed configurations used have displacements desc ribable in" the' form .. ..1 -2 uL = (c0 + c111 + c28) e + ((10 + dlo + d28) e (3.32) where ci and di are arbitrary constants. This displacement vector is linear in o. and 8, and consequently has the property of displacing straight lines into Straight lines. It is the same one used by Turner et al.13 The six arbitrary constants in (3. 32) are determined in terms of the three parameters >11, )1 2, )1 3 called member generalized coordinates. These parameters are respectively defined as the elongations of sides 11’ £2, 13 and are the natural forms given by Argyris.l4 The node diSplacements and constraints, relative to the £1, L2, 43 coordinates are shown in Fig. 3. 5, and given by _. _ —1 —2 _ —1 -2 uL(0,0)—w1(1)e +w2(1)e — coe +doe 71(10):“) ‘8‘ = (c +c)'é'1+(d +11)?2 L ’ 1(3) 0 1 0 1 11(01)=0 =(c+c)E’+(d+d)E2 L ’ 0 2 0 2 (3.33) 50 42 3 .5" A -1 w1(3) e 1113 w -2 2(1)8 4» ’2 l 4 1 —l 2 W1(1)e Fig. 3. 5. Covariant components of node displacements relative to oblique coordinates. ___....g 51 Solving for the six constants in (3. 33) and substituting into (3. 32) we obtain .. .—1 —2 = - - - - 3. 34 u [+fiwl(3)+(l c1 8)w1(l)]e +[(1 a 8)w2(1)]e ( ) The member generalized displacements are then related to the components of node diSplacement by 13 A1: 'tl. W(O,1): ‘1’]: W1(3) 12: -t2° w(0,0) = -1112“) (3.35) )‘3: 't3. W(090): -W1(l) where f , F2, t3 are vectors tangent to the three sides of the triangle and are given in (3. 7). Substituting from (3. 35) into (3. 34) we obtain uL = 71; {81111+(1-a-8)13x3}21 1 _2 + :{U-a-mlzkz} e (3.36) In general when a member is displaced from its initial position (Fig. 3.1a) to its displaced position (Fig. 3.1b) it undergoes a rigid body rotation AER defined below. The coordinate system le, yRZ, yR3, see Fig. 3.1b, is defined so that the rigid body member rotation and node 1 displacements relative to it vanish. Nine deformed configurations are directly related to components of node rotation measured relative to the le, yRZ, YR3 coordinate system. Three deformed configurations are of a type graphically shown in Fig. 3. 6 a, b, c and only have displacement components in 52 1 .. _ ___.). 1 - 51an 881(1) 2 —. 1 e l /y Sintlll 9 2(1) V2 __1_._ . 1 E \ \ \ 1 \ 1 1 >./"' \Y y -—.-— e yl SlnLIJZ 9 1(2) ((1) (e) (f) 3 3 3 Y tY t Y 2 \_Ti:|1_xllgl(l) ’ >1 '5 S1n - ___..Y _:_l____.. 2 ‘flanl iii-(3v)- y2 C 1 — \ T - --.——-—>» e I sin 4J1 "10e1(3) S1an3 11 2(3) / 1 ... ' -¢———>. e \31n 412 10 2(2) \ 1 / k y1 1 )1 e \Y1 Y sin ((12 12 1(2) (g) (h) (i) Fig. 3. 6. Displacement modes associated with node rotations. 53 the plane of Si (see page26). The remaining six only have displace- ments components normal to Si and are shown in Fig. 3. 6 d-i. The nine functions of Hi shown in Fig. 3. 6 are constructed from the equations of the six phnes intersecting the plane of S; as shown in Fig. 3. 7 and given in Table 3. 1. The equation of the planes are normalized so that their principal intersection angle is unity. The set of points defining the domain of a subtriangle as shown in Fig. 3. 8 are designated by Aal’ AaZ’ Aa3’ Abl’ AbZ’ Ab3; i. e. , Aal corresponds to the set of points lying on triangle 1, 2, 4, etc. The symbol Aa when used witha function gal in l the form Aal g implies that ga‘1 is defined over Aal' The form a1 (Aal gal) (AaZ gaZ) = Aal AaZ gal gaZ implles that ga1 ga2 is defined only for points common to both Aal and Aa2° Consider the two functions taken in the form ._ .1. Z12 ‘ 2 {ga2(gal + gb3)Aa3 + gal(ga2 + gb3)A'b3} 13 (3'37) _ 1. 213 _ 7- {gal (ga3 + gbzmaiz + ga3(ga1 + gbZ)AbZ} 12 (3°38) The functions Z12 and Zl3 together with their first partial derivatives are continuous for all points of the triangle. The second derivatlves of Z12 are constant on Aa3 and Ab3’ and have a discontinuity along line 3, 6. The second derivatives of Zl3 are constant on AaZ and AbZ , and have a discontinuity along line 2, 5. It follows that all third derivative of 212 and Zl3 vanish for all points of the triangle and consequently strain compatibility is satisfied for inplane displacement configurations constructed from Fig. 3. 7. Lines of intersection of planes given in Table 3. 1. Fig. 3. 8. Symbols used to designate subdomains of triangle. Table 3.1. Functions Used to Construct Displacement hAodes. LINE FUNCTION EQUATION ‘ 2, 3 gal = 1-0. 43 = O 3, 1 gaZ = a = O 1, 2 ga3 = s = o 1, 4 gbl = a - fl = O 2, 5 gbZ = -1+o.+2[3 = O 3: 6 gb3 = 1-2a -p = o 55 these functions. A graphic description of these functions is shown in Fig. 3. 9. The function Z12 vanishes along lines 1, 3 and 2, 3 as indicated in Fig. 3. 9 and its normal derivativesvanishes along 2, 3. It has a sIOpe of unity in the direction of g1 at point 1. The configuration associated with X 4 is graphically described in Fig. 3. 6a and has the form uL = >\4(Z13 nZ - Z12 n3) (3.39) where 32 and 33 are defined in (3. 6) and Fig. 3. 4. By two consecutive permutations of the indices 1, 2, 3 into 3, 1, 2 and 2, 3, 1 in (3. 36) and (3. 37) we construct four additional functions Z21, Z23, Z31, Z32. On substituting from Table 3.1 we obtain 2 = l {a(2-3a-2p)A + (1-01 433 } 1 12 2 a3 A'b3 3 _ L 213 _ 2 { -(1-o. -13)(1-a -3B)Aaz + [3(2-211 -3s)Ab2} 12 z =1—{s(2a-s)A +02A }1 (3 40) 23 2 a1 bl 1 ° Z21 : 12(32 Aa3 - (l-a -B)(1-3a -[3)Ab3} 13 z31 - 12-{41-(1 -s)(1-a -3[3)Ab2 + 52 Abz} 12 é-{fiZA -a(a-213)Ab1} I, Z al 32 The displacement vector associated with the three additional inplane deformed configurations then has the form 6L ‘ "4(213‘12212113)+ "5‘221’13 ‘ 223111) +k(3(232’11 ' 23an) (3.41) 56 2 \LI Fig. 3. 9. Graphic description of functions Z12 and Z13. (a) (b) Fig. 3.10. Force vector subjected to a. virtual diSplacement. 57 On substituting in (3. 41) for the functions and unit vectors, rearranging terms and combining with (3. 36) we obtain EL = II; { -1311x1 - (1-01 -p)£3x3 -[ (1-11 -mzAaz + (3(2-2a 739M132] M4 2 - [ [3(2c1 -[3)Aa1 + o. Abl] Axs + [ fiZAal + o.(213-c1)Abl - (1-11 pm -11 -3B)Aa2 + sZAbZ] Axé} “e1 + 721: { ..(1 -a -fi)12>~2 + [ o.(2-3o. .2(3)Aa3 + (1-01 -B)2Ab3] A114 -[(12Aa3 + (1-11 —s)(1-3a+s)Ab3 +p(2a -B)Aa1 + azAb1]A)\ 5 + [(321181l + c1(2[3-o.)Ab1] Axé} 32 (3.42) Eq. (3. 42) contains six of the independent vector functions of Hi' The remaining six configurations are taken in the form E =w E3 = { g g X + g )x + g )1 L 3 a2 a3 7 ga3 al 8 gal a2 9 + [ g1:1(gaL3 Aal + gaZ A131” )‘10 + [ gb2(gal AaZ + ga3 A1.-;2)])‘11 -3 +[g (g A H; A )]X e b3 a2 a3 a1 b3 12 (3.43) The 81x functions gaZga3’ . . , [gb3( gaZ Aa3 + g3.l Ab3)] and their first partial derivatives are continuous for all points of the member. The last three of these functions respectively have second derivative discontinuities along lines 1, 5; 2, 6; 3, 4. Substituting for the function in (3. 43) we obtain 58 +[ -(1-0 -2£3)(l -0 -F5)Aa2 - (l -a -ZB)BA~bzl>\11 + [ (1(1 -Zo. -[3)Aa3 + (1 -o —s)(1—2(3-(3)Ab3] x12 33 An alternate form of (3. 44) is given in Appendix D. 3. 5. Member-Shell Transformation Matrix (3. 44) uL = {apx7 + p(1—a -(3)x8 + o(l -o. -mxg + [ (3(a -[3)Aa1 + (1((1 -(3)Ab1] km The member generalized forces 2i and member generalized coordinates Ai are related by the strain energy relation T (SW = FT 6U I. 1 1 1 andif 6A = M. 6U. 1 1 then F 6U = 2:511. 1 1 1 2 2T M. 6U. 1 1 1 (3. (3. (3. (3. If (3.48) must hold for all virtual displacements 6Ui it follows that (3 45) 46) 47) 48) . 49) 59 On taking the transpose of both sides of (3. 49) we obtain F. = MT 2. (3.50) 1 1 1 Here the components of Fi are components of node moment and force vectors, and the components of Ui are components of displacement and rotation vectors. Consider the displacement vector 61} and force vector f shown in Fig. 3. 10b. The work due to f when it undergoes a virtual displacement 65 is defined by the dot product 5w = “f- 613 (3.51) The covariant base vectors (unit vectors) of an oblique Cartesian coordinate system are —el, 32 and the contravariant base vectors (not necessarily unit vectors) are 31, 32. They are related by —l — _ —Z - _ 1 e e1 — e . e2 — —1 _ _ _ . _2 _ o (3.52) e e2 - el e .— The vectors f- and 61-1 in terms of their components are .. 1 _ _ 6n = 6w e1 + 6w e2 (3. 53) = 5w1 El + 5w2 EZ _ 1 _ z _ f ..f e1+f e2 (3.54) : f 1 + f E2 V?“ 60 Substituting from (3. 53) and (3. 54) into (3. 51) we obtain ..1 —2 1— Z— 1 2 6W = (fle + fze )' (6w el + 6w e2) = flow + fzow — — — - 2 =(f1el + f2e1)' (owle1 + 6w2e2) = f1 owl + f 5w2 = (£131 + fZ-e'l) ' (owlel + 6w2E2) : . 12 Sln L); [11 6w1 + fzéwz + c051); (f1 6wz + f26w1)] = (fl; + £2; ) ' ((5le + 6wZ-e ) : f16w1 + f26w2 l 2 1 2 + c0511: (f16w2 + fzéwl) (3.55) where LIJ is the included angle of the oblique coordinates. We note that 6W:[f1 12] 6w1 = [f1 12] 6w1 (3.56) 2 L6w 6w2 but in general 6W )4 [f1 12] 5w1 (3.57) 6W2 6W )4 [f1 12] 6wl (3.58) 6w2 We then conclude that if the generalized coordinates are contravariant components of a displacement vector then the associated generalized forces are covariant components of a force vector and vice versa. 61 From (3. 35) and Fig. 3. 5 we obtain the node displacements relative to the £1, £2, L3 coordinates in the form .. ..1 -2 uL(l) = X3 e(.1‘)+ X2 e“) 5M2) = o (3.59) -2 u13(3) : >‘1‘"3(3) where the base vectors are those of Fig. 3. 2a. It then follows from the above discussion that if the components of node displacements are interpreted as components of generalized coordinates then the associated components of generalized force can be interpreted as components of node forces through the energy relation aw =? -5Ii I (1) (1) +?(2)- 6'13 611 (3.60) + 'f - (2) ‘ (3) (3) To do this we take the variation of (3. 59) and then substitute it into (3. 60) and obtain _ 1 - 2 —. . ...1 -2 6WI - (f(l,e1)el(l)+f(l,el)eZ(l;)) (6X3 8(1) + 6X2e(1)) 1 - Z .. ...2 + ”(3. e3) 61(3)+ f(3, e3) e2(3)" ”"1 6(3)) ‘3' 6” = f1 5X + f2 6R + f2 5X (l,el) 3 (l,el) 2 (3,e3) 1 l 2 where f“: el)’ f“, el)’ etc. are force components referred to the base vectors of Fig. 3.2. If 6114 = 6X5 =... =6Xlz = O (3.62) then u M ”H 6WI . 6A1 0'l 6X1+ 0'2 6112 + U36K3 (3.63) 62 Eq. '8 (3. 61) and (3. 63) can hold for all virtual displacements that also satisfy (3. 62) only if f1(1,e1) 3 f2(l,el) 0'2 (3.64) f2(3,e3) : ‘71 The components of the node forces ?(1)’f(2)'?(3) expressed in terms of 0'1, 0 2, 0'3 can be obtained directly from (3. 64) and the equations of equilibrium. They are shown in Fig. 3.11a and have the form 1'(1) : ”331(1)+"232(1) f(2) : “131(2)+"2ez(2) (3.65) f = (3) “281(3)+"182(2) The member generalized displacements A4, K5, . . . , X12 are related to the components of the node rotation vector measured relative to the yi, yli, y; coordinates as shown in Fig. 3. 6. The node rotation vectors are _ _ 1 _ _ _ ¢(1) ‘ sianl 0‘11 ‘ x8) e1(1) sin—__ITJI “12+ kc9)‘"2(1) + >‘4‘°'3(1) - _ 1 _ 1 _ 15(2) ‘ amp—'7 0‘12 >‘9) e1(2) + simpzo‘m + )‘7’ 82(2) + >‘5‘3 3(2) _ _ _L_ _ _ _ ¢(3) sin¢_—: 0‘10 ‘ >‘7)‘31(3)+1311r1‘_"1(3"3()‘11 + x8) e2(3) + x6630) (3.66) 63 \ /Zel(3) 3 “3 e1(1) 1 2 .63 32(2) 0" — _- 2 82(1) \o-l elm (a) q132(3) q E \ / 2 1(3) 3 q3 e1(1) 1 2 , 'q3 E2(2) 432(1) C11 31(2)\ (b) Slnl.IJ3(O'll + 0-8)E(23) sin¢3(010 - 07);:3) —3 6 8(3) sin1p1(0' l -0’8)El __2 \( Sinlllzfirlo + 0’7) e(2) _3 ‘ —3 “4 911111, \1) “59(2) sinu|41(0'12 + 09%?” sinxpzwlz - 09);:2) (C) Fig. 3.11. Relationship of generalized forces to node variables. 64 If the components of the node moment vectors 5(1), “1(2), 5(3) are interpreted as generalized forces and related to the member by means of an energy relation similar to (3. 48) then by using arguments similar to the above we can obtain (1) = sin1p1(crll - 08);:1) + sian1(0'12 + 09)E(21)+ 0'43?” Bl “1(2) 2 sinipzwlz - 69)-5(2) + sinupzwlo + 67);.(22) + 0' E?“ 5(3) = Sin¢3W1o ' ofire-(3) + Simp3(‘711 + “8):;(23) + 663(3) (3.67) The above are shown in Fig. 3.11b. The node forces, necessary for equilibrium when the member is subjected to node moments given by (3. 67), can be obtained from the six scalar equations of equilibrium. The force components normal to the plane of the triangle are associated only with the in-plane moment components and can be obtained directly, by taking moments about the three sides of the triangle. The in-plane force components are associated only with the moment components normal to the plane of the triangle. Since there are six in-plane force components and only three independent equilibrium equations associated with them, it follows that this system of forces has three redundancies. The redundencies are essentially removed by imposing certain symmetry requirements on the components, i. e. , the six components are interpreted in terms of three parameters ql’ q2, q3 related as shown in Fig. 3.11b. 65 The six equilibrium equations are, after some simplification, givenby 3 _ .1. f(1,.21) ‘ 2(1";"12 '12 0'11) 3 1 1 f(2, e2) 2‘1“; “1o '1; “12) e _ _l_ l f(3,e3) ‘ 2‘22 “11 ”Tl—“10) (3'68) 11 q1 ___ ZX(°4+“5+“6) 12 ‘12 = Z‘K("4+“5+“6) 13 q} : ZA(U4+“5+“6) From (3. 65), (3. 67), (3. 68) we obtain the relationship between node forces and moments referred to the coordinates of Fig. 3. 2 a, b, and the member generalized forces. This relationship is given'in. matrix form by (3. 69). 1 f(1,121) f(1,121) f(1, e1) m1(1, el) m2(1, e1) m3(1, el) f(2,132) f(2,152) f(2. e2) m1(2, e2) m2(2, e2) m3(2, e2) f(3, e3) f(3, e3) f(3, e3) m1(3, e3) m2(3, e3) m3(3,e3) h N N NJ: :> NA >| N4: (11> (3.69) 10 ll or12 67 The coefficient matrix in (3. 69) is designated by ME . The relation- 1 ship between the components of node forces and moments referred to the i-member coordinates (yl, y2, y3) and those in (3. 69) can be obtained directly from (A. 22), (A. 23), and from the definitions of the quantities. The relationship in matrix form is given by (3. 70). The coefficient matrix in (3. 70) is designed by Mr; . 1 The matrix Mi is related by MT = MT MT (3.71) and is given explicitly by (3. 72). 68 Aos.mv “Wo.mvm:s Amo.mv~zs Amo.mvzcs Amm.mw Amo.mr Amm.mr Amm.mvmss Amo.mv~qs Amo.mvzas Amo.~:. Amm.~r ANm.~yw Azo._vmcs Azm.zv~cs Azm.zv_:s Azo.zyw Azo.zyu N m a As aw fl. d mafia mafia uamoo Namoo madam m9GMm Hanan macaw Namoou ~%moo 37:3 madam msmoo a a o fiaswm ~9moo gsswmu o o dsmog. f(11) f(11) f(11) m(11) 2 m(11) m(11) f(21) f(21) f(21) m(21) m(21) m(21) (31) f(31) f(31) m(31) m(31) m(.11) N W 3211 a2 a2 a2 '1’4A1 “4.4.1 ”4'A1 O 0 0 O 3 3 3 0 1 1 1 27— T 21— 3 3 3 o o o o A o 17— 0 0 23 a 1 1 o - _ o 323 2 o 1 o o 1_a1_a1_a1 411.1 2 4.4.1 3 4A1 3 o--—-1---1—-—1— 2123 23 23 o o 0 fl— 1 fi- 0 01% 13 13 8‘2 0 _1_ 6‘2 21113 221112 o o 1 o 1 1 0 f3 _3_ .3. 4A 4A 4A 0 o o o _1 0 O 0 I— 1 73-1—1 01—3- 13 23 13 _32 31 032 21113 2213 2113 o o o 1 (3. 72) 10 ll “12 70 Two of the triangle rigid body rotation components have a direct interpretation, i. e. , they are the two inplane components of rotation of the plane defined by the three node points of the triangle. The component of rigid body rotation (A9111) normal to the plane does not have a direct interpretation and is in fact intrinsically related to the manner in which the redundancies, mentioned above, are removed. In order to obtain this interpretation we let a triangle in equilibrium undergo a virtual rigid body rotation 69; so that the magnitude of the forces and internal strain energy do not change. Only inplane components of node forces and displacements need be considered. The virtual displacements of the node points with respect to the oblique coordinates of the node points (Fig. 3. 2) can be represented in the form .. ..1 —2 5‘10) z 6""1(1)‘*(1) + 5W2<1)e<1) ._ _ _1 -2 5u(2) — 6w1(2) 6(2) + 6w2(2) e(2) (3.73) _ _ _1 -2 611(.1) ‘ 5W1(3)e(3) + 5W2(3)e(3) The node forces can be represented in terms of the six parameters 0'1, 0'2, 0'3, ql, q2, q3 (Fig. 3.11a, b). They are f(1) = (”3 + ‘13) 31(1) + (“2 ‘ q2’2520) 1(2) =(0'1+ (1151(2) + (0'3 - q3)'éz(2) (3.74) 5(3) = (“2 + q2)-‘;1(3) +("1 ' q1)-‘52(3) 71 The components of the node moment vector (0'4, 0'5, 0’6) normal to the plane are also required. The work due to the node forces and moments when the triangle is subjected to a rigid body rotation 69; must vanish. This can be expressed in the form 3 - - - - - - _ (0'. +0'S+0'6)69R +f 6u +f ~6u +f -6u(3)—0 4 (1) (1) (2) (2) (3) (3.75) Since 69;. does not cause the sides to elongate the work due to the equal in magnitude and Oppositely directed force components associated with each of the parameters 0'1, 0'2, 013 must vanish independently. On removing these parameters from (3. 74) and substituting the resulting expression and (3. 73) into (3. 75) we obtain after some simplification 3 (0'4 + 0'5 + 0'6) 69R - q1 (5w2(3) - 6w1(2)) - q2 (5w2(1) - 6w1(3)) - q3 (6w2(3) - 6Wl(l)) = 0 (3-76) Substituting the last three of (3. 68) into (3. 76) and dividing through by (0'4 + 0'5 + 0'6) we obtain 3 I 12 l 59R = 32‘: (5W2(3) - 5W1(2)) + EX (5W 2(1) ' 5W1(3)) I 3 + a (5W2(3) - 5W1(1)) (3.77) From the assumptions of linear theory we have 72 3 £1 £2 AGR = 4A (Aw2(3) - AW1(Z)) + 4A (AWZU) - Awl(3)) £3 + :3. (“’20) -Aw1(l)) (3.78) The rotation of the coordinate system y;, ya, ya relative the coordinate system yl, yz, y3 referred to in Sec. 3.1 is AGR defined above. 3. 6. Member Stiffness Matrix The triangular member stiffness matrix defined by (2.75) is a 12x12 matrix I"i with det Pi > 0. We partition (2. 70) in the form 6 =[115 Hg] .Al. 9.79) i 1 A % wh IETA d tth a 1 t t ° b ere 1i 11 corre8pon s o e 18p acemen vec or given y (3. 42) and HT_ A corresponds to the displacement vector given by 1 2 q (3. 44). If the deformation parameters are taken in the form T_ T T E .[E1 E2] (38m where E1 and E2 are respectively given by (3.16) and (3. 26) then from (3. 27), (2. 70), and (2. 71) it follows that 1 U (3.81) From (3. 27), (3. 42), and (3. 44) it follows that (3. 81) will reduce to FIE; 321in1 o T .1111— b113,] :1. 0 52711112; _A21_ (11111 0- ”All — .0 H223. _AZL (3.82) On substituting the coefficient matrix of (3. 82) into (2. 75) we obtain 1' 3 r1 - r - r 1 1"mi I121 11.111 0 G1 0 H11. 0 : £1. T i, ‘ dA 1“21. I‘ 1 0 I122. O C'1 O 11’22. 1 22i 1 U "' T éfIIll.Gi'fl'lidA o 1 1 11i _J (3.83) The matrices I”11 and F22 are respectively called the member i . o 1£11111. Gi'II' dA membrane stiffness matrix and member bending stiffness matrix. These matrices are given in Tables (3. 2) and (3. 3), and a detailed derivation is given in Appendix C. 74 75:. flsmeN. mmNmmummHmvu N a ~+N Non mmNmN+NmHmN+ :fm: vaWlH—Nlfi m. N IIIII fl 1| H1N+AHM1~3W|1NIH Nd NV .4. Eufldnmmvuwl Wm. gm an H Nu m M .II N+N 2 < “3.3.32 mquwmfi—m vnmhngoz .N ..m oHQmH ._I ||||| 1|. _ H1©~+A~wm .mvu " :mowhuognbnmv _ _ NsNNN-NN.NN- _ _ _ . mmm+NmN+~mmvE_ «FWNIW.<_IIIIII _ Hamms +1fi. 1 dams - _ _ N _ «2.. N _ _ _ NstN+NNHNN-_ NNva-NNzNN- H . . _ N N 2 NF N N 2 Ne _ N I m .m N .m d l _ N N N +N CAN “N N+N N+N NHa. _ lllllllllllllllllll L 1| 11 I | I _T N _ ... _ ... l. _ EN -1 ...- 3w.) (MNTI £23-.”leme _ 2+ Womb _ _ u; .m N_ Nd NN _ 3 11111 .. ...1........|..1.111!..111.-..-...-1 _ 2..-...1N1DiN 1N..- s.Nlm1_:1- NNN1N14. - %)[ c1(2-3c1)Aa3 +(1-a) (l-30.)Ab3] )5} 3"(21) + {o(l-o))\9 + ozA — (l-o)2Aa ). blxlO 2 ll 3 + [ (1(1 -2o)Aa3 + (1 -a)(l-2c1)Ab3] 112} '1'“) (3. 85) 77 Let L1, L2 correspond respectively to the line segments 1-6, 6-2 shown in Fig. 3. 7 so that f(a))L‘.I,. g(a)(Ll+L2) implies that f(a) is defined over L1 and g(o.) is defined over L1 and L2. Along the K14 axis the following then holds A'al = L1 + L2 AbZ = 1..1 + L2 (3.86) Aa3 : I"l Ab3 = L2 From the above it then follows that (3.85) can be expressed in the form E = {%[a(2-3Q)Ll + (14021.21).4 - [11(2.3a)1.1 NIH N2 + (l-o.)(l-3c1)LZ]).5} 1(1) + {o(1-11)(L1 + L2».9 + [11(1—211)L1 +(1-a)(1-20)L2I"1.{2(3-(31) (3. 87) The components of node rotation about the y]: axis (192“ i)’ $2” 1) 1 . and about the yR ax1s #33”, .i.)’ (1)3(2’ i) are related to )1 4, X 5, )1 9 and X by 12 )- -1 F — — '- ¢2(1,1) ° 0 '1 '1 "4 4542’” 0 o 1 -1 )5 = (3.88) 413“,” 1 o 0 o )9 ¢3(2.i) 0 1 O 0 >‘12 ... .1 _ ...l _. 78 The above can be obtained directly from Fig. 3. 6. These rotation components are identical for adjacent members only if the components of rigid body rotation (see Sec. 3. 5) about axes normal to the plane of the trianges are identical for the two members. On solving (3. 88) for )(4, X5, X9, X12 and substituting into (3. 87) we obtain after some simplification Ii = -l-{[o.(2-3a)L +(l-a)ZL ]¢ 2 1 2 3(1.i) N2 -[a(2-3a)Ll + (l-a)(1-3G)LZ]¢3(Z,1)} (1) {-[ a(2-3a)L1 + (l-G)ZL2] 412“, 1) NIH -.3 2 L1 + (l-o.)(l-3o)LZ] ¢Z(Z,i) } J (3.89) - [ -a A graphic description of the manner in which the displacements along the line are related to components of node rotation, is given in Fig. 3.12. From (3. 89) it can be seen that along the line the inplane di Splacements are related to the components of node rotation in a ITrianner identical to displacements normal to the plane. Since the COmponents of rotation for common node points of adjacent members ar e required to be identical it follows that displacement continuity is ITlaintained along the entire line even if the initial intersection angle d~Oes not vanish, provided that the normal components of rigid body 1‘ Otations of the two members are identical. If in addition the normal SIOpes of two adjacent members vary C ontinuously and linearly along their intersection line then normal slope continuity is maintained for all points of this line. In the 79 k v (a) (b) (d) Fig. 3. 12. Relationship between displacement components normal to an edge and node rotation components. 80 formulation used, normal slqae continuity is completely satisfied only under a restrictive set of conditions. These conditions are related to the geometric shape of adjacent members and the type of deformation present along the intersection line. The inplane diSplacements (3. 42) do not result in slope change along the edge. The displacements normal to the plane of the triangle are constructed by polynomials no higher than the second degree. Consequently the normal slopes associated with the deformed configurations of 17, . . . , )1 12 (see Fig. 3. 6) vary linearly along line 1-2, but the variation of the configuration associated with K can have a finite discontinuity 12 at the midpoint. The manner in which the normal slope is related to K12 is shown in Fig. 3.13. Also the relationship of the midpoint discontinuity to the geometric parameters is described there. The magnitude of the discontinuity is designated by d and is related to the components of node and midpoint rotations by 1 d : ¢1(6,i) ' 2(¢1(1,1)+¢1(2,1)) (3'90) From Fig. 3. 13(a) we obtain the relationships between the components in (3.90) and X12. They are .311 <1’1(1,1) z 4:01.111le = '24 )‘12 8‘2 ¢1(2,1) : COtLp2X12 : 2A "12 (3°91) (Ii-1f) ¢1(6,1) : 1:01.116 "12 z " 4A )‘12 81 1 (3.1 x12="’2(2.1) 4’1(1.i)|_4 (1)1159” 1J¢1(2. i) 1* I rlfi 1 f1 3 3 y3 2" (a) 2— ' R 1 >// 12 .1i -'“” 1 V 7‘ YR (b) 1(2,1 ¢ ' M 3 IIIIIl ""' V’xi, I In“ 1’ + ,”’ d (ll ””“'“ m - <1°1(1..1) d’1(2,1) 1 | (C) (d) Fig. 3.13. Variation of normal slope for secondary bending modes. 82 Substituting from (3.91. ) into (3. 90) and simplifying we obtain 1:4? ). T 12 (3.92) d: In order to interpret the conditions necessary for slope continuity to be satisfied we first consider two adjacent members with their middle surfaces on (near to, in the linear sense) a plane as shown in Fig. 3.14. The various parameters of one member are designated in accordance with earlier notation and its adjacent member has its corresponding parameters distinguished by a prime. By examining Fig. '3 3.13 and 3.14 it can be seen that 112 = ")12 (3.93) Then for normal slope continuity to be preserved d Z -d' (30 94) Substituting from (3. 92) into (3. 94) we obtain 2 2 2 1 2 I -1 (1') - (I ) 2 l A = _ 2- 1 k1 (3.95) 2A 12 2A' 12 Substituting from (3. 93) into (3. 95) and dividing through by X12 we obtain 2 2 1 2 I Z 3 . 6 12 -11 _ ”2) -(11) . ( 9) 2A — ZA' If the adjacent member does not lie in the same plane/then (3. 93) is not in general satisfied and normal slope continuity is insured only if d : -d' : O (3.97) 83 3’3 Fig. 3.14. Notation for two adjacent triangles. 84 This can be satisifed in general only if 2 2 12-11 2A ' 2 2 3. 8 (15) 41;) ( 9’ 2A or simply if 1 2 (3. 99) II b -I II If the intersection angle is not very large then 112 z 4312 and (3. 96) will be approximately applicable. If it is large, as for example at a 900 corner, then normal slope continuity can be insured for all points of the intersection line only if (3.99) is satisfied. As already implied, the normal slope discontinuities along side 11, 12, 13 are respectively related only to the member generalized coordinates; X10, X11, X12. These can in turn be interpreted as being related to curvature changes along their respective lines. Consequently large discontinuities can occur along a line only if large curvature gradients are present there and then only if geometric pr0perties discussed above admit to it. For a given triangulation of a shell the slope discontinuities and curvature gradients associated with a solution can be obtained directly. From this information a better triangulation can be obtained by satisfying or more closely satisfying conditions (3. 96) and/or (3. 99) between members with large curvature gradients. 85 We classify the bending displacement configurations into those associated with X 7, >18, X9 which we call the primary bending modes and those associated with 110, X11, 1.12 which we call the secondary bending modes. As already indicated, adjacent member slope discontinuities are associated only with secondary bending modes. We now examine the normal slope continuity between adjacent members with respect to a limiting process that causes the area of the triangle to vanish without materially disturbing the ratio of its side lengths. We assume that the deformed intersection line forms part of a regular curve and in the above indicated limit approaches a straight line. There is no ad hoc reason why this should be the case and to prove such an assertion would undoubtedly require several sufficiency conditions on the limiting process together with a proof demonstrating that it is an inherent requirement for minimizing the potential energy. For this the explicit representation of the bending stiffness matrix should prove of value. In any event this assumption leads to a useful speculation. Consider the displacements along the intersection line of two adjacent members relative to the line joining the common node points. Since the components of displacements tangent to this line do not affect normal slope continuity we need not consider them. The displacements normal to this line can be decomposed into two orthogonal components which, as has already been indicated are identically related to associated node rotation components so that we need only consider one of these components. 86 Consider a curve described by the function f(x) in the interval (x, x + Ax) as shown in Fig. 3.15a. We approximate the curve by summing the two functions associated with the parameters 61, 62 and graphically shown in Fig. 3.15b, c. The approximating function and its first derivative are required to be identical at the end points of the interval. The functions assoc1ated with 61 and 62 respectively vary relative to line AB (Fig. 3.15a) in a manner identical to the normal displacements along an edge of the primary and secondary bending modes. For the alternate bending modes given in Appendix D this is equivalent to a Hermite inter- polation. From the Taylor series expansion we obtain f(x + Ax) = f(x) + f1(x) Ax + .311.17.- + 2 (3.100) f'(x +Ax) = f"(x) + f" (x)Ax + f"'(x) 923‘,— + . .. where primes denote derivatives. From Fig. 3.15 we obtain f'(x) = 0-6 +6 l 2 (3.101) f'(x+Ax) = a +61+62 f(x + Ax) - f(xL Ax Solving for 61 and 62 we obtain 51 - -:—[f'(x + Ax) - f'(x)] (3.102) 1 f +A -f 62 = 31f'(x+Ax) +f'(x)] - (x A? (x) 87 + l : 1m+am . 2 I (b) (C) Fig. 3.15. Type of approximation implicit in bending nodes. 88 Substituting (3.100) into (3.102) and simplifying, we obtain 81 = f"(x)Ax 0%) + f"'(x)Ax2 (é- . %) + f”"‘(x)Ax3 (115° £7) + . . . 62 = 1"'(x)Ax2(271~2-,- - 171-H f""(x)Ax3 (2' 131 - 343-) + (3.103) As Ax approaches zero A 51 z f"(x) 7’5- 2 (3.104) 6 z Ax 2 fIII(x) _1._2___ Therefore, when Ax is small, 62 is small compared with 61 . In its relationship to the limiting process of the triangular element, Ax corresponds to one of the side lengths of the triangle. From (3. 104) it would appear that in the limit the primary bending modes will dominate the behavior, and that the secondary modes and consequently also the slope discontinuities become higher order effects. We note that for the alternate bending modes given in Appendix D the normal slope discontinuities are associated with cubic terms whereas the highest terms in the primary bending nodes are quadratic . IV. FORMULATION FOR SHELLS OF REVOLUTION 4.1. Interpretation of Problem in Terms of Discrete Elements The method presented in Chapters II and III is used to formulate the problem of large deflections of shells of revolution having a small imperfection. The undeformed shell geometry is thus describable by the meridian curve of the middle surface (Fig. 4.1), a function describing the thickness, here assumed to vary only along the meridian, and a function describing the imperfection described later in the Chapter. We limit the solution to one describable by a half period strip. The strip is in turn described by a simplified model consisting of 48 triangular members (see Fig. 4. 2). The pattern used consists of 12 rows of four triangles. In order to accommodate a given geometry the triangles vary in size and shape from one row to the next, but the undeformed triangles of a given row are identical. In the computer program used the location of node points along the meridian can be arbitrarily Specified. If one end of the shell is closed the triangles in that end row degenerate into lines and are consequently removed. In the procedure used it is necessary to invert a matrix of rank equal to the number of unknown generalized coordinates ”for each linear increment. If the generalized Coord-inates'are, interpreted as components of node displacement and rotation vectors, and if the pattern of triangles is. as indicated above (Fig. 4. 2) then the rank of the matrix is between 121 and 167 depending on the boundary conditions used. If for a given elastic system and loading 89 90 \\ ) m [m r _l _ Hi \I _ m z 9 (av 1 Fig. 4.1. Shell geometry. Fig. 4.2. Triangulation. 91 Fig. 4.3. Node coordinates. 92 the incremented solution does not' undergo a region of instability, then it may be possible to obtain acceptable results by using only two to five linear increments; however, if the solution does undergo an unstable region then it is unlikely that acceptable results can be obtained with less than 20 to 30 linear increments. A pattern consisting of a much smaller number of triangles cannot adequately describe the elastic properties of the class of structures being considered. To deal with problems having a region of instability would require an excessive amount of computer time (approximately 10 times the amount used) esPecially for the computer program deve10p- ment which required a considerable amount of testing. For this reason, the degree of freedom of the system is reduced by inter- preting the node displacements and rotations in terms of a smaller set of unknown parameters. The explicit form of the interpretation (equations of constraint) is presented later in this chapter. 4. 2. Coordinate Systems In accordance with the discussion in Sec. 2. 7 three types of right handed orthogonal coordinate systems are used, 1. e. , the general, node, and member coordinate systems. The general and node coordinates systems are fixed, and the member coordinate are'redefined after each linear increment. The O-general coordinates l x ,x2, x3 are defined (see Fig. 4. 2) so that the x3 axis lies on the axis of the shellazrl is directed upward, the x1 axis is directed radially 93 outward and intersects the shell axis and bottom middle node point, and the x2 axis is perpendicular to the x1 and x3 axes. The 0. -noc:e coordinates 23d),z(2a),z(':).are defined (see Fig. 4. 3) so that the 2(0) axis is parallel to the shell axis, the 2:0) axis is directed outward intersecting the shell axis and the c1 -node point, and the Z . . . 1 ,3' Z(G)ax1s 18 perpendicular to the 2(a) and 2(a) axes. Note that the node coordinate systems of node points on a meridian curve have their axes directed in identical directions. The member coordinate systems are defined in accordance with Sec. 2. 7. 4. 3. Reduced Set of Generalized Displacement Parameters For convenience each node point is designated by the integers 0., [3 (Fig. 4. 4). Since the axes of node coordinates associated with a common meridian curve are parallel, we distinguish them by using only the first integer of the associated node point. The following node displacement and rotation components are with re8pect to the node coordinate systems . 1 2 3 1 2 3 (a0. a)’ Au(c113.1)’Au(1113,<1)A9010.<1)’Ae(c111.<1)”"9(c113.c1) Au (4.1) These variables are related to a smaller set of generalized ~ cooridinat'es ‘ as follows: 1 1 2 . l ,- Au(051 a) = Avm) + Ava” 31n[ 2 (a -3) 11.] 3 l - (0.0,0) = AV(5) cos[-2-(a-3)Tr] 3 4 5 . l Au((1 5, c1) = Avm) + Av(m51n[-2—(o. -3) TI' ] 94 l 6 l . A6015: (1) = Avm) cos[:(o-3)1r ] 2 7 8 . 1 A0 (913: (1) : Avm) + Avm) 8111):: ((1-3) 1r. ] 3 _ 9 1 . A9(a(3,a) — AVG-3) cos [-2—(11-3) 1T] (4. 2) By using (4. 2) we in effect relate the 30 discrete variables of the five node points of a row to nine variables. This in the discrete sense implied by (4. 2) constrains the strip to a sinusoidal type of variation. The boundary conditions along the two meridian edges are automatically satisfied when (4. 2) is used. We designate by T _ 1 2 3 ”(110. a) ‘ [ “(1113.11) Au(c111. a)’ Au(c111, a)’ 1 2 3 A90113. a)’ M 0113.01) A9 («113. 1)] (4'3) T __ 1 2 3 4 9 AV(11) ‘Mvmr A"(13) A"(13) A"(13) A"(3)1 (4'4) )- -1 1 sin[-£—(o.-3)fl.] 0 0 0 0 0 0 0 0 cos[ 12((1-3).TT] 0 1 s1n[1§(e-3)r] CR :0 cos[%(043)i17] (.9) 0 1 s1n[1§(d-3)11[ 0 cos[%(o.-3)Itri] (4.5) Then (4. 2) can be represented in the form AU = c AV (4.6) (am) Rm ((3) 95 1,12 2,12 3,12 4,12 1,11 2,11 3,11 4,11 1,10 2,10 3,10 4.10 1.9 2.9 3.9 4.9 1,8 2,8 3,8 4,8 1,7 2,7 3,7 4,7 1,6 2,6 3,6 4,6 1,5 2,5 3,5 4,5 184. 2,4 3,4 4,4 1,3 2,3 3,3 4,3 1,2 2,2 3,2 4,2 1,1 2,1 3,1 4,1 Fig. 4. 4. Ordering of triangular members. 96 Let T _ T T T T AU ‘[AU<11,1)'AU(21,2)' A”(51..5)’AU(1.2,1)’ T T T A = A ..., V [ V“), AV(7)] .., A T ”(57. 5)] (4. 7) (4. 8) The generalized forces associated with AU , AV , AU, and AV (0-13. C1) (F3) are re8pectively cbsigiated by A Fe. 3.11) A (13) elements of AP [32 if i is even 51 < [32 if i is odd (4.19) Substituting (2.102) into (2.109) and using the notation of this section we obtain (a. (a )' )' (a ) z (0) 1 l (0) 2 , (0) 3 Mk [Mk1 Doc) D(0) .Mkz D D<0) .Mk3 Doc) D(0) ] P '1 AU("‘1‘31 ' “1 ) Auwz"? “2) AU ((13133: 0-3) J .. (4. 20) From (4. 6) and (4.19) we obtain ...AU -1 r—C I O - p—Av — (“151' “1) Rail): ($1) I . I (02) L _J I Substituting (4. 21) into (4. 20) we obtain 101 (a ) (a > (a ) (0) 1 (0) 3 9(0) 2 1:“le D(1013(0) CR(G)+M1<3D(1<)D(0) CR R( D(k)D(0) CR 02)] 1(1):3 1 WI) Vail) (4.22) This equation has the form AA = [M M ]_AV 2 (4.23) 1‘ Vkl sz “31) mag) The shell stiffness associated with triangular member k is then given by K = FMT r‘ [M . M ] =M 1‘ M 1‘ vkl k Vkl, sz Vk 1‘ Vk T ka2 if j is odd Kkergq I“kn‘dv Mv1=M$Pka k2 k2 kl k k T kal if j is even h ..J (4. 24) We partition the stiffness matrix in (4. 24) into 9x9 submatrices in the form r- H Kk, L l Kk,1,’ 2 Kk == (4. 25) Kk, 2, 1 Kk, 2,__2J 102 The shell stiffness matrix has the form where r- -1 K11 K12 K21 K22 K23 K K K K = 32 33 34 . K43 K44 ' . . ~ ' , K13,13. K13,14 . K14,13 K14,14 L .4 (4.26) K11= K1,1,1+K2,1,1+°” +K8,1,1 T K12 ‘ K21" K1,1,2“K2.,1,2.+”° +K8,1,2 K22 K1,2’Z+...+K8,2,2+K9,1,l+...+K16’l’l _ T .. K23 " K32 ‘ K9124”“ +K16,1,2 K33 ‘ K9,2,2"'°' +K16,2,2"K17,1,1“"' +K24,1,1 K = K 14,14 41,2,2+"°+K48,2,Z (4. 27) 4. 6. Shell ImPerfection The meridian curve of a shell of revolution can be described by the two functions rm) and z(n) (Fig. 4.1b) where n is a parameter. In the discrete sense used here these functions become 103 rm) = r“), r(z), ..., r(7) (4.29) Z 2(:3) = 25(11' ”'(21' (7) The components of node points relative to the O-general coordinates are then given by l .9. "(mm = "(131 MW 4%] 2 . ’L 3 _ “(mm ' ”(m ‘4'”) The shell surface is assumed to have a slight imperfection with n circumferential periods. The form of the imperfection is taken so that the node components of the position vector of the undeformed surface have the form 1 _ , , .1. xm, 1'3) = (rm) + Egm) sm[ (a-3)—E—]) 6031(0-3lfil’l 2 . L . "(mm = ".(m ”gun sm[“"3)%']’ “““'”f; 3 . .' where 3(a) : g(1)' 00., g7 (4032) called the imperfection function ani scalled the imperfection constant are specified. 104 4. 7. Boundary Conditions and Loading The generalized displacement parameters of AV“) and AV”) describe the deformation respectively along the top and bottom edges of the half period strip. In order to interpret the various types of boundary conditions we give the physical inter- pretation of these parameters below and a graphic description in Fig. 4. 5. v?” = uniform radial displacement 2 . . . . v0) = 8111 varying radial displacement v(31) = cos varying tangential displacement v?” = uniform vertical displacement 5 . . . . . v“) 2 sm varying vertical displacement 6 _ . . . vu) — cos vary1ng radial rotation 7 . . . V“) = uniform tangential rotation 8 . .. . . V“) = Sln varying tangential rotation 9 _ . . . v”) - cos varying vertical rotation (4. 33) A similar interpretation holds for the elements of AV7 . In order to tie down the shell we always required 4 Va) —- 0 (4. 34) Four types of boundary conditions are considered. These boundary conditions and their requirements on AV“) along the bottom edge are as follows: 105 Dis plac ements R otations shell axis I : ; E 6 / v . ‘7 . I / ' / / / , 8 ‘7 v L, Fig. 4. 5. Distributions associated with generalized coordinates v1. 106 Free 4 _ V(1) - 0 Hinged v1=v2 =v3 =v4 =v5 =v6 =v9 =0 (1) (1) (1) (1) (1) (l) (1) Fixed vl =v‘2 =v3 =v4 =v5 =v6 =v7 =v8 =v9 =0 (l) (l) (l) (l) (1) (1) (1) (l) (1) Symmetry 4 5 6 7 8 vm ‘ Va) z "(1) = Vm ' Va) z 0 (4.35) By symmetry boundary conditions we imply that the bottom edge lies on a plane of symmetry with respect to the resulting deformations. The type of load used is limited to an axial load. For convenience we use axial displacement increments so that v37) is always specified. The boundary conditions for the t0p edge are then as follows: Free no zero components Hinged vl =v2 =v3 =v5 =v6 =v9 =0 (7) (7) (7) (7) (7) (7) Fixed vl =v2 =v3 =v5 =v6 =v7 =v8 =v9 =0 (7) (7) (7) (7) (7) (7) (7) (7) Symmetry 5 6 7 8 Vm = Vm z Vm = Vm ‘ 0 (4.36) 107 The axial displacement increment is taken in the form Av4 = p('r ) 5 (4.37) (7) k where 6 is a specified constant called the displacement increment constant, and p (7k) is in general assigned a different value for each linear increment and called the displacement increment function. 4. 8. Computer Pregram A computer prOgram was written for the class of problems with a geometric configuration, boundary conditions, and triangulation as described in this Chapter. It permits a fixed normal surface load and an incremented average axial displacement of the top edge. It is used to obtain resultant axial force versus average axial displace- ment curves, and displacement curves. The program was written in 3600 Fortran source language for the Control Data 3600 computer at Michigan State University. This language contains the features of Fortran-63. The program was written so that a minimum of input data is needed. This data consists of fifteen computer cards describing the geometry, material properties, boundary conditions, surface load, and magnitude of axial displacement increment. Computer program details are given in Appendix E. V. COMPUTATIONAL RESULTS AND CONCLUSIONS 5.1. Some General Remarks This investigation was primarily directed at the development and explicit representations of arbitrary triangular member bending stiffness matrices, and to determining their applicability for geometrically linear and non—linear plate and shell problems. Two bending stiffness matrices were obtained. As already indicated, the displacement modes (3. 44) for one of them satisfies compatibility to a high degree (see Sec. 3. 7) whereas the other (D. 2) relaxes slope compatibility along the edge of adjacent members but appears to give better numerical results. The importance of the coupling between membrane and bending behaviors, particularly in geometrically non-linear problems, is of course contained in the shell formulation and reflected in some of the numerical results, but membrane behavior was not a major point of consideration. Very good comparative results were obtained; however, the numerical results are in no way adequate for the purpose of drawing any general conclusions. For this purpose the solution of problems having regions with various forms of singular behavior and the numerical examination of convergence associated with increased refinement would be very useful. Most of the numerical results were obtained from the computer program given in Appendix E and based on the formulation given in Chapter IV. The prOgram performs almost all interpretive and 108 109 computational Operations internally from a relatively small amount of input data necessary for the description of the problem and for both linear and non-linear problems. 5. 2. Linear Results In order to obtain an approximate estimate of the accumulative roundoff error due to all sources in the sequence of computations, an axially loaded perfect cylinder (Fig. 5.1) was solved. The problem results in a uniform membrane stress state which the modeling can describe exactly and yield results that are exact to within computational erros. The results are given in Table 5.1. The accumulative effect in non-linear problems is contained in these results since 25 increments were performed. Although the cylinder is initially perfect the accumulated round- off error results in a fictitious imperfection. This appears suddenly in the 20th increment and accounts for part of the reduction, in the 20th and 25th increments shown in Table 5.1. The imperfection results in both circumferential and longitudinal waves. The longitudinal waves are, however, too long, due to the node distribution, and consequently the sharp reduction in axial load associated with the bifurcation phenomenon is not revealed. The displacements along a radial line for the linear problem of a concentric plate3 with fixed boundary conditions and axial load are given in Fig. 5. 2 and Table 5. 2. Results were obtained for both bending stiffness matrices. Both matrices gave good results; however, as already indicated, the matrix given in Appendix D gave particularly Fig. 5. l and Table 5.1. 110 u :r I ‘11— V I '1 _11— ’7 1 . 0. 01""! ‘- 6" _J I - +1— I .1 __JL—JL L 3 -2 6 Pa (x10 ) ur x10 a V x10 - Exact C omputed Exact C omputed 1 9. 4248 9.389 .15 .1500 2 18.8495 18.73 .30 .2994 3 28. 2743 28. 03 . 45 . 4484 5 47.1238 46. 48 .75 .7455 10 94. 2477 91. 80 1. 50 l. 482 15 141.3715 136.0 2.25 2.212 20 188.4954 177.4 3.0 -- 25 235. 6192 212.4 3.75 .. Comparison for a uniform stress field. 111 accurate results, that are in fact to within the roundoff error of the c omputations . Some useful results can be obtained by introducing the constants c1 and C2 which we respectively call the primary and secondary bending constants. We use these constants to alter the 6x6 member bending stiffness matrix I"i as follows F of Pi clczri ll 12 2 (5.1) (11 czri czri 21 22 where the F. are the 3x3 submatrices of I". . If c = c = l. 0 1k! 1 l 2 then the matrix is unaltered. If 0 < cl < 1 then the stiffness associated with the primary bending modes is partially nullified and if 0 < c-z < 1 then the stiffness associated with the secondary bending nodes is partially nullified. The node displacements along a radial line for various values of 01 and c 2 of the concentric plate shown in Fig. 5.1 are given in Table 5. 2. These results clearly show the dominant role that the primary bending modes have on this solution. For c = l. 0 and 1 c2 = 0. 6 the stiffness associated with the secondary bending modes is essentially reduced by 64% but the resulting increase in the maximum displacement is only 2. 2%. From the definition of the constants it follows that for c = 1 02 = 0. 6 the displacement will increase by l ( 2 - 1) 100% = 177% (5.2) (0.6) Table 5. 2. 112 Plate (Fig. 5.2).and‘Percent Error. Displacements Along a Radial Line for a Concentric Radius of Node 3. 0 3. 5 4. 0 4. 5 5. 0 5. 5 6. 0 Exact Displacement 2. 852 2. 597 2. 011 l. 307 0. 6525 0. 1796 0. 0 Bending V 7 iiatrix, C1 C2 Table 13.1 1.0 1.c> 2.854 2.606 2.020 1.315 0.6573 0.1820 0.0 ° 0.06 0.35 0.44 0.57 0.73 1.32 Table 3.3 1.0 1.0 2.752 2.516 1.949 1.265 0.6276 0.1694 0.0 3.51 3.10 3.06 3.24 3.97 5.68 Table 3 3 1.0 0.8 2.775 2.536 1.966 1.277 0.6344 0.1723 0.0 ' 2.70 2.33 2.75 2.36 2.77 4.06 T;b3le 1. 0 0. 6 2. 814 2. 569 1. 992 1.295 0. 6465 0-1780 0- 0 ° 1.34 1.45 0.97 0.92 0.92 0.89 T390313 1.0 0.5 2.852 2. 601 2.016 1.314 0.6582 0.1836 0.0 ' 0.02 0.14 0.25 0.47 0.87 2.23 T333” 1.0 0. 62 2. 806 2. 563 1. 986 1. 292 0. 6441 0.1768 0. 0 ' 1.61 1.30 1.25 1.22 1.29 1.56 - T331316 0.992 0.62 2.854 2.607 2. 020 1.314 0.6552 0.1800 0.0 ‘ 0.07 0.38 0.44 0.50 0.41 0. 22 113 E l\ \ \\ \\XR r P 60. rigid-'1' 7 / / 0‘. l'Lh/F" )/ III” II ”III — — Bending matrix Table 3. 3 used 2 Bending matrix Table D.l used and exact solution /7' O / H Vertical displacement, x10 / / / / 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Radial distance from shell axis, in. Fig. 5. 2. Vertical Displacements Along a Radial Line for a Concentric Flat Plate with Fixed Boundary Conditions. 114 5. 3. Non-Linear Results Some numerical results for the large deflection problem of a concentric plate (Fig. 5. 3) and two shallow conical shells (Fig. '3 5. 4, 5. 5) were obtained. These results are compared with the numerical results of Newman and Reiss30 and the experimental results of Almen and Lazyle. 29 Free boundary conditions were assumed and the bending stiffness matrix given in Table 3. 2 was used for all three problems. The results obtained in general compare favorably. For one of the conical shells, both the numerical results of Newman and Reiss and those obtained have a noticable discrepancy with the experimental results of Almen and Lazyle. This could be due to the presence of a partial constraint (friction) on the boundary when the experiment was performed. To demonstrate this, numerical results were obtained for a variety of boundary conditions (Fig. 5. 6). It can be seen from Fig. 5. 7 that for hinged boundary conditions the maximum axial load just before instability is approximately 12 times the value obtained for free boundary conditions. A relatively small radial constraint could account for most of the 10% discrepancy present. Results were obtained for the large displacement problem of the concentric plate with the bending constants c 1 = 0. 992 and C2 = 0. 62. The curve for these results is not shown since there is no discernible difference between it and the curve given by Almen and Lazyle. 115 0. 020" mm | 0.375" | 1. 0" r: :1 Almen and Lazyle test — —- — Present solution —' "" Linear solution 9 8 / l / 7 / / I 6 ‘ I 5 . / fl 3 4 we“ , «3 / / o 1-1 3 /’/ l/ 2 / 1 O / 0 1 2 3 Deflection, 0.01 in. Fig. 5. 3. Load-Deflection Curve for a Concentric Flat Plate. 116 0.139" I 0. 238'I 7 I E Q” I 12. 766" rL Almen and Lazyle test 0 Newman Reiss solution — — — Present solution 9 8 7 6 I ’ 7%.." ’1 53' 5 o o -1 4 13" 8 3 ._l 2 l 0 0 l 2 3 4 Deflection, 0.1 in. Fig. 5. 4. Axial Load-Deflection Curve for a Shallow Cone. Load, 100 lb. 117 0.0 M ‘0.25 N 98 L MM J 10. 25" VI Almen and Lazyle test 0 Newman Reiss solution — — — Present solution Fig. 5. 5. 1 2 3 4: Deflection, 0.1 in. Axial Load-Deflection Curve for a Shallow Cone. 1000 lb. Load, @@®o ll 8 Boundary Conditions top and bottom free top hinged and bottom free t0p and bottom hinged top and bottom symmetry ll 10 \ID )\ \hm O ”A? @/ if). 63 -2 \P -3 -4 \J Fig. 5.6. Axial Load-Deflection Curve of the Shallow Cone Shown in Fig. 5. 5 for Various Boundary Conditions. Vertical displacement, x-lO.l in. Radial displacement, xio'3 in. 31. an O -1 -2 -3 -4 -5 l 2 3 4 5 6 7 \ ' \ \ 5.1250 4. 7396 4.3542 3.9688 3. 5833 3.1979 2. 8125 .7) 11sz (12 f. u r 119 Node Point No. Radial distance from shell axis, in. Fig. 5. 5 with Free Boundary Conditions. 8 6 x 4 \\\ 2 \\\ -2 \\ -4 \\ -6 . 5. 7. Displacements Along Radial Line for the Cone of Node Point No. 120 1 3 4 5 (' 7 1 .3 0 "$2 -1 Q 5 \ \ 8 -2 \ \ .9. in“ s Q -3 ; . '5 \ .2 1:: _21 \ § 1 \ ® -5 5.1250 4.7396 4.3542 3.9688 3.5833 3.1979 2.8125 Radial distance from shell axis, in. -2 / Radial Displacement, 3210'3 in. -3 / Fig. 5. 8. Displacements Along Radial Line for the Cone of Fig. 5. 5 with Hinged Boundary Conditions. 121 Node Point No. (see Fig. 5. 7) 1 2 3 4 5 6 7 l . in. H \\ z \ Vertical displacement, x10 5 5.1250 4.7396 5.342 3.9688 3.5833 3.1979 2.8125 Radial distance from shell axis, in. -2 -4 Radial displacement, x10"3 in. -6 Fig. 5. 9. Displacements of Node Points Along Radial Line for the Cone of Fig. 5. 5 with the Various Boundary Conditions of Fig. 5.6. 122 5.4. Conclusions The usefulness of stiffness matrices associated with arbitrary triangular plate elements is self evident. With the aid of these matrices and modern computing facilities we can formulate and solve many shell problems. The ease with which many geometric configurations can be modeled by triangular plate elements and the inherent ease with which most boundary conditions can be imposed makes this procedure very attractive. The two bending stiffness matrices obtained for arbitrary triangular plate elements appear to give results of sufficient accuracy for many applications. To ascertain the full flexibility of this procedures, however, will require more numerical results than those obtained, especially for problems with various forms of singular behavior. This investigator was particularly interested in the bending stiffness matrix given in Table 3. 3 because it is associated with displacement modes that have regions of constant curvature and twist, and can consequently be more readily adapted to problems with non-linear stress-strain relations. The modes associated with both bending stiffness matrices can be conveniently used to include the effects associated with thermal strains that vary through the thickness of the shell. 10. ll. 12. 13. BIBLIOGRAPHY Langefors, B. "Analysis of Elastic Structures by Matrix Transformation with Special Regard to Monocoque Structures. " Journal of the Aeronautical Sciences, Vol. 19, No. 7, 1952, pp. 451 -458. Argyris, J. H., and Kelsey, 8., Energy Theorems and Structural Analysis, Butterworths, London, 1960. Timoshenko, S. , and Woinowsky-Krieger, S. , Theory of Plates and Shells, McGraw-Hill, New York, 1959. Mushtari, Kh. M., and Galinov, K. Z., Non-Linear Theo_ry of Thin Elastic Shells, Tatknigoizdat, Kazan, 1957. Available in English as NASA-TT-F62, 1961. Novozhilov, V. V. , The Theory of Thin Shell_s, Noordhoff, Groningen, 1959. Langhaar, H. L. , Energy Methods in Applied Mechanics, Wiley, New York, 1962. Courant, R., and Hilbert, D., Methods of Mathematical Physics, Interscience, New York, 1953, pp. 175-176. de Veubeke, B. F. , "Upper and Lower Bounds in Matrix Structural Analysis, " AGARDograph, 72, Pergamon Press, New York, 1964, pp. 115-201. Prager, N., and Synge, J. L., "Approximations in Elasticity Based on the Concept of Function Space, " Quart. Appl. Math., Vol. 5, No. 3, 1947, pp. 241 -269. Hrennikoff, A. P. , "Solution of Problems in Elasticity by the Framework Method," J. Appl. Mech., Vol. 8, No. 4, 1941. McHenry, D. , "A Lattice Analogy for the Solution of Stress Problems, " J. Inst. Civil Eng., December 1943. Parikh, K. S., and Norris, C. H., “Analysis of Shells Using Framework Analogy, " World Conference on Shell Structures, Oct. 1962, pp. 213-222. Turner, M. J., Clough, R. J., Martin, H. C., and TOpp, L. J., ”Stiffness and Deflection Analysis of Complex Structures, " J. Aeron. Sci., Vol. 23, No. 9, 1956, pp. 805-823. 123 14. 15. l6. 17. 18. 19. 20. 21. 22. 23. 24. 124 Argyris, J. H. , ”Recent Advances in Matrix Methods of Structural Analysis, " Progress in Aeron. Sci. , Vol. 4, Macmillan, New York, 1964. Zienkiewicz, O. C., and Holister, G. 5., Editors, Stress Analysis; Clough, R. W. , "The Finite Element Method in Structural Mechanics, " Ch. 7; Zienkiewicz, O. C. , "Finite Element Procedures in the Solution of Plate and Shell Problems, " Ch. 8; de Veubeke, B. F. , "Displacement and Equilibrium Models in the Finite Element Method, " Ch. 9; Wiley, New York, 1965. Melosh, R. J. , “Basis for Derivation of Matrices for the Direct Stiffness Method, " AIAA Journal, Vol. 1, No. 7, 1963, pp. 1631-1637. BOgner, F. K., Mallett, R. H., Minich, M. D., Schmit, L. A., ”Development and Evaluation of Energy Search Methods of Nonlinear Structural Analysis, " AFFDL-TR—65-ll3, Wright-Patterson Air Force Base, Ohio, Nov., 1965. Melosh, R. J. , "A Stiffness Matrix for the Analysis of Thin Plates in Bending, " J. Aerospace Sci. , Vol. 28, 1961, pp. 34-43. Adini, A. , "Analysis of Shell Structures by the Finite Element Method, “ Ph. D., Dissertation, California Univ. , Berkeley, 1961. Tocher, J. L. , “Analysis of Plate Bending Using Triangular Elements, " Ph. D. , Dissertation, California Univ. , Berkeley, 1962. Reissner, E. , "On a Variational Theorem in Elasticity, " J. Math. Phys., Vol. 29, 1950, pp. 90-95. Reissner, E. , "On a Variational Theorem for Finite Elastic Deformations," J. Math. Phys., Vol. 32, 1953, pp. 129-135. Argyris, J. H. , "On the Analysis of Complex Elastic Structures, " Applied Mechanics Reviews, Vol. 11, No. 7, July 1958, pp. 331 -338. Turner, M. J., Dill, E. H., Martin, H. C., and Melosh, R. J., "Large Deflections of Structures Subjected to Heating and External Loads, " J. of Aero/Space Sci. , Vol. 27, No. 2, Feb., 1960. 25. 26. 27. 28. 29. 30. 31. 125 Gallagher, R. H. , and Padlog, J., “Discrete Element Approach to Structural Instability Analysis, " AIAA Journal (TN), Vol. 1, No. 6, June 1963, pp. 1437-1439. Denke, P. H. , "The Matrix Solution of Certain Nonlinear Problems in Structural Analysis, " Journal of the Aeronautical Sciences, Vol. 23, No. 3, March 1956, pp. 231-236. Lansing, W. , Jones, I. W. , and Ratner, P. , "Nonlinear Analysis of Heated, Cambered Wings by the Matrix-Force Method, " AIAA Journal, Vol. 1, No. 7, 1963, pp. 1619-1625. Lansing, W. , Jones, I. W. , and Ratner, P. , “Nonlinear Shallow Shell Analysis by the Matrix Force Method, " NASA TN D-1510, 1962, pp. 753 -761. Almen, J. O., and Laszle, A., “The Uniform-Section Disk Spring," Trans. ASME, Vol. 58, No. 4, 1936, pp. 305- 314. Newman, M., and Reiss, E. L., “Axisymmetric Snap Buckling of Conical Shells, " NASA TND-1510, Dec. 1962, pp. 451- 462. Pian, T. H. H. , ”Derivation of Element Stiffness Matrices by Assumed Stress Distributions, " AIAA Journal, (TN), Vol. H, No. 7, July 1964, pp. 1333-1336. Appendix A. OBLIQUE COOR DINATES A. 1. Some Properties of Oblique Cartesian Coordinate Systems Consider the oblique Cartesian coordinates designated by Lk and the rectangular Cartesian coordinates designated by xk. These coordinate systems are related as shown in Fig. A1. The coordinates Ll, L2 and x1, x2 define a common plane, and the coordinates {,3 and x3 are normal to this plane. The angle {-.I 0 4.2 (= x1 0 L2) is designated by LIJ . Unit base vectors of the Lk and xk coordinates are respectively designed by Ek and j'k . The associated reciprocal base vectors are 3k and jk (= jk), and are defined by the properties 'e'k-El = 1fork=1 (A.1) 'j—k°j'£ = Ofor kfil It then follows that the magnitudes of the base vectors are IE I = l‘e' | - 1 1 2 (A.2) I21) - I'ézl — 1 _ _ sian 71 -_2 1" 1- IJI=IJI=IJ11=IJZI =1 (A...) A position vector '1: can be represented in any me of the following forms - --k k:- 7k r=§ek=§ke ZXJk=XkJ (A.4) where Lk are called contravariant components of the vector, Ck are called covariant components of the vector, 3k are called contravariant 126 127 €31 l 51 ‘L ml N Fig. A1. Coordinates and Base Vectors. Fig. A2. Base Vectors. (b) (bl 128 base vectors, and 3k are called the covariant base vectors. The coordinates Lk and xk are related by xl=§1+cos¢§2 x2: +sin¢§2 x3: +§3 In tensor notation (A. 5) can be represented in the form k I k x =c£§ The inverse to(A.6) is designated by Lk = d]; x1 where the matrices [cg] and [(1%] are related by k k -1 Id, 1 = is, ] The covariant components are related by xk = d §k=c Substituting (A. 6), (A. 7), (A. 9). (A. 10) into (A. 4) we obtain k- la- k 14- ; ek : X J! = g Ck]! -k *1 ke-I Cke = x1] = Lkdfj (A. (A. (A. (A. (A. (A. (A. (A. 5) 6) 7) 9) 10) 11) 12) 13) Eq. '3 (A. 12) and (A. 13) can be satisfied for components of all position vectors 7 only if 129 .. I ... = . 4 — _ k-J ' ek — d1 J(x) (A. 15) From (A. 5) and (A. 9). we obtain 1 l 1-1 11 OT c1 c2 c3 cos I11 1 g 2 2 2 _ . 3C1 c2 c3 — 0 51111.)! 0 (A.l6) to? c: c; 0 0 1 ‘_ _J _. ._ r— — )— —) l 1 1 cos d1 d2 d3 1 - sin L11 0 2 2 2 _ 1 d1 d2 d3 — 0 sin L); 0 (A. 17) 3 3 3 d1 d2 d3 0 0 l For the oblique and Cartesian coordinates with covariant base vectors as shown in Fig. A2 the transformation relationships are as follows: _E - 1203 0) sin (0 0-1 _7— 1 1 1 J1 32 = cos (.02 sin (1)2 0 j-Z (A. 18) E 0 0 1 ; 3' 3 ’ 3 a _ ._ _J __ ... '1- . - . - 1 J1 Fsm 102 81111;.)l 0 e1i .- _ l — J2 - (costsl sin 0.12 - sinw1 cos 002) C08 “)2 C03 (01 O 62; 3'3 0 0 1 E3; 130 A. 2. The Plane Stress Strain Energy Expression Referred to Oblique Cartesian Coordinates The covariant components of stress, strain, and displacement when referred to the x1, x2 coordinates are respectively designated by gij , Sij and ui; and when referred to the L1, CZ coordinates are reapectively designated by O'ij, eij’ and WI. The corresponding contravariant components are designated by 9;”, g”, 111, 0'”, 61‘], and w1 . The stress-strain relationships in plane stress are 11 _ _ E 51 ‘ 511‘ 2 (911+ ”$22) 1-1/ 22 _ _ E E ‘ £22 ‘ 1_ 2 (522 +511) (A10) 12 E E 0 = g = —— g = —- (l-v) g «- 12 l+v 12 (l-VZ) 12 The strain displacement relationships are _ l m ‘ 2 (u1:j+uJ.1) 1 au'i Bu. =7(—_—.—+ i) (A.21) BXJ 8x and — 1— (w + ) 6ij ‘ 2 i,j “3,1 1 3w. SW. = §(—-}-+——J—) (A.22) 82;} at,1L _ k 2 ij ‘ C'i Ci $1.1 ““231 = 5d! 5 (A.24) 131 Eq. (A. 24) is equivalent to the following matrix operations _ I " "- I - P I ‘- 1‘ 1 cos il- 511: E12 1 I '0 611 I E12 1 :'sin0 I l "- _. 1- -- : "Eogl-T—f- --__.|---_ "—-I ------ $21 I E22 'sin0 I sinq; 621 I 622 0 I L. f ._ _ 1 .. _ I __ _ l .1 P E 1 cos 6 1 e 11 I sinL|J 11 sinlp 12 = ————————————— + ——————————————————————— 2 3|, 1 5|, '2ng.). 611 silnq. €12I cos2 611+ 1 " 2cos , sin 1).) sin L): sin I): .. I (21.25) The plane stress strain energy for a thin plate of thickness t and middle surface S is given by 1 ij W = — t 0' 6.. I1 2 ISI“ ’1] (A. 26) Expanding (A. 26) and substituting (A. 20) into the resulting expression we obtain 2 2 1 -22) " 2(1"’)(§~11§22 "£12”dx dx (A. 27) Substituting the strain relationships given by (A. 25) into (A. 27), putting u = 2(1-v), and noting that dxl. dxz = sin 1); dgl déz we obtain 5t 1 2 WI = 2 . H .2 (611‘2C054’612J’622) l 2(1-1/ )Slnlp S Sln 1).: 2 1 2 "“ (611622 '512 ”‘16 dg Eq. (A. 28) is the desired relationship. (A. 28) 2 132 A. 3. The Thin Plate Bending Strain Energy Expression Referred to Oblique Cartesian Coordinates The covariant components of bending moments and curvatures when referred to the xk coordinates are designated by in” and 5%., and when referred to the Lk coordinates are designated by m.. J 13 and Kij . The associated contravariant components are designated identically except that the indices are raised. The bending moment-curvature relationships in rectangular coordinates are mll : m -J(K + VK ) v- ...11 -ll 22 22 m ‘ Jr5322 ‘ "015224r ”511) (11,29) 12 _ .. - r3212 - .O’U-VHSIZ where ’0/ = E 133 12(1 -V 2') The curvature displacement relationships are 8 2W3 51.: 113 .. = -—-—i--—-.- (A.30) J 1 1.1 3X 3x3 3W3 Ki. 2 W31. = ——i—'—-:' (A.31) The components of the curvature tensors are related by 1 .k Kij _ c.i cj 33k! (A.32) _ k 1 :ij - di dj Kkl (A.33) 133 Noting the similarity between (A. 24) and (A. 33) it follows directly from (A. 25) that PK“ K “ i- K I-MK + 1 K .1 «11 -12 11 : sinqu 11 sinip 12 = ' Z _ cos 93 l l cos 1 _ 2cosgg .512 £271 SinLIJ Kll sinL); K12l sin .0 "11 + sinzlp K22 simp “12 _ _. I .. (A.34) The strain energy W due to bending stresses for a thin I 2 plate of thickness t and middle surface S is given by 1 w =—tffmin I2 2 s" 4]. dx1 dxz (A.35) Substituting (A. 34) into (A. 29) and then substituting the resulting expression into (A. 35) we obtain - __1__ 2 2 l 2 WI .- 2 Slnhll ff{ . 2 (K11 '- ZCOS¢K12+K22) ' “(K11K22 'K12)}d5 d; 2 S sin 01 (A.36) Substituting (A. 31) into (A. 36) we obtain 2 2 2 3 w 8 w 3 w w =—¢g:—f { l ( -2c08¢ +-'—73 )2 1‘2 281an sf “1124‘ 61213,; a{‘13,} 3:202; 2 2 3 w 82w 8w 2 . 3 3 3 1 2 -M [7— ' ———g -(-———).]1d§ d4 861061 8628; 2 . 8802 (A. 37) This is the desired relationship. Appendix 13. CONSTRUCTION OF TRANSFORMATION MATRIX (0) D1(i) This matrix relates the base vectors J(i) of member i with node points 1, 2, 3 to the base vectors J(0) of the general coordinate system and has the form 3“1(i) 720) 3'30) 7 _ .. b11 b12 bl3 J1(0) b21 b22 ID23 3 2(0) (3' 1’ b31 b32 b33 J3(0) b _J h. ... where the elements of the coefficient matrix are Obtained from (2. 101) and are given below. 0‘ I 11 12 U ll 13 21 22 23 31 32 33 _1[x1 -xl] 1'; <2) <1) 1 [x2 _x2 ] I; (2) (1) -1- [x3 -x3 l 13 (2) <1) 1 1 1 2 1 EU? ' aL2 x(l) ' a1"(2)+ ‘3 "(3)] 1 [-axz-a x2 ”22] _—I—2A 3 2 (1) 1 (2) 3 "(3) 1 [-a x3 -a x3 +12-x3] "'2A"'13"' 2 (1) 1 (2) 3 (3) 3 2 3 3 2 2 3 i.[x2 x3 -x x2 +x x -x x +x x 2A (1) (2) (1) (2) (2) (3) (Z) (3) (3) (1) 3 2 ] "‘(3)"(1) J-[x3 x1 -x1 x3 +x3 x1 -x1 1:3 +x3 x1 -x1 x3 ] 2A (1) (2) (l) (2) (Z) (3) (Z) (3) (3) (l) (3) (1) i-[x1 x2 -x2 x1 +x1 x‘2 «x2 x1 +xl x2 -x2 x1 ] 2A (1) (2) (l) (2) (2) (3) l2) (3) (3) (l) (3) (l) (B.2) 134 Appendix c. DERIVATION DETAILS FOR MEMBER STIFFNESS MATRICES C. 1. Membrane Stiffness Matrix for Triangular Member The membrane stresses of a member are associated with the inplane displacements given by (3. 42) and have the form L, w151+w222 e 1.1;{-p11).1 - (l-o.-[‘3)13>.3-[(1-.6} Ta— = "1’2: {1212 + 2[(1.3s.(3)Aa3 - (1-0-(3)Ab3]>~4 - 210.4513 + (2-30-(3)Ab3 + 8Aa1+ aAbl] 15 + 2(8-s)A.bl>.6} (c.2) 135 136 8w)l A 7,—5— -I;{ 11311” X3-2[--(1-<113)Aa2+(1'°'3mA'b2])‘4 - 2(e -p)Aa1 x5 + 2[ p Aal + aAbl + (2-2e-38)Aa2 + 8Ab2]).6} awz W: 1-2—{12X2-2[0Aa3 +(1-s-8)Ab3]>.4 - 2[(1..2d..(3)Ab3 + (a-B)Aal] A5 + 2[(3Aa1 + aAbl] 1.6} ((3.2) The deformation parameters 61, 62, 63 are defined by (3, 15) and are related to displacements by (3.16). Substituting (C. 2) into (3.16) we obtain 61 a 1111 + 2A{[ .(1.2s.(3)AaL3 + (1.61—2(3)Ab2]>.4 + [ClAa3 + (1.6...(3)A.D3]>.5 - [(1-O-(3)A.az + pAbz] 1.6} £2 = + 1212 + 2A{f-[ aAa3 + (1-6-p)Ab3]x4 -[(1-2(1-B)Ab3 + (6-0)Aa1115 + [(3.4611 + 6Ab1]>.6} (0.3) 63 £313 + 2A{[(l-s-s)Aa2 + 8Ab2]).4 - [ [BAal + aAbl])\5 +[ -(.l 1.5-216 1116} 1 2 +;—$3{36xf +54x§+5412+36113-721316+36111} 41113 122 gfeg dads: 313x3+ 3T93{-43>.2161>.3+216>. 16} 1 +—4——AZ {5412+361§+5412+3614x +3615). 7211} 2592 4 6 5 6 ‘ 6 4 ffezdads= 1—1212+4A1———3 {216x31 -216>.).} s 3 2 3 3 2592 4 35 1 +—"’——A2 {5412+54>.2+36>.2 7211 +361). +3611} 2592 4 5 6' 45 56 45 21113 1 ff 6162 map: -1313>.l>.3 +——— 339——-3= {216 x 315 - 216 1316} SI 2.41332 3 3 3 +3-—3———931 {-2161 14+ 2161113} +-3——4593{.1814-1815-3613 -7211 +36156416+3611} 4 5 21113 ffeze 3 dads: -3- 13131313 +——- -3——-593 2(216 1413 - 216 1313} Si 21113 +33—9—3 { 216 x 413 + 216 1613} 4A2 2 2 2 +3393-{-36>.4-1815 -1816+36x4>.5 -7213x3+36x4x 6} 1 2111 fS f 63 ldodB = +— 3 11131113 +—-— -3--—593 {216 1 x4 - 216 1113} SI 211113 +35——9-—3 3'{216 1513 - 216 1613 3 312—132 2 2 2 +7—593 {-1814 - 36 15 - 18 16 + 36 1415 + 361516 - 721416} (C. 5) 139 Substituting (C. 5) into (C. 4) we obtain the expanded form of F A (C.6) i . where I"11 , called the member membrane stiffness matrix, is given in Table 3. 2 and A1i = [13 .. 16] (C.7) As shown in (3. 83) the membrane and bending behaviors of the member are uncoupled so that from the principal of virtual work 6W I1 . El. : 8x3 : P11. A1. le, .00, 6 (C08) 1 J 1 1 where El are the generalized forces associated with the generalized i coordinates 1 A . l. 1 C. 2. Bending Stiffness Matrix for Triangular Member The bending stresses of a member are associated with the displacements normal to the plane which have the form EL - W33 '33 = 11(3)“? + (3(l-o-5))\8 + o.(1-11--(3)>\43 + [Lam-1911.31 + run-13mm] x10 - [<1-a-2m11-a-mA32 + [3(1 -o.-ZB)AbZ] X11 + [o.(1-20.-(3)Aa3 + (1 -o.-fl)(l -20.-(3)Ab3] X12 (C.9) 140 Substituting (C. 9) into the expression (3. 23) relating the deformation parameters €4,65,€6 to W3 we obtain 64 - 2X7 + 4(-Aa1 + Ab1))110 + ZAblel - 2Aa3X12 ‘ Ab2”‘11 + 2Ab3"12 -21 -2A 1 +4(-A a2 65 8 a1 10 66 = ' 2"9 + 2Ab1x1o ZAaZXII + 4"Aa3 + A163)"12 (C.10) From (C, 10) and Table Cl we obtain ff e4 dads: 21‘2 + 8150+ 1:1 + 1:3 - 213131+ 217113 -3139111 S]. ”3110112 “311012 f f a: dads: 212 8+1f9+ 81?1 + 133 + 213110 - 213113 - 31—119111 S]. '3‘10112 'gxuhz 1‘er (10.de 21f3+120 +1f3+81f3-219110+219111 -3139113 1 “3110112 éxllxlz éf €465 dads = 21318 + 13110 - 17113 + 21?0 + 2151+ 2119111 i ' )‘loxlz ' x8xll + X8)‘12 ' k11>‘12 2 2 4f €566 dads - 2111+ 2133 + 21819 + 13111- 13110 + 2111113 1 ' )‘llxlo ' 19113 + kc9>‘1o ' xlleo 2 2 131 6661 dadp _ 21133 +110 + 21917 +191133 - 19111+ 2133110 1 “)‘12X11'X7xlo +19111-139131 ((3.11) 141 Substituting (C. 1 1) into (3. 25) we obtain the expanded form of P A (C. 12) where F22 , called the member bending stiffness matrix, is given in i : Table 3. 3 and T _ Azi - [ks], .00, x12] ((3.13) The principle of virtual work requires 3W I2 . 223 = 3343 =Tzzi AZi 3:7, ..., 12 (C.14) where E 2 are the associated member generalized forces. i Appendix D. SOME ADDITIONAL RESULTS ON THE TRIANGULAR PLATE BENDING STIFFNESS MATRIX D. 1. Alternate Form of Member Bending Stiffness Matrix The six independent displacement configurations associated with bending and corresponding to those given by (3. 43) are taken in the form _ _ ..3_ uL — W3 e _ {gaZ ga3 >17 + ga3 ga1 x8 + ga1 ga2 19 . 3 + ga2 ga3 gbl )‘10 + ga3 gal gbZ )‘11 + gal gaZ 363 "12}8 (D.l) These displacement configurations have the same qualitative form given by (3. 43) and shown in Fig. 3. 6; however, displacement and slope continuity are not in general preserved for any triangulation. As a result the element of the stiffness matrix cannot be interpreted as satisfying upper bound requirements. This stiffness matrix appears to give better results whenever large curvature gradients are not present (see Ch. V) and is consequently included here. Substituting from Table 3. 1 into (D. 1) we obtain e = {0.317 + (S(l-o-BMB + (l-o-(B)o. X 9 + ofl(o. 43))» 10 + (3(1-0. -(3)(-1+o.+2(3))\l1 + (1-6-6)(a)(1-26-p)113 (13.2) On taking the second derivatives of (D. 2) and substituting into the expression for the deformation parameters (3. 23) we obtain 142 143 64 : -211X7 - 611(a-B)Xlo + 2I3(1-u-B)Xll- 2£3(l-o. -(3))\12 65 = - 21218 - 2116119 + 612(1-0 -2(3)Xll + 213(1le 66 = - 213139 + 21113119 - 2136111- 613(1-2a -13)113 (D.3) From (D. 3) and Table C1 we obtain 2 _ 2 2 2 2 1 2 2 _1_ 2 2 £f€4 dads _ 21117 + 311110 + 3 13111+ 3 1 113 i 4 4 2 ' 3111317111 + 3111217112 ' 312191912 2 _ 22 22 1_ 22 1_ 22 ffes dads— 21318 + 313131+ 3 13113 + 3 11110 S. 1 4 4 2 ' 3£2£3"8"12 + 3121311910 ' 3131111910 +3121“? +3—1f1f +-1-1"'1‘1Z -11 1 11 2 _ 2 £I‘6dadp"u3" 312 0 3 1 3 31910 1 3121011 WM: \ON 2 I -— 13 1142111 I I X h _ 2 2 2 1 2 2 3f E455 dadfi ‘ “1121718 + 311*7’30 ‘ 3 111313113 +2 13119 1 5 1 1 2 ' 31113110111- 31113139113 + 2 £2111 £1111 11111 122 2 2 1" '3‘2"8"11+3 23812’3’ 231112 ‘6312 _ 2 2 2 _1_ 2 2 Isf65‘6 “dB ' 212‘3‘8’9 + 31211911 " 3 12111910 + 2 13111 1 5 1 3121911112 3‘2 1111119+3 :132 21211 +51111 -1111 1 1 '6! 0'3'3912 331910 3311210 144 2 2 2 1 2 2 = — —- -1 £95664 dadfi 2131119119 + 3 1319113 3 131319111+ 2 3113 1 5 1 1 2 - _. 1 __1 1 —1 313 1113119 3 3 2*12"11+2 1130 1 2 2 2 2 2 1 -.. -_ _ -_ 1 612111 31113119+3 111217311 3£1 2"10"11 (D. 4) Substituting (D. 4) into (3. 25) we obtain the expanded form of _ 1. T WI ‘ 2 A 2. I22. A2. (13.5) 2 1 1 1 where the member bending stiffness matrix P22 is given in Table i D1. D. 2. Member Transformation Matrix for Bending For convenience the member transformation matrix M23 relating the deformation parameters 1 7, . . . , 1 12 , associated1 with bending, to the normal node displacement components and inplane rotation components is given. This can be obtained by taking the transpose of the coefficient matrix in (3. 72), deleting the first six rows, and deleting those columns associated with inplane node displacement components and rotation components normal to the plane. The transformation matrix is given by (D. 6). 145 3 .3 :36 :mvm :mr. :3 o :3 o :35 save 3 i a a :39 t. o 0 RED < u m“ WM 4 fl o o MN LIFO 4 mg IE-6 Aw IN... 3 o o o WAN mp? mm < .w. A o o o a mu - m~_q Na 4 NH : OH 146 fim m A1+Amm- _ . mmmmm + ammwm + _ <0 _ N H m N <0 .m III I Nm+~m+m$zvmqn l||l|||I~I I1 llllll .1 1.: 1 Comm: _ _ x + _ _ 3+mem+mmzm+ _ mmm+mm~mm+mm+ _ _ <0 . <0 N H m..-In <¢ N N MN _ mmm+~£z 0mm . 6+N6320~ _ my my I||l|||ll_ llllllll all! lllll g _ A1+AMMNNI 3 _ A1+Ammmm+mm~m+ L m+mmfim+ _ moddm+mmamm+ L _ <0 _ <0 <0 ......M311mw 01-31430 5.111611% _ _N N N ~ ~ ww qg N0. lllllllllll .lul'llillli llilllli llii'l < .. a. _ a. LI .3. A N 04.2-“ N 0.43- 3?. N 0.42 - E+wl _ adammvmm My _ fimvmd -m~ _ Amoummvmm mew Mm |||||||||| _lililllllfilllllllll < < < < 3.1.3.2 mm VMNML «(WC Ma- _Emgmmkflwlwwwz- ETMNNV Amumvmm N~_ .m Nu. .mww mm llllllllll LIIIIIIIIIJIIIllllllillllllrlllllllI'll < _ < < < < EN. 3 w a 0m~mJ3~+ 2 0mm: 2mm-m£.ml0I¢~oxomvh<£C0m ON «\mOZDOQ 024 mwlozu Zn wfl< mh~ZD JJmu >JJ<~X< Z< k0 WU~Fm~flwFU .AMVDO oawohvaooa0ohv>OoAOJFVZO ..0000.AN~.N~.¢0wm.“MVUO£.AnvaI.aflon0F~DZ~FZOU DmeJUIQOPIHN.XmN.FF~DZ~FZOU OMmOJUIZOFFOGI¢NomeuF<£UOk anonymG.XO-FINoX~u o®>INoX~uob>+Noxmnofl>IN.X-om>Imoxuu.¢>IN.X-.M>INoX—~.N>IN.X- ou>INox-oxon\mh2mZOQZOU F2m2w0ahmii>leOhIN~.XMN.F<2&O& awMQEIQObIQoXWN0Fahm2£>m120FhomImu.XmN.h<£&Om AwwflklZOFPOOIu~oXWN~F<£¢Ok amzonhnozou >Ga n~v~2m0U*.fivflfluafio~.Xd bonifiOnuOO Monl~on~00 JJwIm OMZQOKMOZD act GOPUM> ZO~F~wOQ mDOZ ZUOk aN.QZ.DZUow.ON hznafl on hznflfl mnhzufln ahonuuoauvaiwkvom OUF~.OO¢ME QZ.QO(UG ZOUZ—Q.O Dtma 0(OJmU.0 Odwfl UZ~>¢0 Oa<0230m no Zo_bu w>_ba«.n.~.>a .~a~2mou*.~.sz.w*.n.az~u¢mam+.1.acxau.n._.xa .Acmzz_m*.o.az_u#mnm+.n.~.Na u.n._.Na b.~un mag 00 m.~u_ mo" 00 zoeuwuawdza 4<~k~za don dehuw> zo~h~moa kumauou tzuhzou adnw JJU sziwdw Guzzu kfldbm Ooounfio~0230 Occuofiouvuo Oouflfi MuN OO roman muN 0o Ooouafiouvxw noonufi WON OD moonun MON 00 MJU>Uh~o~flZ mmm OD wJU>U ZO~F OoOla5o~0>o Ooouafiouuoo Oa0fl55o~v~0 nNN ONN mnN WON CON 159 cm hzaaa .r._us..o.~uz..z.4.>...¢¢ hzaaa nc kzaan maox<.mc hzaan o.m.xa.taauo.zoo*oo:.o+z.asym+amzoou.mzoo .z.m4.>o*.z.4.xm+_0*aJ.m.m+a¢.I.>o#.4.¢.muan.bo .m.:.>0*.4.n.mu.N.Do .N.t.>0¥.J.N.m+n~.£.>0*.J.~.0ua”.30 b.~u£ own 00 monuJ 0mm 00 mackum> zonhnmon humakou .2.J.>O+.!.J~>uxt.4.> oouut Nnn OD b.~u4 Nnn oo mackum> szzwuo+mnox«nm~0x< .ZCZOUZ~Q*UZ~>un¢.b.>o “Juou.~t.~dv>o gal—JutGladvumAuuz ~+0\.~lad.um~.u_4 m¥02.~u4 0mm 00 A~.O.m¥02.m¥w~>2~b*.wm.t.¥w0.~£.~J.OGuAJ.o .~I—J.*0laJvum~u.I ~+0\.~0~J.Um~.u~J «Jaumuuz m¥62.~u4 Gun 00 a—£.~J.¥0u.t.4.mxo atvum~unt .J.Um~u~4 m¥02.—ut con 00 myuzo~n4 com on mwuaaqu >k~J~m~mem 024 mmwzuuAFm OWUDoma tack “J.z0230*o<04mu+2003#.4.x.~01u.4.t.00 h.~nz CON OD mounJ com on hzwtdeZ~ OUP~.u~ .h.~uJ..oo~mtoaz.J.~o...¢¢ hzuan NNN Oflm O—n OOH DON 0.N\aa-OWUO~mIaN+—.0mwoum+an+—vomwouwwla—ezoumou ONQ monu— 0N6 00 Nita—uwonmuauvcmwanw 0n¢ woman Oucoo anywO—muafl+~uwbnm 00¢ hogan 00¢ 00 mNU~¢Pan.mv.m_.>a.+ «#1..an.a..xaunmw.m..xa..u.n.wo_m m.o¢*. mm:..nu.n~.Nan.~1.a..Na. + m**..mu.ma.>au.aw.-.>a.+ mm*.am1.n_.xau.an.-.xa..u.m.wo~m _. m.o¢*. m**a.mn.m~.Nac.nw.n~.Na. + m**.amw.m~.>anxnw.m~.>a.+ m:¢..mw.m~.xauans.n_.xa..u...wo_m a n a X~KF c.0030 oaoobvoooaoofiv>0o50oh020 ..novco“NuoNuocvwk.”MCUOannvmoz.aflomvm .AbomuxaoaNno0.N~.Uo “Moovmo“mono094mo“Momvwmo.NoOoN~u<\¥m\ZOIIOU Edam wzanOCDDm 02m wDZ—PZOU "UZOU+AZ.J~NCIAIoJuNfl ~0200+.I.J~>¢uatodv>fl «<20U+atodvxauatoduxa a“2030*AJ.~2ofl~ ~.~u~z car on .o._u.z our on ~—._u¢~ our on .5._.hzmzwsm aqsooz<2ak «on < xaabaa._a.._.>a*.nw.n_.xa+ m anw.n~.xa*.mw.~_.>mn.n1.m_.>a¢.mw.macxa+ . .mw.m_.xa*.av.-.>au.m1.m_.>um.an.~_.xa.*.w.m.amuan.n.mm ..an.._.~a*.nw.m_.xas._1..asxa*anw.n~.Na+ m .nw.n_.Na*.~5.~_.xan.nw.m_.xam.mw.ma.wa+ a .m1.m..Na*a.1."—.xan.mn.m~.xa*.aw.~_.Na.m.w.m.amuam.n.um ..an._~c>a*.nw.n_.Nau._n.~..Na*.nn.m~.>m+ N .nw.n_.>a*.mw.ma.Nan.nw.m_.Nat.N1.m..>a+ _ .Nw.ma.>a*..1.__.Nauamn.m_.Nam."5.._.>a.*.u.m.amu.~.n.wm .nu.n~.Na*.n.c.am+.mn.~_2Na*.w.m.aw+a.3.~_.~mm.5.m.amu.n.m.Mm .mn.n_.>a*.w.¢.am+amn.m_.>a¢.n.m.am+.an.-.>a*.n.mcamu.m.m.wm .nn.n_.xa*.5.¢.aw+.mw.ma.xa*.n.m.aw+._1.__.xa*.n.N.aMu.~.N.wm ..as.aacNac.mn.m..Na.*.w.a.ama.m.~.wm ...5.-.>at.mn.~_.>a.*.w.~.awu.m._.mm ..as.a—.xat.mw.m_.xa.*.w...awua~.~.wm x2apd oafiom0XGoaNno0oN-Uo amoovmo«Monofl04moamomuwmo“NoOoNu0<\¥m\ZO£lOU ¥U¥m¢ mzubboambm 02w ZEDFUQ AOQVZOUZ~Q.MGBon3o~03on£3oN23o~23oac~m3o50oF~ZDOoah0QI~ko G aNoflonwhmoOJowJoO¥o0¥oan.n~mQOoQZUoD!Uom~OX< otflOZZDO. m a0.N~.FU.aM000£oam.MO!ono.002.0 .anvDO onoohvooonGoF~>OoAOoF~ZG G n N u 166 nmzouu~o+z. Joym udzouua E. J.¥m 00!.Z.O+Z.<*~fi.2.4.mm~km+~OZOUu«0200 onr .m¥.z.o+z.<*ao.z.4.mu~km+auzouu«Uzou .0¥.2.z.<*.n.2.4.zu~bm+~mzouu.mzou .m¥.z.z.«t.5.Z.J.zm_km+udzouu~~n zosoo_u.m.~.xmoz_ zoaaua~.~.xwoz~ aqzmu.4.zosooa.m .4.zosoo~.mn.4.3oa..m .4.zoa~.mun<3m 2 ._us omm on cam .oom .oom .z.u_ a<3mu24.zosoo_.< .4.xosoo~._n won or maoa mozaa_u.zosoo_.po>.a_ wozakzoo mozakzoo .x.5.~a~. ma 2.2ux co" oo oo .moa .oo “as.w.po>_aa. ma 2.2ua mo" 00 o.oax~a mom roa_a~ 2.2u3 cm 00 c.2n2awsmo zo_»_a .am.mo.xmoza ..a.nosm ..no.mo.<...no.ho>~na zoamzwzao .zamkwo.z.m.z.<.>z.»_nnzoz monoma .acho>aa\.4.zosoo..mn.4.:osoo_.m z.an4 orm 00 com .oom .omn .z.u_ aa.»o>_a\.4.zonoo_.«u.4.zosoo_.< z._u4 omn oo o.~u.zoaooa.zosoo_.< kzmzmsm ko>~a >m 30a ho>aa mo~>_o hacko>~a¢zawsm0usamhmo 0404. OmF 04F O~F mOF OOF OFO 000 Omo 040 000 ONO Ono 000 0mm 00m 004 mmv Omv 004 0N4 004 000 000 OFm OOM mmm 0mm 040 000 ONm 0.0 0.0 0.0 oo~ 0.0 0.0 010 OO+mOOOoO ON 0— 0 7. C~O~C~C~C~C"C~C~0mO~O~O~C—C~C~CHC~O~O~040~O~OHC~O~O~O~O~C~O~O~O~O~0~O~O~O~O~anu 1. oao.c_o~c~o~caoao—oaoaoacacac.oaoaoacaoao_o~o~o.o~o_oao~o_o~o~o_o~o.o~oaoao—o~o~ . oo+wooo.o moum mb.m| oo+woom.~ "and moo. moo. moo. moo. moo. moo. moo. moo. moo. moo. moo. moo. ooomm.o mmmom.o hoom~.o oomm..o mmnmo.o omaco.o ooooo.o ommam.m mobo_.n mnmmm.m mommo.m rucmm.¢ mmom>.¢ oomma.m oo+moom.o+mo+moo.om+ 0o." o0.” oo.~ co." maaumaommemmmaa .~I.~+.~+.~+.an._+.—+.a+.~u ombmmmmaommamb mmoccm.o+.a+.an.~n.o+.~+.~+.~u.o+ mmmcmcmcmmmmmmmm mmmmm"mamcmmmmmmmomemmmmmmmammmmwm.cmcmcammmmmammmmma"Na.cmmamammcaenmmmamaamwam ooom.o_.zoa. 171 E. 4. Binary Deck Structure (a) Job card (b) Binary object program (c) Run card (d) 6 data cards used for all problems (e) 15 data cards used to describe problem as follows: Card (all listing are from bottom to tOp) Formate (l) 4 stiffness constants 4F5. 2 Primary membrane Secondary membrane Primary bending Secondary bending (2) Young's modulus, Poisson's ratio E10. 2, E10.3 (3) Radial component of node points 7F10 (4) Vertical component of node points 7FlO (5) Thickness of triangles 12F5 (6) Boundary conditions are 411 described by 4 integers A, B, C, D. A, B are for the top edge and B, C are for the bottom edge. A, C = 0 or 1 if radial displacement is 0 or unconstrained B, D = O, or 1 if rotation tangent to the edge is O or unconstrained (7) Relaxation constant E10. 3 172 (8) Axial displacement increment (9) Internal pressure (10) and (11) Axial displacement increment function (12) Number of periods (13) Number of axial displacement increments (14) Imperfection constant (15) Imperfection function E10.3 E10.3 80F2.l I2 12 E10.3 7FlO