A STIFFNESS MODEL FOR PLANE STRESS ANALYSIS by PAUL D. HIBBERT B.A.Sc. (Civil Eng.) The University of British Columbia, 1958 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF ' • MASTER OF APPLIED SCIENCE in the Department of CIVIL ENGINEERING We accept this thesis as conforming to the required standard The University of British Columbia April 1965 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make i t freely available for reference and study. I further agree that per-mission for extensive copying of this thesis for scholarly pur-poses may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of CIVIL ENGIKEERIKG The University of British Columbia, Vancouver 8, Canada Date A p r i l . 1965 - i i -Abstract The s t i f f n e s s properties of a f i n i t e sized square model element of an isotro p i c plate under plane stress are presented. The model i s a piece of the material i t s e l f and not a configuration of bars that replace the material. I t i s shown that the use of an extra kinematic condition instead of an a r b i t r a r y displacement function usually produces better r e s u l t s . An example shows that t h i s extra condition greatly increases the accuracy of the model and so reduces the number of degrees of freedom required to approximate the continuous system. A system of bounding the solution i s shown by assigning values to free parameters. Thus the maxi-mum error i s known and can be reduced without going to a f i n e g r i d . The derivation of the is o t r o p i c square i s then expanded to obtain the s t i f f n e s s matrix of an elemental orthotropic rectangle. The l i m i t s of s t a b i l i t y f o r both the is o t r o p i c and orthotropic models are investigated. - i i i -ACKNOWLEDGEMENTS The author wishes to express his gratitude, to Dr. R.F. Hooley for his invaluable guidance during the research and in the preparation of this thesis. Also, the author is grateful to the National Research Council of Canada for financial support in the form of an assistantship, and to the staff of the Computing Centre at the University of British Columbia for their help. April, 1965 Vancouver, Br i t i s h Columbia - i v -TABLE OF CONTENTS Abstract A cknowle dgement s Contents Notation Page I. Introduction 1 I I . Isotropic Square 3 2.1 Derivation 3 2.2 C r i t i c a l Value of a for Isotropy 7 2.3 Bounding Technique 8 2.4 Example 10 2.5 The Determination of Stresses from Known Deflections 11 I I I . Orthotropic Rectangle 14 3.1 Derivation 14 3.2 C r i t i c a l Values of a and P for Orthotropy 16 IV. Conclusions 17 Table I 20 14 Figures 21 Bibliography 35 Appendix A. l Description of the Computer Program, P0RT I 36 A.2 Si m p l i f i e d Flow Diagram 40 A.3 Fortran IV L i s t i n g 41 V Notation In some cases the computer program symbols vary from those defined here. The major differences are pointed out i n the description of the program, which i s to be found i n Section V. The word "entry" i s used throughout the text to refer to a matrix element i n order to avoid confusion with a model element. b,c,d,e,f,g,h i,k,l,m,n,o,p * off-diagonal entries of S E ' » i s o t r o p i c modulus of e l a s t i c i t y E , E - modulus of e l a s t i c i t y i n subscripted d i r e c t i o n x y F «• force vector F i (i=l,2,3,...8) - nodal force G <= shear modulus I » moment of i n e r t i a L « horizontal dimension of model M = bending moment R » aspect r a t i o ( v e r t i c a l dimension to horizontal) S » s t i f f n e s s matrix of element T. • transformation matrix t thickness a, 0 diagonal entries of S shear s t r a i n /x Gt/4 R RGt/4 A SEE deflection vector t r i a l d eflection vector A' A transposed - v i -A. (i-1,2,3,. A6, A3 ( X Ax ^xy* ^yx o~" , o x' y ,8) = nodal deflection - angle changes - s t r a i n vector ^x <y « s t r a i n i n subscripted d i r e c t i o n - Et / 2 ( l - n 2 ) - B E x t A ( l W y X ) • Poisson*s r a t i o f o r isotropy • Poisson*s r a t i o f o r s t r a i n i n f i r s t subscripted d i r e c t i o n , stress i n second stress vector cr-stress i n subscripted d i r e c t i o n shear stress Rt E /12 x' t E y/12R A Stiffness Model For Plane Stress Analysis I . Introduction Various models have been proposed for the approximate solution of the plane stress problem, notably those of Hrennikoff ( l ) a , McCormick (2) , Melosh (3), Clough ( 4 ) , and Turner et a l . ( 5 ) . I n a l l cases the plate i s divided into a g r i d of a f i n i t e number of elements. In the analysis i t i s considered that the displacement of only the nodes i s s u f f i c i e n t to describe the deformation anywhere within the element. The models of Hrennikoff and McCormick consist of various configurations of bars proportioned to d u p l i -cate the e l a s t i c properties of the material they replace, whereas Melosh, Clough, Turner et al.consider the element to be a piece of the material i t s e l f , and with the aid of energy c r i t e r i a they obtain more accurate s t i f f -ness matrices than were available from the l a t t i c e concept. The s t i f f n e s s matrix presented i n t h i s thesis follows the procedure outlined by Turner et a l . of applying unit deformations successively to the eight degrees-of-freedom (two t r a n s l a t i o n a l deflections are assumed at each node of a rectangle). As Turner pointed out, the same resu l t can be obtained from s t r a i n energy considerations. In either case, the general d i s t o r t i o n of the rectangular element i s considered to be a superposition of the d i s t o r -tions due to tension, shear, and bending. In ( 4 ) Clough presented i n d e t a i l the s t r a i n energy derivation, based on superposition, of a s t i f f n e s s matrix of an isotropic rectangle and an i s o t r o p i c t r i a n g l e ; the s t i f f n e s s matrix of the t r i a n g l e however was based on superposition of only the cases of tension and shear. Later (6) he obtained a s t i f f n e s s matrix of an ortho-tropic, triangle but did not present r e s u l t s . Numbers i n brackets refer to the bibliography. Instead of obtaining the general d i s t o r t i o n of the element by super-position, Melosh generated a displacement function. The function gave the displacement of points within the element as ordinates to a hyperbolic paraboloid. He presented a s t i f f n e s s matrix based on the above displacement function and found that i t s a t i s f i e d macroscopic equilibrium and the require-ments of monotonic convergence ( i . e . as the c e l l size i s reduced the solution i s never farther from the exact solution). The use of Melosh's matrix, as w i l l be shown, leads to results which are farther from the exact solution than Clough's, although i t i s nearer to the exact than Hrennikoffs or McCormick's. A l l of the aforementioned matrices y i e l d solutions which converge on the actual solution from below, that i s , the deflections are too small. Therefore i f a s t i f f n e s s matrix i s used which y i e l d s deflections which are too large, the solutions w i l l converge on the actual solution from above. I f both of the curves converge monotonically then the true solution must l i e between them and thus the maximum error i s known for the p a r t i c u l a r g r i d size used. I t is shown i n t h i s thesis that i f the bending c r i t e r i o n is deleted and i f instead some of the diagonal entries of the s t i f f n e s s matrix are considered as free parameters, matrices y i e l d i n g any magnitude of deflec-t i o n can be obtained; hence the true solution can be bracketed, or bounded from above and below. - 3 -I I . Isotropic Square 2.1 Derivation Consider a homogeneous plate stressed i n i t s plane as shown i n F i g . 1(a). Suppose that an elemental square Q i s removed and that i t has the same boundary stresses acting as were present i n the loaded plate, then the general stress d i s t r i b u t i o n can be represented by F i g . 1(b). Furthermore, i f the element i s small, the stress d i s t r i b u t i o n can be approximated as l i n e a r as shown i n F i g . 1(c). Moreover, the stress d i s t r i b u t i o n of F i g . 1(c) i s replaced by the s t a t i c a l l y equivalent set of eight forces shown i n F i g . 1(d). I t should be noted that when the forces have been gathered at the nodes i t i s impossible to determine i f the forces arise from shear or normal stress, or from a combination of the two. S i m i l a r l y the actual deflections of the square Q as depicted i n F i g . 1(e) can be approximated as l i n e a r , as shown i n F i g . 1 ( f ) . The deflection of only the nodes, as shown i n F i g . 1(g), i s s u f f i c i e n t to describe t h i s l i n e a r d i s t o r t i o n . An 8 x 8 s t i f f n e s s matrix S w i l l govern the behaviour of t h i s element such that • SA - F (1) Where F and A are eight-component vectors which represent the forces and deflections of Fi g s . 1(d) and 1(g). The f i r s t column of S represents the forces existing at the nodes when A^ = 1 and a l l other displacements are •zero. F i g . 2(a) shows the element distorted under these forces. S i m i l a r l y the second column of S contains the forces shown i n F i g . 2(b) when A 2 = 1 and a l l the rest are zero. These l a t t e r forces are obtained by rotating the '. element of F i g . 2(a) about i t s Z - axis. The remaining s i x columns can be derived from the f i r s t two by ro t a t i n g the elements about an axis normal to the paper. The entries i n S are shown i n F i g . 3. I t should be noted that • - 4 -columns 2 to 6 are merely permutations of column 1 and that there are only eight unknowns. Various relationships between the entries of S w i l l now be established. Because S must be symmetrical about i t s main diagonal c - g (2) Three more relationships are obtained from equilibrium of the forces of F i g . 2(a) • 2V = 0 b + c = f + g ( 3 ) EH = 0 a + h = d + e ( 4 ) XM = 0 cc - d ' - b - g ( 5 ) In order that the model represent the r e a l material, i t must exhibit the same e l a s t i c properties. Consider f i r s t a model element of length L subjected to a uniform t e n s i l e stress o— . F i g . 4 shows t h i s stress gathered y to the nodes as concentrated forces, and the nodal deflections which are derived from Hooke's Law. The following equations are obtained when the rec i p r o c a l theorem i s written between F i g . 2(a) and F i g . 4 and between F i g . 2(b) and F i g . 4 o- L o- L -(a + h) u - f - + (b + c) —Z- - 0 ( 6 ) o- L o-L o - L t -(b + c) u. + (a + h) -f- = — | — (7) 'Where E and > are the e l a s t i c modulus and Poisson*s r a t i o respectively and t i s the plate thickness at element Q. Consider next the model element under the action of a uniformly d i s t r i -buted shear stress, T , as shown i n F i g . 5 with the stresses gathered to form forces at the nodes. In t h i s case A. = - A. = T L / G , where G i s the shear modulus. Equation (8) i s obtained when the reciprocal theorem i s written between the conditions of F i g . 2(a) and F i g . 5 f - <*) « 2Jj± (8) Although the model element i s i n equilibrium and has the same e l a s t i c properties as an element of the r e a l material, there are, only seven equations l i n k i n g the.eight unknowns of S. However, a solution can be obtained i n terms of one unknown, say a. I f t h i s i s done, then (1 + u.) (9) X 4 X 4 X 4 X 4 ( 3 u ' - l ) (10) (2|A - 2) + a (11) (6 - 2u.) - a (12) f - b / (13) g - c (14) h = •X - a (15) where X = — ^ 2 ( 1 ^ 2 ) Hrennikoff's a r b i t r a r y but ingenious choice of a s p e c i f i c configuration of bars f o r his model provided the remaining l i n k between the eight unknoxvns and t h i s y i e l d s a = X ( 3 0 - 18u,)/24. S i m i l a r l y , McCormick's choice of bars y i e l d s a = A ( 3 3 - 27u.)/24. Neither of these two models, however, can accurately represent the actual deformed state of the element because no combination of a square distorted to a rectangle (no angle changes at the corners) and a square distorted to a - 6 -parallogram (equal angle changes at the corners) can form the general deflected state shown i n F i g . 6, where A6 £ A(3. In other words, i n order to define a generally deflected element i t i s necessary to define f i v e independent displacements. These may be: (a) four length changes and one nodal angle change, or (b) three length changes and two nodal angle changes, or (c) two length changes and three nodal angle changes. Extension i n the x direc t i o n y i e l d s one independent length change. Extension i n the y d i r e c t i o n gives another, but because of the symmetry due to isotropy i t i s i d e n t i c a l to that of x extension. Shear d i s t o r t i o n produces one inde-pendent angle change. Two displacements remain to be found. These kinematic deficiencies may be removed by means of a t h i r d case of the element distorted i n bending. Bending stresses applied separately to the v e r t i c a l and h o r i -zontal sides produce two independent changes, but due to isotropy they are equal. To derive the bending case refer to F i g . 7, which shows that under an applied moment, M, the angle change between the v e r t i c a l sides i s ML/EI, where I i s the moment of i n e r t i a . Then _^ = A^ = -A_ = - A g = ML /4EI and the remaining deflection are zero. As w e l l , F^ * F^ = - F^ = - Fg => M/L and the others are zero. Equation (16) i s obtained when the reciprocal theorem i s written between the conditions of F i g . 2(a) and F i g . 7. ' (a + d - e - h) g - - I (16) Note that t h i s extra equation i s necessary and s u f f i c i e n t to define the s t i f f -ness matrix of the model element. The solution of the f u l l set of eight equations i s set forth i n the second column of Table I and the results are compared with those of Hrennikoff, McCormick, and Melosh. Later examples 3how that the value of a derived from the bending c r i t e r i o n produces s i g -n i f i c a n t l y better results than the other three. 2.2 C r i t i c a l Value of a for Isotropy Another way to remove the kinematic deficiencies i s to consider o as a free parameter rather than f i x i t s value by a bending c r i t e r i o n . A lower bound on a can be found by demanding that energy be stored rather than released as the model i s strained. This requires that f o r any t r i a l vector A, the following holds true A ' S A ^ 0 (17) In other words, S must be positive semi-definite. This implies that a l l eight upper l e f t i x i subdeterminants must be greater than or equal to zero. This s t i p u l a t i o n applied to the upper l e f t l x l subdeterminant y i e l d s a ^ 0 and to the upper l e f t 2 x 2 subdeterminant yields a > b. To avoid the evaluation of higher order subdeterminants, i t i s convenient to change, temporarily, the vector numbering convention established i n F i g . 1(d) and 1(g). A different a r b i t r a r y sequence of numbering the force and deflec-t i o n vectors i n no way changes the result of the derivation, provided that consistency i s observed. Except for the purpose of s t a b i l i t y investigations, the numbering convention of F i g . 1(d) and 1(g) i s observed throughout the theory. Consider vectors 2 and 3 interchanged. I f the derivation proceeds on t h i s basis, S remains the same except that rows and columns 2 and 3 are interchanged. This manipulation, however, has no effect on the algebraic sign of any of the subdeterminants. The new matrix i s governed by the same s t a b i l i t y c r i t e r i a as the old; hence the upper l e f t 2 x 2 subdeterminant must be ^ 0, therefore a ^ c. Thus be interchanging the remaining vector - 8 -numbers (4 to 8) with vector 2 and evaluating the upper l e f t 2 x 2 subdeter-minant i n each case, a l l eight s t a b i l i t y l i m i t s are obtained, two of which are duplicated. The solution i s : a>Q (18) a>b a > X ( l + u)/4 (19) a > c a 5*. X(3ji - l ) / 4 (20) a > d a ^ X(2fi - 2)/4 + a (n < l ) (21) a > e a ^ X(3 - u)/4 (22) a > h a > X /2 (23) For | i « 1 the above l i m i t a t i o n s are a l l contained i n a ^ X/2 and for any smaller positive value of u. a l l l i m i t s are contained i n a ^ e. Therefore the most severe r e s t r i c t i o n i s cc ^ M 3 - n)/4 and a must be the largest entry i n S. I f a = X (3 - u)/4 then the system i s unstable and the deflection approaches i n f i n i t y . I t i s obvious from i t s d e f i n i t i o n that as a approaches i n f i n i t y , the deflection approaches zero. Therefore, some value of a greater than X(3 - u.)/4 must produce the exact deflection. 2.3 Bounding Technique I f the error decreases as the size of the c e l l s decrease, t h i s error can be bounded by solving the problem f o r several sizes of c e l l and two values of a. I n doing so, the s t i f f n e s s matrix must be formed using equations (9) to (15) i n c l u s i v e . A l l that i s necessary i s f o r one value of a to be greater than that required f o r an exact solution and the other to be l e s s . Those solutions f o r a larger than that of the exact w i l l exhibit deflections which increase with a decrease i n c e l l s i z e , and so approach the exact solu-t i o n from below. Solutions f o r a smaller than the exact solution w i l l show deflections which decrease with a decrease i n c e l l s i z e and approach the true solution from above. Since the exact solution i s approached from above and below, the deflection i s bounded and the maximum error i s known. I f both of the two values of a used produce deflections which approach the true from the same side, then one value must be changed so as to approach from the other side. An obvious advantage i n t h i s system i s that there i s no need to keep decreasing the g r i d size i n order to gain f a i t h i n accuracy. There i s nothing i n the above argument to show that the same value of a produces the exact deflection at a l l points. In f a c t , a different value of a i s required to produce the exact solution at each point. Having shown that deflections may be bounded, i t i s easy to show that strains may also be bounded. Consider, f o r example, a system i n which the deflection of node number 1 i s bounded by two values of a and that of node 2 by two d i f f e r e n t values. Suppose that the maximum errors, as found by bounding, i n the deflections of points 1 and 2 are 6^ and 6^ respectively, then the maximum error i n s t r a i n between the two points i s (16^1 + |62I)/L. When expressed as percentage errors, the error i n s t r a i n may thus be greater than either error i n deflection. The s t r a i n , nevertheless, i s bounded. By a s l i g h t l y more complicated argument which considers the deflections of four points instead of two, i t can be shown that stresses may also be bounded. Melosh (3) gives a procedure f o r bounding deflections at a point based on the fact that the s t r a i n energy of the exact solution l i e s between the minimum potential energy and the minimum complementary energy. To obtain the maximum error using t h i s concept i t i s necessary to obtain solutions f o r r e a l and influence loads f o r two models; one model must be based on continuous deflections over the body and the other on continuous stresses. Melosh 1 s technique implies that the solution from the continuous deflection model produces a d e f l e c t i o n which i s too small, whereas that from the continuous stress model produces a deflection which i s too large. Therefore the true - 10 -deflection i s bounded from above and below. In order to reduce the error found by Melosh's method i t i s necessary to reduce the grid s i z e . This method of bounding deflections at a point, although complete i n theory, i s r e s t r i c t e d due to the lack of stress consistent matrices and so does not yet have p r a c t i c a l application. A search of the l i t e r a t u r e shows no reported stress consistent matrices or attempts at bounding by Melosh's technique. 2.4 Example A cantilever beam of rectangular cross-section with a length-to-depth r a t i o of three was run with the computer program. PORT I , which i s described i n the Appendix. F i g . 8 shows the end deflection of t h i s beam vs. the number of c e l l s high f o r various values of the parameter a. The value of a deduced from bending deformation gives much better accuracy than the other two models but there i s s t i l l some error, when compared with the e l a s t i c i t y solution f o r no f i x e d end warping. This i s due primarily to the r e l a t i v e l y crude approximation of the required parabolic shear d i s t r i b u t i o n between nodes. The two curves f o r McCormick are for the c e l l s pinned together at the nodes and f o r the c e l l s f i x e d together at the nodes. This model was derived on the assumption that no moment exists at the node and i t i s i n t e r e s t i n g to note that when run i n t h i s manner the re s u l t s are s l i g h t l y better. Further, there i s one degree of freedom less per j o i n t when run t h i s way. The curves show that a/X m .861 i s too large f o r f i n d i n g the end deflection since t h i s curve approaches the true from below whereas a/X = .825 i s too small since the approach i s from above. For four c e l l s high the solu-t i o n i s bounded then between 117.05 and 115.54. Cross p l o t t i n g of the curves of F i g . 8 shows that o/X « .85 should produce the exact solution of 116.64 f o r any number of c e l l s high, as shown i n F i g . 9 . With the exception of 1 c e l l high, computer runs v e r i f y t h i s fact to within 0.55$. F i g . 10 shows a plot of the extreme f i b r e bending stress half way along the beam vs. the number of c e l l s high f o r various values of a/\ . The curves are of the same shape and show that the stress i s bounded above and below by different values of a//\ . In th i s case the value of a/X = .861 from the bending solution i s required to produce the exact r e s u l t . This again i s independent of the number of c e l l s high. Also the parameter a required f o r exact stress i s different from that required f o r exact de f l e c t i o n , as would be expected. 2.5 The Determination of Stresses from Known Deflections In order to determine the stresses from known corner deflections, as i n F i g . 11(a) f o r example, i t i s f i r s t necessary to investigate s t r a i n s . From the unit cases of tension, shear, and bending used i n the model derivation, i t was found, over the area of the element, that (. and 6 ' • x y vary l i n e a r l y under bending moments applied to the v e r t i c a l and horizontal sides respectively, that they are constant under tension, and that V i s constant under shear. Accordingly then, the strains at any point within the element may be defined by plane surfaces as shown i n F i g . 11(b), 11(c), and 11(d). I t i s assumed that the element under consideration i s i n the i n t e r i o r of .a structure; the s t r a i n surfaces of three adjacent members are shown i n the diagrams. Note should be taken of the fact that ( and £ are constant * y on l i n e s p a r a l l e l to the x and y axes respectively. In F i g . 11 there are disc o n t i n u i t i e s of the s t r a i n surfaces at the e l e -ment boundaries. I t i s deemed advisable therefore to evaluate strains (and stresses) at the centre of the element i n order to avoid the discontinuity e f f e c t s , and at the same time, give equal weight to the four corner values. - 12 -When the strains have been evaluated, stresses follow from the w e l l -known relationships: o- x - tfx+ n £ y ) E / ( l - j x 2 ) o- (< y + ^ e x ) E / ( i -(24) (25) (26) The stresses at a boundary, however, can be evaluated more d i r e c t l y . The normal stress perpendicular to the boundary and the shear stress are equal to the applied stress. The stress p a r a l l e l to the boundary i s found from o-' X E i + u o— x •* y (27) for a boundary running i n the x d i r e c t i o n , f o r example. Having established the basic concepts, the transformations r e l a t i n g deflections, s t r a i n s , and stresses at the c e l l centre can now be derived. F i r s t l y , i f the relationship between deflections and strains i s £x and T^ i s an 8 x 3 transformation matrix, ^ => T 1 A, where € then by geometry: •y <x - k ( A 1 + A 4 + A 5 + V (28) 2L" ( A 2 + A3 + A6 + A 7 } (29) therefore V " k ( A 1 + A 2 - A 3 - \ + A 5 + A 6 " A 7 " V ( 3 0 ) 2L 0 0 1 1 1 0 1 -1 -1 1 0 1 0 0 1 1 1 0 1 -1 -1 ( 3 D - 13 -Secondly, i f the relationship between strains and deflections i s given o-X by o- = A, where cr-matrix, then by equations (24), (25), (26) E and i s an 8 x 3 transformation 2L(1V E y 2L(l-u.2) - [ ( A x + A 4 + A 5 + Ag) + ( A 2 + A 3 + A 6 • A ?) U.J (32) [ ( A 2 + A 3 + A 6 .+ A ?) + (A 1 + A^ + A 5 + A g) uT] (33) E 4(l+u)L (34) therefore 2L(l«si ) 1 1 u. 1 H l x n (j. 1 x |i 2 2 2 2 2 2 2 2 (35) Thirdly, because £ ** T^ A and o- = T^ A then the relationship between stress and s t r a i n may be written as o- = T^6 , where T^ i s a 3 x 3 trans-formation matrix and i s defined by T 2 « T^ T^. I t i s found that 1 u 0 E 1-u* 0 o i f (36) - 14 -I I I . Ortho-tropic Rectangle 3.1 Derivation The foregoing approach for finding the s t i f f n e s s matrix of a square i s o t r o p i c model can be extended to allow the analysis of a plane orthotropic continuum approximated by rectangles (see F i g . 12). The associated s t i f f -ness matrix, S,; i s shown i n F i g . 13. Since there i s no longer any i s o t r o -pic symmetry, the second column i s not a permutation of the f i r s t , and so, i n general, there are sixteen unknowns i n S. Note should be taken, though, of the fact that columns 4, 5, and 8 are permutations of column 1 and that columns 3, 6, and 7 are permutations of column 2. Four equations are obtained from matrix symmetry, s i x from s t a t i c s , and four from the e l a s t i c laws of tension and shear. The solution of the fourteen equations mentioned above, i n terms of the two diagonal e n t r i e s , a and f$, i s : b - ^ L ^ + R Y x (37) c - Y " X - R V x (38) d - a - 2 Y x (39) e .•- 2 X x + 2 V x - a . (40) (41) g - c (42) h - 2 X x - a (43) i = b (44) k = 2 A - p (45) 1 - c (46) - 15 -m' - b . (47) n - 2 X + 2 V - p- (48) y y ' V P - 2 V y (49) p - c (50) R E t E t where >v * -rrz r X T_7T~^ r *x 4R * y 4 The f i r s t subscript of n refers to the d i r e c t i o n of s t r a i n and the second to the d i r e c t i o n of stress. R i s the aspect r a t i o and E i s subscripted to show the appropriate d i r e c t i o n . As i n the i s o t r o p i c case, two kinematic deficiencies remain which may be removed by applying bending stresses separately to the v e r t i c a l and horizontal sides. I f t h i s i s done then the complete solution i s : b - ^ p £ + R Y x (52) c - ^ - x - R V x (53) d " \ + *x~ V x (54) e - V - ' ^ x * Y x > -<55):. f - b (56) g - c (57) h - X - r - Y (58) X X X 1 • - b (59) p •-" X y + * y + V y ( 6 0 ) where - 16 -k = A - r - Y (61) y y y 1 - c (62) m" - b (63) n = A y - V Y y . ( 6 4 ) ° ~ X T * * 7 - * 7 ' ( 6 5 ) p = c (66) Rt E t E x 12 7 12 R 3.2 C r i t i c a l Values of a and 3 f o r Orthotropy I f the bending c r i t e r i o n i s omitted, and. a and 3 considered as free parameters, as was done i n the is o t r o p i c case, then c r i t i c a l values of them can be obtained. The use of the same method of renumbering, previously des-cribed, to avoid the evaluation of subdeterminants of higher order submatrices than 2 x 2 to determine positive seml-definiteness, r e s u l t s i n the following threeinequalities to define the. lower l i m i t s of a and p: a p > £ s £ ^ R V x ) 2 (67) a A x + V x . (68) 3 > X y + V y (69) Also i t i s found that G> 0. A plot of the c r i t i c a l values f o r a Douglas F i r element of unit aspect r a t i o , and thickness of one inch i s shown i n F i g . 14. Values of a l e s s than 545.3 kips/inch and/or 3 less than 64.3 kips/ inch produce an unstable element, whereas i n f i n i t e values of a and 3 cause zero deflection of the element under load. I t should be noted that the l o c a -t i o n of a and 3 derived from bending i s i n the stable region. IV* Conclusions The author concludes f o r t h i s work that the most accurate s t i f f n e s s matrix derived so f a r i s that o r i g i n a l l y derived by Turner et a l . by super-position of the three basic cases of tension, shear, and bending. Hrennikoffs and McCormick's arbitrary choice of bars w i l l i n general v i o l a t e some kinematic conditions, which leads to a l e s s accurate model. The accuracy of Melosh's model, which i s based on the a r b i t r a r y selection of a displacement function that s a t i s f i e s the boundary conditions on displace-ment i s l i m i t e d to the accuracy of the displacement function i n representing the true deformation. Melosh chose a very simple function but the i d e a l one i s probably much more complicated. I t appears that i n order to extend Melosh's work i t would be necessary to generate displacement functions of a progressively more complicated form, derive the s t i f f n e s s matrix for each one, and compare computer results against known solutions u n t i l a better model was obtained. Nevertheless, an analyst would be confronted with an unknown error, or i f Melosh's bounding technique was used, the error could not be reduced except by a reduction of the c e l l s i z e . _ With the computer program f o r the bounding technique developed herein, the free parameter, a, f o r isotropy, can be chosen so that the models of Turner, Hrennikoff, McCormick, or Melosh are obtained; or two values of a can be chosen so as to bound the true solution from above and below. The pro-cedure f o r bounding the deflection at a node i s as follows: 1. Select a value of a which i s known to produce a deflection which i s too small. For a f i r s t t r i a l i t i s recommended that a be chosen to produce the model of Turner et a l . because Turner's model i s known to be too s t i f f yet gives better accuracy than any other. 2. Obtain computer solutions f o r two grid sizes. Because the value of a chosen i n step 1 produces deflections which are too small, the l i n e joining the two solutions on a plot of nodal deflection vs. g r i d size w i l l have positive slope (e.g. the Clough and Hibbert curve of F i g . 8). 3. Select another value of a so that a deflection which i s too large w i l l be produced. 4. Obtain computer solutions f o r the same two g r i d sizes a3 i n step 2. In t h i s case the l i n e j o i n i n g the two solutions must have negative slope to be v a l i d . I f the slope i s positive then the value of a selected i n step 3 was too large and the deflections are smaller than the exact. In t h i s case repeat step 4 with a smaller value of a u n t i l a negative slope i s obtained. Two curves have then been obtained, one of which has positive slope, the other negative (e.g. F i g . 8, the Clough and Hibbert curve and the a/X *= 0.82$ curve). 5. I f the selected values of a are known to produce solutions which con-verge monotonically, then the actual nodal deflection i s known to l i e between the two curves ( i . e . 114.4 < -p- < 120.6 f o r two c e l l s high, i n F i g . 8). 6. I f the maximum error so obtained i s unacceptable, show the results on a plot of nodal deflection vs. a with the grid size as a parameter, as i n F i g . 9. Select a more precise value of a on the basis of the i n t e r -section of the two c e l l - s i z e curves. Obtain solutions f o r the new a and obtain a smaller maximum error. In t h i s manner the true solution i s always contained between the two curves of smallest positive and smallest negative slope and the error can be reduced further. I f the s t r a i n along a l i n e joining two nodes i s required, then the values of a must.be selected to bound the s t r a i n ; or i f the stress at the centroid of an element i s required, then the a's must be chosen to bound the stress. - 1 9 -The bounding technique described above w i l l reduce the number of degrees of freedom required to approximate a continuous system and so increase the complexity of problem which present day computers can solve. - 20 -TABLE I . ENTRIES OF S a Entry Hibbert, . Hrennikoff McCormick Melosh Turner, Clough (1) (2) (3) (4) (5) a 22 - 6u - 4u 2 30 - 18u 33 - 27n 24 - 8u b 6 + 6\x 6 + 6u. 6 + 6\i c - 6 + 18^ - 6 + 18u - 6 + 18u - 6 + 18u d 10 + 6u - V 2 18 - 6u 21 - 15u 12 + 4|i e 14 - 6|i + 4n 6 + 6|i 3 + 15H- 12 - > f 6 + 6^ 6 + 6\i 6 + 6u. 6 + 6(j. g - 6 + 18u - 6 +• 18JA - 6 + 18u - 6 + 18|x h 2 + 6^ + 4u 2 - 6 + 18n - 9 + 27^ 8u a M u l t i p l y a l l entries by A/24 - 21 -- 22 -C e <-/ / 7 1 v f f <-( a ) Fig. 2 FORCES FOR COLUMNS I AND 2 OF 23 -S = a b g d e f c h b a h c f e d g c h a. b g d e f d g b a h c f e e f c h a b g d f e d g b a h c g d e f c h" a b h c f e d g b a F ig. 3 ST I FFNESS MATRIX OF AN ISOTROPIC SQUARE -2k -<ryL t A c r y L t F ' o-yLt " 7 F cry L ~ E ~ V i 1 i L 1 1 • • 1 y L c r y L t Fig. 4 UNIFORM T E N S I O N -25 -r y 2 T L H i Ay T L t T L t 2 I k TLt ? / / / T L t G K 1^ " T L t 2 T L t Fig. 5 UNIFORM .SHEAR - 26 -Fig. 6 GENERAL DISTORTION - 2 7 -Fig. 7 BENDING 28 -3 x depth 0 - 29 -- 3 0 -NUMBER OF CELLS HIGH Fig. 10 BOUNDING THE E X T R E M E FIBRE STRESS . - .31 -(d) T H E VARIATION OF STRAINS OVER T H E E L E M E N T - 3 2 -R L A y L Fig.,12 ORTHOTROPIC R E C T A N G L E - 33 -S = a i 1 d e m P h b /3 k c f n 0 g c k ,8 b g 0 n f d 1 i a h P m e e m p h a i 1 d f n o g b k e g o n f c k b h p m e d I i a Fig. 13 S T I F F N E S S MATRIX OF A N ORTHOTROPIC R E C T A N G L E - 34 Fig. 14 CR IT ICAL V A L U E S OF a AND /3 FOR D. FIR - 35 -Bibliography 1. Hrennikoff, A., Solution of Problems of E l a s t i c i t y by the Framework Method, ASME, December, 1941, p. A169. 2. McCormick, G.W., Plane Stress Analysis, Journal of the Structural D i v i s i o n , ASCE, V o l . 89, No. ST4, Proc. Paper 3581, August, 1963, pp. 37-54. 3. Melosh, R.J., Development of the S t i f f n e s s Method to Bound E l a s t i c Behaviour of Structures, Ph.D. Thesis, University of Washington, at Seattle, Wash., August, 1962. 4. Clough, R.W., The F i n i t e Element Method i n Plane Stress Analysis, 2nd Conference on Electronic Computation, ASCE, September, i960. 5. Turner, M.J., Clough, R.W., Martin, H.C., Topp, L.J., St i f f n e s s and Deflection Analysis of Complex Structures, Journal of the Aeronautical Sciences, V o l . 23, No. 9, September, 1956. 6. Clough, R.W., Stress Analysis of a Gravity Dam by the F i n i t e Element Method, B u l l e t i n Rilero, No. 19, June, 1963. : : '. 7. Chpleski, Contributions to the Solution of Systems of Linear Equations and the Determination of Eigen Values, Edited by Olga Taussky, National Bureau of Standards, U.S. P r i n t i n g O f f i c e , Applied Mathematics Series No. 39, September, 1954, pp. 31. - 36 -Appendix A A . l Description of the Computer Program, PORT I Dr. R.F. Hooley's computer program, ALP I I , was rewritten from Fortran I I language to Fortran IV and modified to solve s p e c i f i c a l l y the problem of plane e l a s t i c i t y . The program and data are stored e n t i r e l y within the core of an IBM 7040; therefore tapes are not required and the operating time i s very short. B r i e f l y , the description of the sequence of operations and the presenta-t i o n of data required i s as follows: 1. A card i s read which bears the structure data, which consists of: (a) an ar b i t r a r y structure number, NRS, (b) the number of j o i n t s , NRJ, i n the system, (c) the number of rectangular members, NRM, (d) the number of e l a s t i c r e s t r a i n t s , NRE, imposed, and (e) the number of sets of loads, NLC, to be analysed. 2. Next, a card i s read which contains the e l a s t i c properties of the material ( i . e . E , u , E , G) and then u i s calculated by means of the x *xy y yx reciprocal theorem. 3. The next data consist of a series of cards (one card i s required f o r each j o i n t ) ; each bears: (a) the joint number, JN, (b) either a figure 1 i f the joint i s permitted to move horizontally or a figure 0 i f not, (c) either a figure 1 i f the jo i n t i s permitted to move v e r t i c a l l y or a figure 0 i f not. I t i s desirable to obtain a structure s t i f f n e s s matrix with as narrow a diagonal band as possible i n order to minimize the computing time and - 37 to allow the analysis of large structures. • To r e a l i z e a narrow band width, careful consideration must be given to the a r b i t r a r y j o i n t -numbering scheme. A displacement number, ND, i s assigned to the degrees of freedom after each card i s read and the t o t a l number of degrees of freedom, NU, to date i s computed. When a l l the j o i n t data cards have been read, NU i s printed. 4. Next, the free parameters d and 0 (A and AJ i n the program) f o r unit plate thickness are read from a single card. 5. A sequence of cards which contain the member data i s read now. Each card bears: (a) the member (element) number, MN, (b) the upper r i g h t j o i n t number, JUR, (c) the upper l e f t j o i n t number, JUL, (d) the lower l e f t j o i n t number, JLL, (e) the lower right j o i n t number, JLR, (f ) the aspect r a t i o (R i n the theory; ASP i n the program), (g) the thickness of the plate, T, and (h) the horizontal length of the element (L i n the theory; HORL i n the program). After each card i s read, the s t i f f n e s s matrix of each element i s com-puted and stored, since, i n general, there may be a v a r i a t i o n i n the aspect r a t i o or i n the thickness of elements. Also,, position numbers, NP, are assigned. The band width i s calculated f o r each member i n turn but only the maximum value for the preceding members i s retained. After a l l the cards have been read, the computer prints the band width. 6. . E l a s t i c r e s t r a i n t data i s now read. One card i s required f o r each r e s t r a i n t and i t contains: (a) the j o i n t number, N22, where an e l a s t i c r estraint e x i s t s , (b) the spring type, N23, (a figure 1 f o r a horizontal spring; a figure 2 for a v e r t i c a l spring), and (c) the spring constant, ELAST. 7. At t h i s stage the computer f i l l s i n the structure s t i f f n e s s matrix (S i n the program) member by member with the a i d of code numbers, NPS. S i s stored as a singly subscripted variable. 8. The Choleski (3) or Square Root method of solution i s used, whereby the symmetrical band matrix i s decomposed into upper and lower triangles such that S = ZZ'. Z i s the lower t r i a n g l e and Z« the upper (Z» i s the transpose of Z). One advantage of the Choleski system i s that i t i s necessary to store only one t r i a n g l e ( s l i g h t l y more than half) of S. As w e l l , the zeros outside the band are not stored and i t i s the most precise of a l l direct systems. The decomposition i s accomplished at t h i s stage. 9. Next, the computer examines Z for s t a b i l i t y by using the c r i t e r i o n that states that a l l of the diagonal entries f o r a stable structure must be p o s i t i v e . Zero or imaginary entries cause the computer to print the l o c a -t i o n of the f i r s t c r i t i c a l entry encountered and to reject the structure. 10. Because bounding technique requires several runs with a v a r i a t i o n i n the free parameters, i t may happen, when elements have different aspect r a t i o s , that some of the elements are unstable, even though the structure i s not. I f so, i t i s i n t e r e s t i n g to know which of the values of a, p, or af3 . are too small. The computer prints the number of the unstable members and the cause of i n s t a b i l i t y but does not reject the structure. 11. The computer now reads the number of loaded j o i n t s , NJL, and the load data which consists of NJL cards; each bears: (a) the number of the loaded j o i n t , JNU, - 39 -(b) the force i n the x dir e c t i o n , F F ( l ) , and (c) the force i n the y di r e c t i o n FF(2). 12. Back substitution one i s now accomplished, i n which the computer solves the matrix equation ZV = F f o r V, where . F i s the load vector. 13. A i s now found from the equation V » Z'A, where A i s the deflection vector. This operation i s called back substitution two. 14. Next, back substitution three i s performed. Joint deflections (available at the end of back substitution two) are printed; the strains on a l l four sides, the eight nodal forces, and the three centroid stresses are calcu-lated and printed f o r each member. 15. I f more than one set of loads was specified i n step 1, the computer returns to step 11 to process the next set. The looping i s repeated u n t i l the l a s t load condition i s analysed. 16. The computer now either returns to step 1 to read data for the next struc-ture or proceeds to another program. There i s no conversion of units within the program; therefore a l l data presented must be.of consistent u n i t s . The printed output w i l l , of course, be of the same units as were input. I t i s important to note that i t i s necessary to multiply rows and columns hi 5, 6, and 7 of each member matrix, as derived i n the theory, by - 1 to create a consistent set of directions f o r the structure. In the derivation i t was convenient to define vectors as positive when they pointed outward from the element but t h i s convention must be abandoned when a structure matrix i s to be formed. A.2 SIMPLIFIED FLOW DIAGRAM FOR PORT I Start Read Structure and P r i n t Data i : ^ ; ( Read Material Properties) Calculate ^yx Q Print Material Properties^ Ji. Read • Joint and P r i n t Data Assign . Displacement Nos. t Calculate No. of Deg. of ' Freedom.to Date ( Prin t No. of Deg. of FreedonQ ( a l f p r i n t ^ a n d " ^ 4 Read Member and Print Data i Calculate S for Member Assign P o s i t i o n Nos. Calculate Band Width to Date ( Pr i n t Band Width Read E l a s t i c Con-and P r i n t , s t r a i n t Data Compile Structure S t i f f -ness Matrix Calculation Convert S to Z A Calculation Check Structure ' S t a b i l i t y Calculation Check Element S t a b i l i t y -> Read . Load Data and P r i n t J Calculation-Back Subst. I e I I Pri n t Joint Deflections Calculate Member and Print Strains Calculate . Nodal Forces and P r i n t of Member ~5> Calculate Member and Print Stresses - uo -o o C PORT I P.DoHIBBERT C PLANE ORTHOTROPIC RECTANGLES C INPUT AND OUTPUT CONSISTENT UNITS. C HANDLES NEGATIVE DUPLICATE NODE DEFLECTIONS. 999 CONTINUE DIMENSION JN( 13 5 ) » N D ( 1 3 5 » 2 ) » N P ( 2 0 0 , 8 ) ,SM(36)>S(80Q0) >SMM(200»36) DIMENSION E L A S T ( 2 0 ) » N N 2 2 ( 2 0 ) » F K ( 2 7 0 ) » F F ( 2 ) » A S P ( 2 0 0 ) , H O R L ( 2 0 0 ) DIMENSION .6AMC200)»EXT(200)»EXB(200).EYL(200)»EYR(200)»T(2001 DI MFNSI ON SX(200) >5Y(200),SXY(200) »D(270) >F(8) .SA(8»8) »DM(8) DIMENSION NPS(2 00»8) 20 READ l.NRS»NRJ»NRMfNRE»NLC 1 FORMAT (51 10) : , : PRINT 3 » NRS ? NRJ » NRM FORMAT (1H1 ,13HSTRUCTURE 1S ) NO»I3»6H HAS 15»9H J O INTS"» I 5 ,9H MEMBER > - T l -O 3-- 3-- co c+ H* TO P R I N T 5 » N RE » N L C F O R M A T ( 1 X » 3 H A N D I 5 i 2 7 H 1 I T I O N S . ) ELASTIC CONSTRAINTS WITH I 5 i l 8 H LOAD COND h-1 i 1 2 . J i •10 9 READ 6»EX»UXY*EY»G FORMAT UF10.3)' UYX = UXY*EX/EY PRINT 7» E X , U X Y , E Y » U Y X » G FORMAT (1H0,4HEX = F 10 . 3 » 5H UX Y = F'10 » 3 » 4H EY = F 10 . 3 » 5H UYX = F10.3»3H G 1 = F10.3) , : NU = 1 PRINT 8 FORMAT ( l'HO .35HJOINT NUMBER X DISP. Y DISP.) 799 DO 11 I=1»NRJ READ 799, JN( I )»ND( I ,1)*ND( I»2) FORMAT ( 3 I 10 ) ' PRINT 1 0 » J N ( I ) » N D ( I » 1 ) » N D ( I » 2 ) TO FORMAT ( I X » 3 I 1 0 ) _• DO 11 J = l i 2 _ NUM = ND( I» J) NUMA= IABS(NUM) I F ( N U M A - 1 ) 1 1 , 1 2 , 1 3 1 2 N D ( I , J ) = N U N U = N U + 1 ; ' GO TO 1 1 1 3 I F ( N U M A - I ) 1 4 » 1 5 , 1 5 1 5 P R I NT 1 6 » J N J J J - _ _ _ 1 6 F O R M A T ( 1 X , 2 3 H D I S P L A C E M E N T F O R J O I N T I 4 , 1 9 H I S I N WRONG O R D E R ) S T O P 1 4 I F ( N U M ) 1 7 T 1 5 t 1 8 . •• 1 7 N D ( I , J ) = - N D ( N U M A » J ) GO TO 1 1 1 8 N D ( I , J ) = ND ( NUMA > J ) " l l C O N T I N U E N U = N U - 1 N 3 3 = N U_- 2 7 0 . • _ ' I F ( N 3 3 ) 2 3 , 2 3 , 2 1 2 1 P R I N T 2 2 , N 3 3 2 2 F O RM AT ( 1 X , 3 Q H T 0 0 M A N Y D E G R E E S O F F R E E D O M BY I 4 ) . ' S T O P 2 3 P R I N T 2 4 , N U 2 4 F O R M A T ( IH0__1 Q H T H E R E A R E I 5 , 2 Q H D E G R E E S O F F R E E D O M ) R E A D 2 5 » A , A J P R I N T 2 6 , A , A J 2 5 F O R M A T ( 2 F 1 0 . 3 ) hit T 2 6 F O R M A T ( 1 H 0 , 2 H A = F 1 0 . 3 , 8 H A J = F 1 0 . 3 ) . 1 • P R I N T 2 7 1 12 2 7 F O R M A T ( 1 H 1 ,61H'MEMBER UR U L L L L R A S P E C T T H I C K N E S J.L I S L E N G T H ) 2 " 1 0 P R I N T 2 8 9 __ 28__ F O R M A T ( 1 H + , 4 2 H N U M B E R J O I N T J O I N T J O I N T J O I N T R A T 1 0 ) 3 8 7 N B = 0 DO 4 7 K = 1 , N R M 4 G DO 3 3 1 = 1 , 3 6 5 3 3 S M ( I ) = 0 . 0 5 . 4 R E A D 9 , M N » J U R , J U L , J L L , J L R » A S P ( K ) , T ( K ) » H O R L ( K ) o o PRINT 2 »MN, JUR, JUL, J L L , JLR,ASP(K) »T (K ) »H0RL( K ) 9 FORMAT <5I5,3F10o3) 2 F O R M A T ( 1 X » I 5 » 4 I 7 » 3 F 1 0 . 3 ) PARAMETERS S'MA = A#T(K) SMB = UXY*EX*T(K-)/(4.*(lo-UXY*UYX) ) + G*T(K)/4» I SMC = SMB - G*T ( K ) /2o. SMD = SME =-SMA - G.*T(K')/(2o*ASP(K) ) SMD + ASP(K)'*EX*T(K)/(2«*ll.-UXY*UYX) ) I SMH = SME - G*T(K')/(2o*ASP(K) ) ! SMJ = AJ*T(K) j SMK =-SMJ + FY*T(K)/(2.*ASP(K)*(1.-UXY*UYX) ) SMN = SMK + ASP(K)#G*T(K)/2o SMO = SMJ - ASP(K)*G*T(K)/2o C STIFFNESS MATRIX FOR RECTANGLES 2 9 SM( 1) = SMA i SM ( 2 ) = SMB i SM( 3) = SMC SMC 4) = -SMD SM( 5) = -SME SM( 6 ) = -SMB SM( 7) = -SMC SM( 8 ) = SMH . SM( 9 ) = SMJ I SM(10) = SMK SM(11) = -SMC 1 12 SM(12) = -SMB 11 SM(13 ) = -SMN ' 2 " " 1 0 SM(14 ) = -SMO . • c SMC15) = SMC 3 ! i SM(16 ) = SMJ 7 SM(17 ) = -SMB 4 fi SM(18 ) = -SMC 5 SMC 19 ) = -SMO 5 ' 4 3 SMC 20 ) = -SMN o SM(21) .= SM(22) = SM(23) = SMB SMA SMH SM124) = SM(25) = SM( 26) = SMC SMB -SME SM(27) = SM(28) = SM(29) = SMA. SMB SMC SM(30 ) = SM( 31 ) =• SM(32) = -SMD SMJ SMK SM(33) = SM(34) = SM(35) = -SMC SMJ -SMB " ' i c SM(36) = CALCULATE NP(-K»1) = SMA POSITION NUMBERS FOR RECTANGLES ND(JUR,1) p NP(K,2) = NP(K»3) = NP(K,4) = ND(JUR 9 2 ) ND(JUL,2) ND(JUL,1) NP(K,5) = NP-(K»6J = NP(K»7) = N D ( J L L t l ) ND(JLL,2) ND(JLR 9 2) c NP(K»8) = CALCULATE MAX = 0 ND(JLR 9 1 ) BAND WIDTH 9 NB 50 MIN = 500 DO 39 I = 1 »8 N =IABS(NP(K,I)) 40 42 IF (N)39,39,40 IF (N-MAX)41,41*42 • MAX = N 41 4 3. IF (N-MIN)43 ,39,39 MIN = N o o . •" . . . . . -..v.-1.---. . . - - • J 39 CONTINUE NBl = MAX - MIN IF ( NBl - NB)44,44,45 45 NB = NBl 44 . DO 46 1 = 1,36 -21 -. 46 SMM(K»I) = SM(I) 4 7 CONTINUE NB = NB + 1 MIKF = ?*NB - 1 • PRINT 48,NB,MIKE 48 FORMAT (1H0»3HNB=I 5»23H TOTAL BAND WIDTH = I 5 ) IF (NRE) 51 , 51 » 52 52 PRINT 53 53 FORMAT (1H0,39H JOINT TYPE ELASTIC CONSTRAINT) i DO 54 I =1,NRE READ 55,N22,N23,ELAST(I) vn PRINT 55,N22,N23,ELAST( I ) 1 55 FORMAT ( I X , I 10,I 10,E10.3 ) 54 NN22( I ) = ND(N22,N23) 51 CONTINUE C FI L L IN STRUCTURE STIFFNESS MATRIX • NV = (NU*NB) + NU + NB - 2 NO = NV - 8000 IF (NO) 56 , 56 , 57 57 PRINT 58, NO 5 8 FORMAT (1X,11H0VERLAP OF I4,20H IN STIFFNESS MATRIX) STOP 56 MS = NU*NB DO 59 I=1,MS 59 S ( I ) = 0.0 J4 = ((NB+1)*NU)+1 J 5 = J 4 + N B - 3 DO 60 I=J4,J5 60 S ( I ) = 0.0 NBl = NB - 1 o - ^ :> ' ••• •." : .7 . ... . • . ^ 7 . • ; . • '7 , DO 471 I = 1 • NRM DO 471 J = l * 8 471 N P S ( I , J ) = N P ( I . J ) ' DO 61 L=1»NRM J 62 64 DO 64 K = 1» 8 NP(L»K) = IABS(NPS(L »K) ) DO 61 J=1» 8 -;) : 65 •IF(NPS(L»J) )65»61»66 AB= -1.0 GO TO 67 *) 66 AB= 1.0 6 7 J l = ( J - l ) * ( 1 6 - J ) / 2 DO 61 I=J»8 76 I F ( N P S ( L » I ) > 6 8 , 6 1 , 6 9 68 B = -AB GO TO 70 69 70 B = AB IF(NP(L»J)-NP(L,I) )71,72,73 1 72 I F ( I - J ) 7 4 . 7 1 , 7 4 74 K = ( N P ( L , I ) - 1 ) * N B 1 + NP(L,J) IJ1 = I + J l S(K) = S(K) + (2.*B*SMM(L,IJ1 ) ) GO TO 61 71 K = (NP(L»J)-1)*NB1 + NP(L,I) w T ' GO TO 75 7 73. K = (NP(L * I)-1)*NB1 + NP(L,J) 12 75 I J l = I + J l 11 S(K) '= S(K) + B*SMM(L, U I ) " 1 0 61 CONTINUE 9 IF (NRE) 121,121,470 ' 8 470 DO 469 I = 1,NRE 7 K = ((NN22(I) -1)*NB) + 1 C 469 SCK) = S(K) + E L A S T ( I ) 5 . 121 CONTINUE 4 3 C SOLVE o o KASE =1 L1=NB+1 K3=NB-1 L2=(NB*NU)-.K3 NBS = K3*K3 NB2= NB-2 L = 1 IF(S(1>)212,212,216 216 S(1)=SQRT ( S ( l ) ) DO 202 I =2 »NB 202 S(I) = S ( I)/S<1 ) DO 203 L=L1,L2»NB K1=L-NBS K11=K1 K2=L-K3 i IF(K1)204»204,205 r— 204 Kl= ( ( L - l ) / N B ) + 1 -o 205 . DO 206 K=K1,K2,K3 i 206 S(LJ= S ( L ) - S ( K ) * S ( K ) I F ( S ( L ) ) 2 1 2 , 2 1 2 , 2 0 7 207 S(L)= SORT ( S ( L ) ) DO 208 N=l,NR2 M1=K11+(K3*N) IF ( M l ) 2 0 9 , 2 0 9 ,210 209 M1 = K1 210 LN=L+N DO 211 M=M1,K2,K3 MN=M+N 211 S(LN)= S ( L N ) - ( S ( M ) * S ( M N ) ) 2 08 S(LN)= S ( L N ) / S ( L ) 203 S(LN+1)=S(LN+1)/S(L) GO TO 214 212 PRINT 213,L 213 FORMAT ( IX,30HSTRUCTURE IS UNSTABLE AT DISP U ) STOP 5 4 3 6 214 CONTI NUE 228 CHECK EACH RECTANGULAR ELEMENT FOR STABILITY DO 222 K= 1 » NRM XL AM = ASP ( K)#EX / ( lo-UXY*UYX ) XGAM = G/ASP(K) YLAM = EY/(ASP(K)*<l.-UXY*UYX)) YGAM = ASP(K)*G 217 D I AA = ( XL AM + XGAM) *T ( K ) M . IF( A * T ( K ) - D I A A ) 217,217,218 PRINT 219, K 219 FORMAT(IX,13HMEMBER NUMBERI5,32H IS UNSTABLE DUE TO A TOO SMALL) 218 DIAJ = (YLAM + YGAM)*T(K)/4o L F _ i A J * T ( K ) ~ DIAJ) 220,220 ,221 : ; _ : 220 PRINT 223,K 223 FORMAT(IX,13HMEMBER NUMBERI5,32H IS UNSTABLE DUE TO J TOO SMALL) 221 D I A A J = ( ( UXY*XLAM/ASP(K ) +ASP ( K ) * XGAM)* T ( K ) M « ) **2 ; 224 AAJ = A*AJ IF(AAJ#T(K)#T(K)-DIAAJ) 224,224,222 PRINT 225,K \ ...p-_ . 03-225 FORMAT(IX,13HMEMBER NUMBERI5,33H IS UNSTABLE DUE TO AJ TOO SMALL) 2 22 CONTINUE PRINT 2 26 ._. Y 12 11 -10 9 3 7 6 . 5 226 FORMAT(IX,52HALL RECTANGULAR ELEMENTS NOT STATED ABOVE ARE STABLE 227 CONTINUE PRINT 536,NRS . 536 FOR~MATI: 1H1,37HSTR _ESS ANALYSIS FOR STRUCTURE NUMBER 14) DO 457 KASE = 1,NLC DO 537 L=1,NU 537 472 538 "539" FK(L)=0.0 READ 472, NJL FORMAT (110) _.. PRINT ~ 5 3 8 » K A S E , N J L FORMAT(IX,15HLOAD CONDITION 14,6H HAS I4,16H LOADED JOINTS.) PRINT 539 .._ FORMAT(IX,32H JOINT NO. X FORCE Y FORCE) DO 540 J=1»NJL • " _ : READ 550 9 JNIUFF ( 1 ) »FF ( 2 ) 550 FORMAT ( I10»2F.10o3) 541 FORMATt 1X9 I1Q92F10.3) - ; PRINT 5419 J N L U F F d ) »FF(2) DO 540 L = l 9 2 I = ND(JNU.L) ; . - ' I F ( I ) 5429542,543 542 I F ( F F ( L ) ) 5 4 4 , 5 4 0 * 5 4 4 544 PRINT _ _ L i L _ J j _ J _ , — : -545 FORMAT(IX 9 33H STRESS FROM FORCE I2.11H AT JOINT 14. 115H NOT INCLUDED.) . • . GO TO 540 - . __ 543 F K ( I ) = FK( I )+FF(L) 540 CONTINUE C BACK SUBSTITUTE ONE ; ' 310 J2 = (NB*NU)+1 . J3=(NB+1)*NU DO 311 I=1»NU ; : ; . KK= J2+I-1 311 S(KK)=FK(I) S ( J2 ) =S ( J2 )VS( 1 ) ; • - DO 302 L=Ll9L29NB J l = J 2 + ( ( L - l ) / N B ) J = J1-K3 : : : K1=L-NBS IF(K1) 30493049303 304 K1 = 1+ ( ( L - 1 ) / N B ) J = J2 303 K2=L-K3 DO 301 K = K19_K2_9j<3__ : S ( J l ) = S ( J l ) - S ( K ) * S ( J ) 301 J=J+1 3_o_2 _sj. J3 JL f__ (Ji) /s( L > '_ ; ; C BACK SUBSTITUTE TWO L3=L2+1 S( J 3 ) = S ( J 3 ) / S ( L 2 ) DO 308 LR=L1,L2»NB L= L3-LR J l = J2+( ( L - D / N B ) J =J1+ 1 K1=L+1 K2=L+K3 DO 309 K=K1,K2 S ( J 1 ) = S ( J 1 ) - S ( K ) * S ( J ) 309 308 J = J+1 S ( J 1 ) = S ( J l ) / S ( L ) DO 312 1 = 1, Nil 312 C KK=J2+I-1 D(I)=S(KK) BACK SUBSTITUTE THREE 400 PRINT 400 FORMAT(1H1,17HJ0INT DEFLECTIONS) PRINT 401 o 1 401 C FORMAT(IX,50HJOINT NUMBER X-DEFLECTION JOINT DEFLECTIONS DO 501 K=1,NRJ Y-DEFLECTION) N2 = ND(K»1) N3 = ND(K,2) IF(N2) 408,409,410 T 12 408 N4 = -N2 DX = -D(N4) 60 TO 500 JL ••10 9 409 410 DX = 0.0 60 TO 500 DX = D(N2) 3 7 6 500 411 IF(N3) 411,412,413 N5 = -N3 DY = -D(N5) 5 60 TO 501 4 3 412 DY = 0 . 0 GO TO 501 413 DY = D(N3) 501 PRINT 414,K,DX,DY 414 FORMAT(IX,I7,9X,1PE14.7,4X,1PE14.7) PRINT 415 415 FORMAT(1H1 ,14HMEMBER STRAINS) PRINT 416 416 FORMAT(IX,83HMEMBER SHEAR STRAIN 1 ' STRAIN STRAIN) STRAIN PRINT 417 417 FORMAT(IX,83HNUMBER STRAIN AT TOP 1 AT RIGHT AT LEFT ) AT BOTTOM C MEMBER STRAINS DO 418 K=l,NRM DO 421 1 = 1 ,8 i IF(NP(K,I ) )423,424,425 424 DM(I) = 0.0 GO TO 421 1 423 NPI = -NP(K,I) DM(I) = -D(NPI) GO TO 421 425 NPI = NP(K, I) DM(I) = D(NPI ) 421 CONTINUE 426 GAM(K) = 0.5*((DM(1)+DM(4)-DM(5)-DM(8))/ASP(K)+DM(2 1-DM(6))/HORL(K) EXT(K) = ( D M ( D - DM(4))/HORL(K) +DM(7)-DM(3) EXB(K) =(DM(8) - DM(5))/HORL(K) EYR(K) -=(DM(2) - DM(7))/(ASP(K)#HORL(K)) EYL(K) =(DM(3) - D M ( 6 ) ) / ( A S P ( K ) * H O R L ( K ) ) PRINT 428,K,GAM(K) »EXT(K)»EXB(K) »EYR ( K ) » E YL ( K )• 428 FORMAT(IX,14 ,4X,1PE14.7,1P4E16.7) GO TO 418 418 CONTINUE • i ) ,3 i i PRINT 440 I • 1 I 1 440 441 FORMAT(1H1,12HN0DAL FORCES) PRINT 441 0FORMAT(1X»126HMEMBER F l F2 F3 • 1 F4 F5 . F6 F7 2 F8) PRINT 442 442 C FORMAT ( I X , 6HNUMBER) MEMBER FORCES DO 448 K=l , NRM t 494 DO 492 1=1,8 ' IF(NP(K,I)1493,494,495 • DM(I) = OoO 493 GO TO 492 NPI = -NP(K,I) DM(I) = -D(NPI ) i 495 GO TO 492 NPI = NP(K , I ) DM(T) = D(NPI ) N) 1 -492 CONTINUE DO 446 J = l , 8 N = ( ( J * ( 1 7 - J ) ) / 2 ) - 8 446 DO 446 I=J,8 IN = I+N S A ( I , J ) =• SMM( K, IN ) 1 T 12 447 DO 447 1=1,8 DO 447 J = l , 8 S A ( I , J ) = S A ( J , I ) 2 . . _ _ ! ! ' iO 9 DO 452 1=1,8 F ( I ) = 0.0 DO 452 J = l , 8 3 4 3 7 6 452 448 453 F ( I ) = F ( I ) + S A ( I • J ) * D M ( J ) PRINT 45 3 , K , F ( 1 ) , F ( 2 ) , F ( 3 ) , F ( 4 ) , F ( 5 ) , F ( 6 ) » F ( 7 ) , F ( 8 ) FORMAT(IX,14,1PE15.7,1P7El 6.7 ) 5 5 4 3 PRINT 455 o 455 FORMAT ( 1H 1 » 27HMEMBER STRESSES AT CENTROID) PRINT 456 456 FORMAT(1X»61HMEMBER NUMBER X-STRESS Y-STRESS IR STRESS) SHEA PAR = l.-UXY*UYX .1 DO 457 K=1»NRM SXY(K) = GAM(K)*G 458 SX(K) = EX*(EXT(K)+EXB(K)+UXY*(EYL(K)+EYR(K) ) )/(2o*PAR ) SY(K) = EY * ( E Y L ( K ) + E Y R ( K ) + U Y X * ( E X T ( K ) + E X B ( K ) ) ) / ( 2 o * P A R ) 466 PRINT 467»K*SX(K)» SY(K)? SXY(K) 467 FORMAT( 1 X » I 7 » l P E 2 1 o 7 » l P 2 E 1 7 o 7 ) f 457 CONTINUE GO TO 999 STOP END , i ; vn — 1 12, 11 ' 2 10 9 . 3 8 . 7 4 6 . 5 5 4 3 . 6
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A stiffness model for plane stress analysis
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A stiffness model for plane stress analysis Hibbert, Paul D. 1965
pdf
Page Metadata
Item Metadata
Title | A stiffness model for plane stress analysis |
Creator |
Hibbert, Paul D. |
Publisher | University of British Columbia |
Date Issued | 1965 |
Description | The stiffness properties of a finite sized square model element of an isotropic plate under plane stress are presented. The model is a piece of the material itself and not a configuration of bars that replace the material. It is shown that the use of an extra kinematic condition instead of an arbitrary displacement function usually produces better results. An example shows that this extra condition greatly increases the accuracy of the model and so reduces the number of degrees of freedom required to approximate the continuous system. A system of bounding the solution is shown by assigning values to free parameters. Thus the maximum error is known and can be reduced without going to a fine grid. The derivation of the isotropic square is then expanded to obtain the stiffness matrix of an elemental orthotropic rectangle. The limits of stability for both the isotropic and orthotropic models are investigated. |
Subject |
Strains and stresses |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-09-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050614 |
URI | http://hdl.handle.net/2429/37548 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1965_A7 H5.pdf [ 2.93MB ]
- Metadata
- JSON: 831-1.0050614.json
- JSON-LD: 831-1.0050614-ld.json
- RDF/XML (Pretty): 831-1.0050614-rdf.xml
- RDF/JSON: 831-1.0050614-rdf.json
- Turtle: 831-1.0050614-turtle.txt
- N-Triples: 831-1.0050614-rdf-ntriples.txt
- Original Record: 831-1.0050614-source.json
- Full Text
- 831-1.0050614-fulltext.txt
- Citation
- 831-1.0050614.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050614/manifest