A STIFFNESS MODEL FOR PLANE STRESS ANALYSIS by PAUL D. HIBBERT B.A.Sc. ( C i v i l Eng.) The University of British Columbia, 1958 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF ' • MASTER OF APPLIED SCIENCE i n 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 i n 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 purposes may be granted by the Head of my Department or by his representatives. I t 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 s i z e d square model element of an i s o t r o p i c p l a t e under plane s t r e s s are presented. The model i s a piece of the m a t e r i a l i t s e l f and not a configuration of bars that the m a t e r i a l . replace I t i s shown that the use of an e x t r a kinematic c o n d i t i o n i n s t e a d of an a r b i t r a r y displacement f u n c t i o n u s u a l l y produces b e t t e r results. An example shows that t h i s e x t r a c o n d i t i o n g r e a t l y increases the accuracy of the model and so reduces the number o f degrees of freedom required to approximate the continuous system. A system of bounding the s o l u t i o n i s shown by assigning values t o free parameters. Thus the maxi- mum e r r o r i s known and can be reduced without going to a f i n e g r i d . The d e r i v a t i o n of the i s o t r o p i c square i s then expanded t o obtain the s t i f f n e s s m a t r i x of an elemental o r t h o t r o p i c rectangle. The l i m i t s of s t a b i l i t y f o r both the i s o t r o p i c and orthotropic models are i n v e s t i g a t e d . - iii - ACKNOWLEDGEMENTS The author wishes to express his gratitude, to Dr. R.F. Hooley for his invaluable guidance during the research and i n the preparation of this thesis. Also, the author is grateful to the National Research Council of Canada for financial support i n the form of an assistantship, and to the staff of the Computing Centre at the University of British Columbia for their help. A p r i l , 1965 Vancouver, B r i t i s h Columbia - ivTABLE OF CONTENTS Abstract A cknowle dgement s Contents Notation Page I. II. Introduction 1 I s o t r o p i c Square 3 2.1 2.2 2.3 2.4 2.5 III. IV. Derivation C r i t i c a l Value of a f o r Isotropy Bounding Technique Example The Determination of Stresses from Known D e f l e c t i o n s 3 7 8 10 11 Orthotropic Rectangle 14 3.1 Derivation 14 3.2 C r i t i c a l Values o f a and P f o r Orthotropy 16 Conclusions 17 Table I 20 14 Figures 21 Bibliography 35 Appendix A . l D e s c r i p t i o n of the Computer Program, P0RT I A.2 S i m p l i f i e d Flow Diagram A.3 F o r t r a n IV L i s t i n g 36 40 41 V Notation I n some cases the computer program symbols vary from those defined here. The major d i f f e r e n c e s are pointed out i n the d e s c r i p t i o n of the program, which i s t o be found i n S e c t i o n V. The word "entry" i s used throughout the t e x t t o r e f e r t o a matrix element i n order t o avoid confusion w i t h 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 x y - modulus of e l a s t i c i t y i n subscripted d i r e c t i o n F «• force vector F i (i=l,2,3,...8) - nodal force G <= shear modulus I » moment o f i n e r t i a L « M = bending moment R » aspect r a t i o ( v e r t i c a l dimension t o h o r i z o n t a l ) S » s t i f f n e s s matrix of element T. • transformation matrix h o r i z o n t a l dimension of model t thickness a, 0 diagonal e n t r i e s o f S shear s t r a i n Gt/4 R /x RGt/4 A SEE d e f l e c t i o n vector t r i a l d e f l e c t i o n vector A' A transposed - viA. (i-1,2,3,. ,8) = nodal d e f l e c t i o n A6, - angle changes - s t r a i n vector « s t r a i n i n subscripted d i r e c t i o n X - Et/2(l-n ) Ax - A3 ( ^xy* ^yx 2 B E x t A ( l W y X ) • Poisson*s r a t i o f o r i s o t r o p y • 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 , s t r e s s i n second s t r e s s vector o~" , o x' y ^x <y cr- s t r e s s i n subscripted d i r e c t i o n shear s t r e s s Rt E /12 x' t E /12R y A S t i f f n e s s Model For Plane Stress A n a l y s i s I. Introduction Various models have been proposed f o r the approximate s o l u t i o n of the plane stress problem, notably those of Hrennikoff ( l ) , McCormick a Melosh (3), Clough (4), and Turner et a l . ( 5 ) . (2), I n a l l cases the p l a t e i s divided i n t o a g r i d of a f i n i t e number of elements. I n the a n a l y s i s 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 w i t h i n 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 m a t e r i a l they replace, whereas Melosh, Clough, Turner et al.consider the element t o be a piece of the m a t e r i a l i t s e l f , and with the a i d of energy c r i t e r i a they obtain more accurate stiff- ness matrices than were a v a i l a b l e 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 t h e s i s f o l l o w s the procedure o u t l i n e d by Turner et a l . of applying u n i t deformations successively to the eight degrees-of-freedom (two t r a n s l a t i o n a l d e f l e c t i o n s are assumed at each node of a r e c t a n g l e ) . As Turner pointed out, the same r e s u l t can be from s t r a i n energy considerations. obtained I n e i t h e r 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 t i o n s due to t e n s i o n , shear, and bending. In (4) Clough presented i n d e t a i l the s t r a i n energy d e r i v a t i o n , based on s u p e r p o s i t i o n , of a s t i f f n e s s matrix of an i s o t r o p i c 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 t e n s i o n and shear. L a t e r (6) he obtained a s t i f f n e s s matrix of an ortho- t r o p i c , t r i a n g l e but d i d not present r e s u l t s . Numbers i n brackets r e f e r to the bibliography. Instead of obtaining the general d i s t o r t i o n of the element by superp o s i t i o n , Melosh generated a displacement f u n c t i o n . The f u n c t i o n gave the displacement of points w i t h i n the element as ordinates to a hyperbolic paraboloid. He presented a s t i f f n e s s m a t r i x based on the above displacement f u n c t i o n and found that i t s a t i s f i e d macroscopic e q u i l i b r i u m and the r e q u i r e ments of monotonic convergence ( i . e . as the c e l l s i z e i s reduced the s o l u t i o n i s never f a r t h e r from the exact s o l u t i o n ) . The use of Melosh's matrix, as w i l l be shown, leads to r e s u l t s which are f a r t h e r from the exact s o l u t i o n than Clough's, although i t i s nearer t o the exact than H r e n n i k o f f s or McCormick's. A l l of the aforementioned matrices y i e l d solutions which converge on the actual s o l u t i o n from below, that i s , the d e f l e c t i o n s are too s m a l l . Therefore i f a s t i f f n e s s m a t r i x i s used which y i e l d s d e f l e c t i o n s which are too l a r g e , the solutions w i l l converge on the a c t u a l s o l u t i o n from above. I f both of the curves converge monotonically then the true s o l u t i o n must l i e between them and thus the maximum e r r o r i s known f o r the p a r t i c u l a r g r i d s i z e used. I t is shown i n t h i s t h e s i s that i f the bending c r i t e r i o n is deleted and i f instead some of the diagonal e n t r i e s of the s t i f f n e s s m a t r i x are considered as free parameters, matrices y i e l d i n g any magnitude of d e f l e c t i o n can be obtained; hence the true s o l u t i o n can be bracketed, or bounded from above and below. - 3 II. I s o t r o p i c 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 a c t i n g as were present i n the loaded p l a t e , 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 s m a l l , t h e s t r e s s 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 o f F i g . 1(c) i s replaced by the s t a t i c a l l y equivalent s e t 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 t o determine i f the forces a r i s e from shear or normal s t r e s s , or from a combination of the two. S i m i l a r l y the actual d e f l e c t i o n s o f the square Q as depicted i n F i g . i n Fig. 1(f). 1(e) can be approximated as l i n e a r , as shown The d e f l e c t i o n 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 m a t r i x S w i l l govern the behaviour of t h i s element such that • Where SA F and A - F (1) are eight-component vectors which represent the forces and d e f l e c t i o n s of F i g s . 1(d) and 1 ( g ) . The f i r s t column of S forces e x i s t i n g at the nodes when A^ = 1 •zero. the represents the and a l l other displacements are F i g . 2(a) shows the element d i s t o r t e d under these forces. second column of S contains the forces shown i n F i g . and a l l the r e s t are zero. element of F i g . Similarly 2(b) when A 2 = 1 These l a t t e r forces are obtained by r o t a t i n g the '. 2(a) about i t s Z - a x i s . The remaining s i x columns can be derived from the f i r s t two by r o t a t i n g the elements about an a x i s normal t o the paper. The e n t r i e s 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 eight unknowns. be Various r e l a t i o n s h i p s between the e n t r i e s of S only will now established. Because S must be symmetrical about i t s main diagonal c - (2) g Three more r e l a t i o n s h i p s are obtained from e q u i l i b r i u m 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) I n order that the model represent the r e a l m a t e r i a l , i t must e x h i b i t the same e l a s t i c properties. Consider f i r s t a model element of length subjected to a uniform t e n s i l e stress o— . y L F i g . 4 shows t h i s stress gathered t o the nodes as concentrated forces, and the nodal deflections which are derived from Hooke's Law. The f o l l o w i n g equations are obtained when the r e c i p r o c a l theorem i s w r i t t e n 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) —Zo- L -(b + c) u. 'Where E and t and > o-L + (a + h) -f- - 0 = —|— (6) o-Lt (7) are the e l a s t i c modulus and Poisson*s r a t i o r e s p e c t i v e l y i s the plate thickness at element Q. Consider next the model element under the a c t i o n of a uniformly buted shear s t r e s s , T, distri- as shown i n F i g . 5 with the stresses gathered to form f o r c e s at the nodes. I n t h i s case A. = - A. = T L / G , where G i s the shear modulus. Equation (8) i s obtained when the r e c i p r o c a l theorem i s w r i t t e n 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 e q u i l i b r i u m and has the same e l a s t i c properties as an element of the r e a l m a t e r i a l , there are, only seven equations l i n k i n g the.eight unknowns o f S. However, a s o l u t i o n can be obtained i n terms o f one unknown, say a. I f t h i s i s done, then X 4 X 4 X 4 X 4 where (1 + u.) (9) (3u'-l) (10) (2|A - 2) + a (11) (6 - 2u.) - a (12) f - b g - c h = •X - / (13) (14) a (15) 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 c o n f i g u r a t i o n of bars f o r h i s model provided the remaining l i n k between the eight unknoxvns and t h i s y i e l d s yields a = X ( 3 0 - 18u,)/24. S i m i l a r l y , McCormick's choice of bars a = A ( 3 3 - 27u.)/24. Neither of these two models, however, can a c c u r a t e l y represent the a c t u a l deformed state of the element because no combination of a square d i s t o r t e d t o a rectangle (no angle changes at the corners) and a square d i s t o r t e d t o a - 6 parallogram (equal angle changes at the corners) can form the general d e f l e c t e d s t a t e shown i n F i g . 6, where A6 £ A(3. I n other words, i n order to define a generally d e f l e c t e d element i t i s necessary t o 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 l e n g t h changes and three nodal angle changes. Extension i n the x d i r e c 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 t o i s o t r o p y i t i s i d e n t i c a l t o t h a t of x extension. pendent angle change. Shear d i s t o r t i o n produces one inde- Two displacements remain to be found. These kinematic d e f i c i e n c i e s may be removed by means of a t h i r d case of the element d i s t o r t e d i n bending. Bending stresses a p p l i e d separately t o the v e r t i c a l and h o r i - zontal sides produce two independent changes, but due t o i s o t r o p y they are equal. To derive the bending case r e f e r t o F i g . 7, which shows that under an a p p l i e d moment, where I M, the angle change between the v e r t i c a l sides i s ML/EI, i s the moment of i n e r t i a . Then and the remaining d e f l e c t i o n are zero. and the others are zero. _^ = A^ = -A_ = - A g = ML /4EI As w e l l , F^ * F^ = - F^ = - Fg => M/L Equation (16) i s obtained when the r e c i p r o c a l theorem i s w r i t t e n between t h e conditions of F i g . 2(a) and F i g . 7. ' (a + d - e - h) g - - I (16) Note t h a t t h i s e x t r a equation i s necessary and s u f f i c i e n t t o define the s t i f f ness matrix o f the model element. The s o l u t i o n o f the f u l l s e t of eight equations i s set f o r t h i n the second column of Table I and the r e s u l t s are compared w i t h those of Hrennikoff, McCormick, and Melosh. L a t e r 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 b e t t e r r e s u l t s than the other three. 2.2 C r i t i c a l Value of a f o r Isotropy Another way t o remove the kinematic d e f i c i e n c i e s i s t o consider a f r e e parameter rather than f i x i t s value by a bending c r i t e r i o n . bound on a o as A lower can be found by demanding that energy be stored rather than released as the model i s s t r a i n e d . This requires that f o r any t r i a l vector A, the f o l l o w i n g holds t r u e A 'S A I n other words, S ^ 0 (17) must be p o s i t i v e s e m i - d e f i n i t e . This i m p l i e s t h a t a l l eight upper l e f t i x i subdeterminants must be greater than or equal t o zero. This s t i p u l a t i o n a p p l i e d t o the upper l e f t l x l subdeterminant y i e l d s and t o the upper l e f t 2 x 2 subdeterminant y i e l d s a ^ 0 a > b. To avoid the e v a l u a t i o n of higher order subdeterminants, i t i s convenient t o change, t e m p o r a r i l y , the vector numbering convention e s t a b l i s h e d i n F i g . 1(d) and 1 ( g ) . A d i f f e r e n t a r b i t r a r y sequence of numbering the f o r c e and d e f l e c - t i o n vectors i n no way changes the r e s u l t of the d e r i v a t i o n , provided that consistency i s observed. Except f o r the purpose of s t a b i l i t y i n v e s t i g a t i o n s , the numbering convention of F i g . 1(d) and 1(g) i s observed throughout the theory. Consider vectors 2 and 3 interchanged. this basis, S interchanged. I f t h e d e r i v a t i o n proceeds on remains the same except that rows and columns 2 and 3 are This manipulation, however, has no e f f e c t on the a l g e b r a i c s i g n o f any of t h e subdeterminants. s t a b i l i t y c r i t e r i a as t h e o l d ; must be ^ 0, therefore The new matrix i s governed by t h e same hence the upper l e f t 2 x 2 subdeterminant a ^ c. Thus be interchanging the remaining vector - 8 numbers (4 t o 8) w i t h vector 2 and evaluating the upper l e f t 2 x 2 subdeterminant i n each case, a l l e i g h t s t a b i l i t y l i m i t s are obtained, two of which are duplicated. The s o l u t i o n i s : (18) a>Q For |i « 1 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 a > e a ^ X ( 3 - u)/4 (22) a > h a > (23) (n < l ) (21) X /2 the above l i m i t a t i o n s are a l l contained i n a ^ X/2 and f o r any smaller p o s i t i v e value o f u. a l l l i m i t s are contained i n a ^ e. the most severe r e s t r i c t i o n i s cc ^ M 3 - n)/4 Therefore and a must be the l a r g e s t entry i n S. If a = X (3 - u)/4 approaches i n f i n i t y . then the system i s unstable and the d e f l e c t i o n I t i s obvious from i t s d e f i n i t i o n t h a t as a approaches i n f i n i t y , the d e f l e c t i o n approaches zero. Therefore, some value of a greater than X ( 3 - u.)/4 2.3 must produce t h e exact d e f l e c t i o n . Bounding Technique I f the e r r o r decreases as the s i z e of the c e l l s decrease, t h i s e r r o r can be bounded by s o l v i n g the problem f o r several s i z e s of c e l l and two values of a. I n doing so, the s t i f f n e s s m a t r i x must be formed using equations (9) t o (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 t o be greater than t h a t r e q u i r e d f o r an exact s o l u t i o n and the other t o be l e s s . Those solutions f o r a l a r g e r than that of the exact w i l l e x h i b i t d e f l e c t i o n s which increase w i t h a decrease i n c e l l s i z e , and so approach the exact s o l u t i o n from below. Solutions f o r a smaller than the exact s o l u t i o n w i l l show d e f l e c t i o n s which decrease w i t h a decrease i n c e l l s i z e and approach the t r u e s o l u t i o n from above. Since the exact s o l u t i o n i s approached from above and below, the d e f l e c t i o n i s bounded and the maximum e r r o r i s known. the two values of a I f both of used produce d e f l e c t i o n s which approach the t r u e from the same s i d e , then one value must be changed so as to approach from the other s i d e . An obvious advantage i n t h i s system i s that there i s no need t o keep decreasing the g r i d s i z e i n order t o 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 produces the exact d e f l e c t i o n at a l l p o i n t s . a a I n f a c t , a d i f f e r e n t value of i s required t o produce the exact s o l u t i o n at each point. Having shown t h a t d e f l e c t i o n s may be bounded, i t i s easy t o show t h a t s t r a i n s may also be bounded. Consider, f o r example, a system i n which the d e f l e c t i o n of node number 1 i s bounded by two values of 2 by two d i f f e r e n t values. a and that of node Suppose that the maximum e r r o r s , as found by bounding, i n the d e f l e c t i o n s of points 1 and 2 are 6^ and 6^ respectively, then the maximum e r r o r i n s t r a i n between the two points i s (16^1 + |6 I)/L. 2 When expressed as percentage e r r o r s , the e r r o r i n s t r a i n may thus be greater than e i t h e r e r r o r i n d e f l e c t i o n . 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 d e f l e c t i o n s of f o u r points i n s t e a d o f two, i t can be shown t h a t stresses may also be bounded. Melosh (3) gives a procedure f o r bounding d e f l e c t i o n s at a point based on the f a c t that the s t r a i n energy o f the exact s o l u t i o n l i e s between the minimum p o t e n t i a l energy and the minimum complementary energy. To o b t a i n the maximum e r r o r using t h i s concept i t i s necessary to o b t a i n s o l u t i o n s f o r r e a l and influence loads f o r two models; one model must be based on continuous d e f l e c t i o n s over the body and the other on continuous stresses. Melosh s 1 technique implies t h a t the s o l u t i o n from the continuous d e f l e c t i o n model produces a d e f l e c t i o n which i s too s m a l l , whereas that from the continuous s t r e s s model produces a d e f l e c t i o n which i s too l a r g e . Therefore the t r u e - 10 d e f l e c t i o n i s bounded from above and below. I n order to reduce the e r r o r found by Melosh's method i t i s necessary t o reduce the g r i d s i z e . This method of bounding d e f l e c t i o n s at a p o i n t , although complete i n theory, i s r e s t r i c t e d due to the l a c k of s t r e s s consistent matrices and so does not yet have p r a c t i c a l a p p l i c a t i o n . A search of the l i t e r a t u r e shows no reported s t r e s s consistent matrices or attempts at bounding by Melosh's technique. 2.4 Example A c a n t i l e v e r beam of rectangular c r o s s - s e c t i o n with a length-to-depth r a t i o of three was run w i t h the computer program. PORT I , which i s described i n the Appendix. F i g . 8 shows the end d e f l e c t i o n of t h i s beam v s . 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 b e t t e r accuracy than the other two models but there i s s t i l l some e r r o r , when compared w i t h the e l a s t i c i t y s o l u t i o n f o r no f i x e d end warping. This i s due p r i m a r i l y to the r e l a t i v e l y crude approximation of the r e q u i r e d parabolic shear d i s t r i b u t i o n between nodes. The two curves f o r McCormick are f o r 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 e x i s t s at the node and i t i s i n t e r e s t i n g to note t h a t when run i n t h i s manner the r e s u l t s are s l i g h t l y b e t t e r . F u r t h e r , there i s one degree of freedom l e s s per j o i n t when run t h i s The curves show t h a t a/X m .861 way. i s too l a r g e f o r f i n d i n g the end d e f l e c t i o n 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 s o l u - 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 t h a t o/X « .85 should produce the exact s o l u t i o n 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 f a c t t o w i t h i n 0.55$. F i g . 10 shows a p l o t of the extreme f i b r e bending s t r e s s h a l f 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 t h a t the s t r e s s i s bounded above and below by d i f f e r e n t values o f a//\ . I n t h i s case the value of a/X bending s o l u t i o n i s required to produce the exact r e s u l t . independent of the number of c e l l s h i g h . = .861 from the This again i s Also the parameter a required f o r exact s t r e s s i s d i f f e r e n t from t h a t required f o r exact d e f l e c t i o n , as would be expected. 2.5 The Determination of Stresses from Known D e f l e c t i o n s I n order t o determine t h e stresses from known corner d e f l e c t i o n s , as i n F i g . 11(a) f o r example, i t i s f i r s t necessary t o i n v e s t i g a t e s t r a i n s . From the u n i t cases of t e n s i o n , shear, and bending used i n the model d e r i v a t i o n , i t was found, over the area of the element, t h a t (. and 6 ' • x y vary l i n e a r l y under bending moments a p p l i e d to the v e r t i c a l and h o r i z o n t a l sides r e s p e c t i v e l y , t h a t they are constant under t e n s i o n , and that constant under shear. V is Accordingly then, the s t r a i n s a t any point w i t h i n the element may be defined by plane surfaces as shown i n F i g . 11(b), 1 1 ( c ) , and 11(d). I t i s assumed t h a t t h e element under c o n s i d e r a t i o n i s i n the i n t e r i o r of .a s t r u c t u r e ; the diagrams. the s t r a i n surfaces of three adjacent members are shown i n Note should be taken of the f a c t t h a t ( and £ * on l i n e s p a r a l l e l t o the x and y are constant y axes r e s p e c t i v e l y . I n F i g . 11 there are d i s c o n t i n u i t i e s o f the s t r a i n surfaces at the e l e ment boundaries. I t i s deemed advisable therefore t o evaluate s t r a i n s (and stresses) at the centre of the element i n order t o avoid the d i s c o n t i n u i t y e f f e c t s , and a t the same time, give equal weight t o the four corner v a l u e s . - 12 When the s t r a i n s have been evaluated, stresses f o l l o w from the w e l l known r e l a t i o n s h i p s : o- tf + - x n£ ) E/(l -jx ) (< o- (24) 2 x y (25) + ^e )E/(i - y x (26) The stresses a t a boundary, however, can be evaluated more d i r e c t l y . The normal stress perpendicular t o the boundary and the shear stress are equal t o the applied s t r e s s . The stress p a r a l l e l t o the boundary i s found from E i o-' X + u o— •* y x (27) f o r 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 d e f l e c t i o n s , 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 r e l a t i o n s h i p between d e f l e c t i o n s and s t r a i n s i s £x ^ => T A, where € and T^ i s an 8 x 3 transformation matrix, •y 1 then by geometry: <x V - k " ( A 2L" ( A k ( A 1 2 1 + A + A + A 4 3 + + 5 A A 6 + (28) V + A 7 (29) } 2 - 3 - \ A + A 5 + A 6 " 7 " V A ( 3 0 ) therefore 0 2L 0 1 1 1 -1 1 0 -1 1 0 0 0 1 1 1 1 -1 1 0 (3D -1 - 13 - by Secondly, i f the r e l a t i o n s h i p between s t r a i n s and d e f l e c t i o n s i s given oX and i s an 8 x 3 transformation o- = A, where cr- m a t r i x , then by equations (24), (25), (26) E - 2L(1V E y [(A + A x [(A 2 + A 4 3 + A 5 + Ag) + ( A + A + A 6 • A ) U.J (32) + A 6 .+ A ) + ( A + A^ + A 5 + A ) uT] (33) 2 ? 3 1 ? g 2L(l-u. ) 2 E (34) 4(l+u)L therefore 1 u. H l 1x n (j. 1 1x |i 2 2 2 2 2 2 2 (35) 2L(l«si ) T h i r d l y , because £ ** T^ A and o- = T^ A s t r e s s and s t r a i n may be w r i t t e n as formation matrix and i s defined by 1 o- = T^6 , T u 2 « T^ T^. then the r e l a t i o n s h i p between where T^ i s a 3 x 3 trans- I t i s found that 0 E 1-u* 2 (36) 0 o i f - 14 III. Ortho-tropic Rectangle 3.1 D e r i v a t i o n The foregoing approach f o r f i n d i n g 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 a n a l y s i s of a plane o r t h o t r o p i c continuum approximated by rectangles (see F i g . 12). The associated s t i f f ness matrix, pic S, ; i s shown i n F i g . 13. Since there i s no longer any i s o t r o - symmetry, the second column i s not a permutation of the f i r s t , and so, i n general, there are s i x t e e n unknowns i n S. Note should be taken, though, of the f a c t that columns 4, 5, and 8 are permutations of column 1 and that columns 3, 6, and 7 are permutations o f 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 t e n s i o n and shear. The s o l u t i o n o f the fourteen equations mentioned above, i n terms of the two diagonal e n t r i e s , b - ^ L ^ c - Y " d - X + - a - 2Y e .•- 2 X x a R Y R V and f$, (37) x x (38) (39) x + 2 V is: x - a . (40) (41) g - c h - 2 X i = b k = 2 A 1 - c (42) x - a (43) (44) - p (45) (46) - 15 - m' - b n - . 2X y + 2 V V P - 2 V p where >v *x - RE * -rrz (47) y - p- (49) y c (50) t r E t T_7T~^ X 4R The f i r s t subscript (48) ' *y of n r 4 r e f e r s t o the d i r e c t i o n of s t r a i n and the second t o the d i r e c t i o n of s t r e s s . R i s the aspect r a t i o and E i s subscripted t o 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 d e f i c i e n c i e s remain which may be removed by applying bending stresses separately t o the v e r t i c a l and h o r i z o n t a l sides. I f t h i s i s done then the complete s o l u t i o n i s : b - ^ p £ c - ^ - d " \ e - V-'^x* x f - b (56) g - (57) h - X + x *x~ V (52) x - R V (53) x x (54) > Y c - r X 1 • p R Y + •-" y X - Y X -<55). : (58) X b + (59) * y + V y ( 6 0 ) - 16 - k =A 1 where - r y - Y (61) y c (62) m" - b (63) n - y = y - V A ° ~ p = Y y . ( T * * 7 - * 7 ' X 3.2 4 ) ( c 6 5 ) (66) Rt E x 6 t E 12 12 R 7 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 f r e e parameters, as was done i n the i s 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 o f the same method of renumbering, previously desc r i b e d , t o avoid the e v a l u a t i o n of subdeterminants o f higher order submatrices than 2 x 2 t o determine p o s i t i v e seml-definiteness, r e s u l t s i n the f o l l o w i n g t h r e e i n e q u a l i t i e s t o define the. lower l i m i t s of a a p > £ s £ ^ R a 3 A l s o i t i s found t h a t A +V x > X +V y x y V ) and p: (67) 2 x . (68) (69) G > 0. A p l o t of the c r i t i c a l values f o r a Douglas F i r element of u n i t aspect r a t i o , and t h i c k n e s s of one i n c h i s shown i n F i g . 14. Values of a l e s s than 545.3 k i p s / i n c h and/or 3 l e s s than 64.3 k i p s / inch produce an unstable element, whereas i n f i n i t e values o f a zero d e f l e c t i o n of the element under load. and 3 cause I t should be noted t h a t 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 m a t r i x 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 superp o s i t i o n of the three b a s i c cases of t e n s i o n , shear, and bending. H r e n n i k o f f s and McCormick's a r b i t r a r y 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 s e l e c t i o n of a displacement f u n c t i o n t h a t s a t i s f i e s the boundary conditions on d i s p l a c e ment i s l i m i t e d to the accuracy of the displacement f u n c t i o n i n representing the t r u e deformation. Melosh chose a very simple f u n c t i o n but the i d e a l one i s probably much more complicated. I t appears that i n order t o extend Melosh's work i t would be necessary to generate displacement f u n c t i o n s of a progressively more complicated form, derive the s t i f f n e s s matrix f o r each one, and compare computer r e s u l t s against known s o l u t i o n s u n t i l a b e t t e r model was obtained. Nevertheless, an analyst would be confronted w i t h an unknown e r r o r , or i f Melosh's bounding technique was used, the e r r o r 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 h e r e i n , the free parameter, a, f o r i s o t r o p y , can be chosen so that the models of Turner, H r e n n i k o f f , McCormick, or Melosh are obtained; or two values of can be chosen so as t o bound the t r u e s o l u t i o n from above and below. a The pro- cedure f o r bounding the d e f l e c t i o n a t a node i s as f o l l o w s : 1. S e l e c t a value of small. a which i s known to produce a d e f l e c t i o n which i s too For a f i r s t t r i a l i t i s recommended t h a t a be chosen to produce the model of Turner et a l . because Turner's model i s known t o be too s t i f f yet gives b e t t e r accuracy than any other. 2. Obtain computer solutions f o r two g r i d s i z e s . Because the value of a chosen i n step 1 produces d e f l e c t i o n s which are too s m a l l , the l i n e j o i n i n g the two s o l u t i o n s on a p l o t of nodal d e f l e c t i o n v s . g r i d s i z e w i l l have p o s i t i v e slope (e.g. the Clough and Hibbert curve of F i g . 8). 3. S e l e c t another value of a so that a d e f l e c t i o n which i s too large w i l l be produced. 4. Obtain computer s o l u t i o n s f o r the same two g r i d s i z e s a 3 i n step 2. In t h i s case the l i n e j o i n i n g the two s o l u t i o n s must have negative slope t o be v a l i d . I f the slope i s p o s i t i v e then the value of a selected i n step 3 was too l a r g e and the d e f l e c t i o n s are smaller than the exact. t h i s case repeat step 4 w i t h a smaller value of i s obtained. a In u n t i l a negative slope Two curves have then been obtained, one of which has p o s i t i v e slope, the other negative (e.g. F i g . 8, the Clough and Hibbert curve and the 5. a/X *= 0.82$ curve). I f the s e l e c t e d values of a are known t o produce s o l u t i o n s which con- verge monotonically, then the a c t u a l nodal d e f l e c t i o n i s known t o l i e between the two curves ( i . e . 114.4 < -p- < 120.6 f o r two c e l l s h i g h , i n Fig. 6. 8). I f the maximum e r r o r so obtained i s unacceptable, show the r e s u l t s on a p l o t of nodal d e f l e c t i o n v s . Fig. 9. a w i t h the g r i d s i z e as a parameter, as i n S e l e c t a more precise value of s e c t i o n of the two c e l l - s i z e curves. obtain a smaller maximum e r r o r . a on the b a s i s of the i n t e r - Obtain s o l u t i o n s f o r the new a and I n t h i s manner the t r u e s o l u t i o n i s always contained between the two curves of smallest p o s i t i v e and smallest negative slope and the e r r o r can be reduced f u r t h e r . I f the s t r a i n along a l i n e j o i n i n g two nodes i s r e q u i r e d , then the values of a must.be s e l e c t e d to bound the s t r a i n ; the c e n t r o i d of an element i s r e q u i r e d , then the bound the s t r e s s . or i f the s t r e s s at a's must be chosen t o - 19 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 Entry Hibbert, Turner, Clough (1) (2) a 22 - 6u - 4u b 6 + 6\x c - 6 + 18^ S a . Hrennikoff McCormick Melosh (3) (4) (5) 33 - 27n 24 - 8u 2 30 - 18u 6 + 6u. 6 + 6\i - 6 + 18u - 6 + 18u - 6 + 18u 18 - 6u 21 - 15u 12 + 4|i d 10 + 6u - e 14 - 6|i + 4n 6 + 6|i 3 + 15H- f 6 + 6^ 6 + 6\i 6 + 6u. g - 6 + 18u h 2 + 6^ + 4 u a V 2 - 6 +• 2 M u l t i p l y a l l e n t r i e s by 18JA - 6 + 18n A/24 12 - > 6 + 6(j. - 6 + 18u - 6 + 18|x - 9 + 27^ 8u - 21 - - 22 - C / / 7 1 e <- f <vf (a) Fig. 2 FORCES FOR COLUMNS I AND 2 OF 23 - S = a b g b a h c h a. d g b e f d e f c h c f e d g b g d e f a h c f e 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 Fig. 3 STIFFNESS MATRIX OF A N ISOTROPIC SQUARE -2k - AcryLt <ryL t "7F F' cry L ~E~ V i 1 1 1 • i L yL • 1 o-yLt Fig. 4 UNIFORM cryLt TENSION -25 - TL H i TLt Ik r y 2 TLt Ay G? / TLt / 1^ TLt TLt 2 Fig. 5 / UNIFORM .SHEAR K "TLt 2 - 26 - Fig. 6 GENERAL DISTORTION - 27 - Fig. 7 BENDING 28 - 3 x depth 0 - 29 - - 30 - NUMBER OF C E L L S Fig. 10 BOUNDING HIGH THE EXTREME FIBRE STRESS. - .31 - (d) THE VARIATION OF STRAINS OVER THE ELEMENT - 32 - RL Ay L Fig.,12 ORTHOTROPIC RECTANGLE - 33 - 1 d e m P h b /3 k c f n 0 g c k ,8 b g 0 n f a S = i d 1 i a h P m e e m p h a i 1 d f n o g b g o n f c k h p m e d I Fig. 13 STIFFNESS k e b i a M A T R I X OF A N RECTANGLE ORTHOTROPIC - 34 Fig. 14 CRITICAL VALUES OF a AND /3 FOR D. FIR - 35 - Bibliography 1. Hrennikoff, A., S o l u t i o n 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 A n a l y s i s , Journal of the S t r u c t u r a l 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 t o Bound E l a s t i c Behaviour of S t r u c t u r e s , Ph.D. Thesis, U n i v e r s i t y of Washington, at S e a t t l e , Wash., August, 1962. 4. Clough, R.W., The F i n i t e Element Method i n Plane Stress A n a l y s i s , 2nd Conference on E l e c t r o n i c Computation, ASCE, September, i960. 5. Turner, M.J., Clough, R.W., M a r t i n , H.C., Topp, L . J . , S t i f f n e s s and D e f l e c t i o n A n a l y s i s of Complex S t r u c t u r e s , 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. : C h p l e s k i , Contributions to the S o l u t i o n of Systems of L i n e a r Equations and the Determination of Eigen Values, E d i t e d by Olga Taussky, N a t i o n a l Bureau of Standards, U.S. P r i n t i n g O f f i c e , Applied Mathematics S e r i e s No. 39, September, 1954, pp. 31. - 36 Appendix A A.l D e s c r i p t i o n of the Computer Program, PORT I Dr. R.F. Hooley's computer program, ALP I I , was r e w r i t t e n from F o r t r a n I I language t o F o r t r a n IV and modified t o solve s p e c i f i c a l l y the problem of plane elasticity. IBM 7040; The program and data are stored e n t i r e l y w i t h i n the core of an therefore tapes are not required and the operating time i s very short. B r i e f l y , the d e s c r i p t i o n of the sequence of operations and the presentat i o n of data required i s as f o l l o w s : 1. 2. A card i s read which bears the s t r u c t u r e data, which c o n s i s t s o f : (a) an a r b i t r a r y s t r u c t u r e 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, t o be analysed. Next, a card i s read which contains the e l a s t i c properties of the m a t e r i a l (i.e. E , u , E , G) x *xy y and then u yx i s c a l c u l a t e d by means of the r e c i p r o c a l theorem. 3. The next data c o n s i s t of a s e r i e s of cards (one card i s r e q u i r e d f o r each joint); each bears: (a) the j o i n t number, JN, (b) e i t h e r a f i g u r e 1 i f the j o i n t i s permitted t o move h o r i z o n t a l l y o r a figure 0 i f not, (c) e i t h e r a f i g u r e 1 i f t h e j o i n t i s permitted to move v e r t i c a l l y or a f i g u r e 0 i f not. I t i s d e s i r a b l e t o obtain a s t r u c t u r e s t i f f n e s s m a t r i x with as narrow a diagonal band as p o s s i b l e i n order t o minimize the computing time and - 37 to allow the analysis of l a r g e s t r u c t u r e s . • To r e a l i z e a narrow band width, c a r e f u l consideration must be given t o 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 a f t e r 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 f r e e parameters d and 0 (A and A J i n the program) f o r u n i t plate thickness are read from a s i n g l e 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, J L L , (e) the lower r i g h t j o i n t number, JLR, (f) the aspect r a t i o (R i n the theory; (g) the thickness of the p l a t e , (h) the h o r i z o n t a l length of the element ASP i n the program), T, and (L i n the theory; HORL i n the program). A f t e r each card i s read, the s t i f f n e s s matrix of each element i s computed and stored, s i n c e , 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 o r i n the thickness of elements. NP, are assigned. Also,, p o s i t i o n numbers, The band width i s c a l c u l a t e d f o r each member i n t u r n but only the maximum value f o r the preceding members i s r e t a i n e d . A f t e r a l l the cards have been read, the computer p r i n t s 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. and i t contains: One card i s r e q u i r e d f o r each r e s t r a i n t (a) the j o i n t number, N22, where an e l a s t i c r e s t r a i n t e x i s t s , (b) the spring type, N23, (a f i g u r e 1 f o r a h o r i z o n t a l s p r i n g ; a figure 2 f o r a v e r t i c a l s p r i n g ) , and (c) 7. the spring constant, ELAST. 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 the program) member by member with the a i d of code numbers, NPS. (S S in is stored as a s i n g l y subscripted v a r i a b l e . 8. The Choleski (3) or Square Root method of s o l u t i o n i s used, whereby the symmetrical band matrix i s decomposed i n t o upper and lower t r i a n g l e s such that S = ZZ'. transpose of Z Z). i s the lower t r i a n g l e and One Z« the upper (Z» i s the 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 h a l f ) 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 d i r e c t systems. 9. The decomposition i s accomplished at t h i s stage. Next, the computer examines Z f o r 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 e n t r i e s f o r a stable s t r u c t u r e must be positive. Zero or imaginary e n t r i e s cause the computer to p r i n t 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 r e j e c t the s t r u c t u r e . 10. Because bounding technique requires several runs with a v a r i a t i o n i n the f r e e parameters, i t may happen, when elements have d i f f e r e n t 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 t o know which of the values of are too s m a l l . a, p, or af3 . The computer p r i n t s the number of the unstable members and the cause of i n s t a b i l i t y but does not r e j e c t the s t r u c t u r e . 11. The computer now reads the number of loaded j o i n t s , NJL, and the l o a d data which consists of NJL cards; (a) each bears: the number of the loaded j o i n t , JNU, - 39 - 12. (b) the force i n the x d i r e c t i o n , F F ( l ) , and (c) the f o r c e i n the y d i r e c t i o n FF(2). Back s u b s t i t u t i o n one i s now accomplished, i n which the computer solves the m a t r i x equation 13. A f o r V, where . F i s now found from the equation V » Z'A, vector. 14. ZV = F i s the load vector. where A i s the d e f l e c t i o n This operation i s c a l l e d back s u b s t i t u t i o n two. Next, back s u b s t i t u t i o n three i s performed. at the end of back s u b s t i t u t i o n two) Joint d e f l e c t i o n s ( a v a i l a b l e are p r i n t e d ; the s t r a i n s on a l l four s i d e s , the eight nodal forces, and the three centroid stresses are c a l c u l a t e d and p r i n t e d f o r each member. 15. I f more than one set of loads was s p e c i f i e d i n step 1, the computer returns t o step 11 t o process the next set. The looping i s repeated u n t i l the l a s t l o a d condition i s analysed. 16. The computer now e i t h e r returns t o step 1 to read data f o r the next s t r u c ture or proceeds t o another program. There i s no conversion of u n i t s w i t h i n the program; presented must be.of consistent u n i t s . therefore a l l data The p r i n t e d output w i l l , of course, be of the same u n i t s as were input. I t i s important t o note that i t i s necessary t o m u l t i p l y rows and columns hi 5, 6, and 7 o f each member matrix, as derived i n the theory, by - 1 t o create a consistent set o f d i r e c t i o n s f o r the s t r u c t u r e . I n the d e r i v a t i o n i t was convenient t o define vectors as p o s i t i v e when they pointed outward from the element but t h i s convention must be abandoned when a structure m a t r i x i s t o be formed. A.2 SIMPLIFIED FLOW DIAGRAM FOR PORT I Start Read and P r i n t ( Q Structure Data i: ^ Read Material Calculate ^yx Print Material Compile Structure S t i f f ness M a t r i x ; Properties) Calculation Convert Properties^ Calculation ' Calculation Ji. Read • Joint and P r i n t Data Assign Calculate ( ( Check Structure Stability Check Element Stability -> . Load Data . Displacement Nos. t C a l c u l a t i o n - B a c k Subst. I e I I Print Joint Deflections Calculate and P r i n t Member Strains No. of Deg. of FreedonQ alfprint^ a n d " ^ 4 Read Member and P r i n t Data i S f o r Member Assign P o s i t i o n Nos. Calculate Band Width to Date Calculate . Nodal Forces and P r i n t of Member Calculate and P r i n t Print J No. of Deg. of ' Freedom.to Date ~5> ( to Z Read and P r i n t Print Calculate A S Band Width Read E l a s t i c Conand P r i n t , s t r a i n t Data - uo - Member Stresses o o > C C C C 999 20 1 PORT I P.DoHIBBERT PLANE ORTHOTROPIC RECTANGLES INPUT AND OUTPUT C O N S I S T E N T U N I T S . HANDLES N E G A T I V E D U P L I C A T E NODE D E F L E C T I O N S . CONTINUE DIMENSION JN( 1 3 5 ) » N D ( 1 3 5 » 2 ) » N P ( 2 0 0 , 8 ) , S M ( 3 6 ) > S ( 8 0 Q 0 ) > S M M ( 2 0 0 » 3 6 ) 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 . 6 A M C 2 0 0 ) » E X T ( 2 0 0 ) » E X B ( 2 0 0 ) . E Y L ( 2 0 0 ) » E Y R ( 2 0 0 ) » T ( 2 0 0 1 DI MFNSI ON S X ( 2 0 0 ) > 5 Y ( 2 0 0 ) , S X Y ( 2 0 0 ) » D ( 2 7 0 ) > F ( 8 ) . S A ( 8 » 8 ) » D M ( 8 ) DIMENSION N P S ( 2 0 0 » 8 ) READ l.NRS»NRJ»NRMfNRE»NLC FORMAT (51 10) , : PRINT 3 » NRS ? NRJ » NRM HAS 1 5 » 9 H J O INTS"» I 5 ,9H MEMBER FORMAT (1H1 ,13HSTRUCTURE N O » I 3 » 6 H 1S ) (1X»3HANDI5i27H ELASTIC CONSTRAINTS WITH I5il8H LOAD COND 1 ITIONS.) •10 9 - 3- - co c+ H* TO hi 1 5 »NRE »N L C FORMAT Ji 3- : PRINT 12. - TlO READ 6 » E X » U X Y * E Y » G FORMAT UF10.3)' UYX = UXY*EX/EY PRINT 7 » E X , U X Y , E Y » U Y X » G FORMAT ( 1 H 0 , 4 H E X = F 10 . 3 » 5H UX Y = F'10 » 3 » 4H 1=F10.3) , NU = 1 PRINT 8 FORMAT ( l'HO .35HJOINT NUMBER X DISP. DO 11 I = 1 » N R J READ 7 9 9 , JN( I )»ND( I ,1)*ND( I » 2 ) 799 FORMAT ( 3 I 10 ) ' PRINT 1 0 » J N ( I ) » N D ( I » 1 ) » N D ( I » 2 ) TO FORMAT (IX»3I10) _• DO 11 J = l i 2 _ NUM = ND( I» J ) NUMA= IABS(NUM) EY = F 10 . 3 » 5H UYX = F 1 0 . 3 » 3 H : Y DISP.) G IF ND(I,J)=NU 13 NU=NU+1 GO TO 11 IF (NUMA-I) 15 16 14 17 18 l l 21 22 2 3 24 25 2 6 hit T . 1 1 • 27 12 J.L 9 3 __ 5 4 , 13 ' ; 14 » 15 , 15 P R I NT 16 »J N J J J ___ FORMAT (1X,23HDISPLACEMENT FOR JOINT I4,19H I S I N WRONG STOP I F (NUM) 17 15 t 18 . •• N D ( I , J ) = -ND(NUMA »J ) GO TO 11 ND(I,J) = ND ( NUMA > J ) " CONTINUE NU = NU-1 N 3 3 = N U_- 2 7 0 . • _ ' IF (N33 ) 23, 23,21 PRINT 22, N33 F O RM A T ( 1 X , 3 Q H T 0 0 M A N Y D E G R E E S OF F R E E D O M BY I4 ) . STOP PRINT 24,NU FORMAT ( IH0__1 Q H T H E R E ARE I5,2QH D E G R E E S OF F R E E D O M ) READ 25»A,AJ PRINT 26,A,AJ FORMAT (2F10.3) FORMAT (1H0,2HA=F10.3,8H AJ=F10.3) PRINT 27 ORDER) T F O R M A T ( 1 H 1 ,61H'MEMBER 33 SM(I) READ G . 12 FORMAT(1H+,42HNUMBER NB = 0 DO 4 7 K=1,NRM DO 3 3 1=1,36 7 5 , 28__ 8 4 11 UR UL LL LR ASPECT IS LENGTH) PRINT 28 "10 2 (NUMA-1) 12 = JOINT JOINT JOINT 0.0 9,MN»JUR,JUL,JLL,JLR»ASP(K),T(K)» JOINT HORL(K) RAT 10) ' THICKNES o o 9 2 I I ! j C 29 i I 1 12 ' 2 ""10 11 c 3 !i 7 4 fi 5 5 '4 3 PRINT 2 » M N , JUR, J U L , J L L , J L R , A S P ( K ) »T (K ) » H 0 R L ( K ) FORMAT < 5 I 5 , 3 F 1 0 o 3 ) FORMAT (1X»I5»4I7»3F10.3) PARAMETERS S'MA = A # T ( K ) SMB = U X Y * E X * T ( K - ) / ( 4 . * ( l o - U X Y * U Y X ) ) + G * T ( K ) / 4 » SMC = SMB - G*T ( K ) /2o. SMD = SMA - G.*T(K')/(2o*ASP(K) ) SME =- SMD + A S P ( K ) ' * E X * T ( K ) / ( 2 « * l l . - U X Y * U Y X ) ) SMH = SME - G * T ( K ' ) / ( 2 o * A S P ( K ) ) SMJ = A J * T ( K ) SMK =-SMJ + F Y * T ( K ) / ( 2 . * A S P ( K ) * ( 1 . - U X Y * U Y X ) ) SMN = SMK + A S P ( K ) # G * T ( K ) / 2 o SMO = SMJ - A S P ( K ) * G * T ( K ) / 2 o S T I F F N E S S MATRIX FOR RECTANGLES SM( 1) = SMA SM ( 2 ) = SMB SM( 3) = SMC SMC 4) = -SMD SM( 5) = -SME SM( 6 ) = -SMB SM( 7) = -SMC SM( 8 ) = SMH . SMJ SM( 9 ) = SMK SM(10) = S M ( 1 1 ) = -SMC S M ( 1 2 ) = -SMB SM(13 ) = -SMN SM(14 ) = -SMO . • SMC SMC15) = SMJ SM(16 ) = SM(17 ) = -SMB SM(18 ) = -SMC SMC 19 ) = -SMO SMC 20 ) = -SMN i o c c 50 40 42 41 4 3. SMB S M ( 2 1 ) .= SM(22) = SMA SM(23) = SMH SM124) = SMC SM(25) = SMB SM( 2 6 ) = -SME SM(27) = SMA. SM(28) = SMB SM(29) = SMC SM(30 ) = -SMD SM( 31 ) =• SMJ SM(32) = SMK S M ( 3 3 ) = -SMC SM(34) = SMJ S M ( 3 5 ) = -SMB SM(36) = SMA C A L C U L A T E P O S I T I O N NUMBERS ND(JUR,1) NP(-K»1) = ND(JUR 9 2 ) NP(K,2) = ND(JUL,2) NP(K»3) = ND(JUL,1) NP(K,4) = NP(K,5) = ND(JLLtl) ND(JLL,2) NP-(K»6J = N D ( J L R 9 2) NP(K»7) = ND(JLR 91 ) NP(K»8) = C A L C U L A T E BAND WIDTH 9 NB MAX = 0 MIN = 500 DO 39 I = 1 »8 N =IABS(NP(K,I)) IF ( N ) 3 9 , 3 9 , 4 0 IF ( N - M A X ) 4 1 , 4 1 * 4 2 • MAX = N IF ( N - M I N ) 4 3 ,39,39 MIN = N "' FOR RECTANGLES i p o o . •" . .... -..v.-.---. . 1 . - -• J 39 CONTINUE N B l = MAX - MIN IF ( N B l - N B ) 4 4 , 4 4 , 4 5 45 NB = N B l DO 46 1 = 1,36 -21 -. 44 . 46 SMM(K»I) = SM(I) CONTINUE 47 NB = NB + 1 MIKF = ?*NB - 1 • PRINT 4 8 , N B , M I K E FORMAT ( 1 H 0 » 3 H N B = I 5 » 2 3 H TOTAL BAND WIDTH = I 5 ) 48 IF (NRE) 51 , 51 » 52 PRINT 53 52 FORMAT ( 1 H 0 , 3 9 H JOINT TYPE E L A S T I C CONSTRAINT) 53 DO 54 I =1,NRE READ 5 5 , N 2 2 , N 2 3 , E L A S T ( I ) PRINT 5 5 , N 2 2 , N 2 3 , E L A S T ( I ) FORMAT ( I X , I 1 0 , I 1 0 , E 1 0 . 3 ) 55 54 NN22( I ) = N D ( N 2 2 , N 2 3 ) CONTINUE 51 F I L L IN STRUCTURE S T I F F N E S S MATRIX C • NV = (NU*NB) + NU + NB - 2 NO = NV - 8000 IF (NO) 56 , 56 , 57 57 PRINT 5 8 , NO 58 FORMAT ( 1 X , 1 1 H 0 V E R L A P OF I 4 , 2 0 H IN S T I F F N E S S MATRIX) STOP MS = NU*NB 56 DO 59 I=1,MS 59 S ( I ) = 0.0 J4 = ((NB+1)*NU)+1 J5= J 4 + N B - 3 DO 60 I = J 4 , J 5 S ( I ) = 0.0 60 N B l = NB - 1 i vn 1 o - ^ :> ' ••• •." 471 62 64 J -;) : 65 66 67 *) 76 68 69 70 72 74 w T ' 12 71 7 73. 75 11 "10 61 9 ' 8 470 7 C 5 . 4 3 469 121 C : .7 . ... . • . DO 471 I = 1 • NRM DO 4 7 1 J = l * 8 NPS(I,J) = NP(I.J) ' DO 61 L = 1 » N R M DO 64 K = 1» 8 N P ( L » K ) = I A B S ( N P S ( L »K) ) DO 61 J = 1 » 8 •IF(NPS(L»J) )65»61»66 AB= - 1 . 0 GO TO 67 AB= 1.0 Jl = (J-l)*(16-J)/2 DO 61 I = J » 8 IF(NPS(L»I)>68,61,69 B = -AB GO TO 70 B = AB I F ( N P ( L » J ) - N P ( L , I ) )71,72,73 IF(I-J)74.71,74 K = (NP(L,I)-1)*NB1 + NP(L,J) IJ1 = I + J l S(K) = S ( K ) + (2.*B*SMM(L,IJ1 ) ) GO TO 61 K = (NP(L»J)-1)*NB1 + NP(L,I) GO TO 75 K = (NP(L * I)-1)*NB1 + NP(L,J) IJl = I + J l S ( K ) '= S ( K ) + B*SMM(L, U I ) CONTINUE IF (NRE) 1 2 1 , 1 2 1 , 4 7 0 DO 4 6 9 I = 1,NRE K = ( ( N N 2 2 ( I ) -1)*NB) +1 SCK) = S ( K ) + E L A S T ( I ) CONTINUE SOLVE ^ 7 . •; . • , '7 1 o o 216 202 204 205 206 207 209 210 211 2 08 203 212 213 5 4 3 6 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 S(1)=SQRT (S(l)) DO 202 I =2 »NB S ( I ) = S ( I)/S<1 ) DO 203 L = L 1 , L 2 » N B K1=L-NBS K11=K1 K2=L-K3 IF(K1)204»204,205 Kl= ((L-l)/NB)+1 . DO 2 0 6 K=K1,K2,K3 S(LJ= S(L)-S(K)*S(K) IF(S(L))212,212,207 S ( L ) = SORT ( S ( L ) ) DO 208 N=l,NR2 M1=K11+(K3*N) I F ( M l ) 2 0 9 , 2 0 9 ,210 M1 = K1 LN=L+N DO 211 M=M1,K2,K3 MN=M+N S(LN)= S(LN)-(S(M)*S(MN)) S(LN)= S(LN)/S(L) S(LN+1)=S(LN+1)/S(L) GO TO 214 PRINT 2 1 3 , L FORMAT ( I X , 3 0 H S T R U C T U R E IS U N S T A B L E AT STOP i r— -o i DISP U) 214 228 217 219 218 220 223 221 224 225 2 22 226 227 536 Y 12 11 -10 9 3 7 6 . 5 537 472 538 "539" CONTI NUE CHECK EACH RECTANGULAR ELEMENT FOR S T A B I L I T Y DO 222 K= 1 » NRM XL AM = ASP ( K)#EX / ( l o - U X Y * U Y X ) XGAM = G / A S P ( K ) YLAM = E Y / ( A S P ( K ) * < l . - U X Y * U Y X ) ) YGAM = A S P ( K ) * G D I AA = ( XL AM + XGAM) *T ( K ) M . IF(A*T(K)-DIAA) 217,217,218 PRINT 2 1 9 , K FORMAT(IX,13HMEMBER NUMBERI5,32H IS U N S T A B L E DUE TO A TOO S M A L L ) D I A J = (YLAM + Y G A M ) * T ( K ) / 4 o L F _ i A J * T ( K ) ~ D I A J ) 2 2 0 , 2 2 0 ,221 : ; _ : PRINT 2 2 3 , K FORMAT(IX,13HMEMBER NUMBERI5,32H I S U N S T A B L E DUE TO J TOO S M A L L ) D I A A J = ( ( U X Y * X L A M / A S P ( K ) +ASP ( K ) * XGAM)* T ( K ) M « ) **2 ; AAJ = A * A J I F ( A A J # T ( K ) # T ( K ) - D I A A J ) 224,224,222 PRINT 2 2 5 , K \ FORMAT(IX,13HMEMBER NUMBERI5,33H I S U N S T A B L E DUE TO A J TOO SMALL) CONTINUE PRINT 2 26 ._. F O R M A T ( I X , 5 2 H A L L RECTANGULAR ELEMENTS NOT S T A T E D ABOVE ARE S T A B L E CONTINUE PRINT 536,NRS . FOR~MATI: 1H1,37HSTR ESS A N A L Y S I S FOR S T R U C T U R E NUMBER 14) DO 4 5 7 K A S E = 1,NLC DO 5 3 7 L=1,NU FK(L)=0.0 READ 4 7 2 , N J L FORMAT ( 1 1 0 ) _.. PRINT~538»KASE,NJL F O R M A T ( I X , 1 5 H L O A D CONDITION 14,6H HAS I 4 , 1 6 H LOADED J O I N T S . ) PRINT 539 .._ J O I N T NO. X FORCE Y FORCE) FORMAT(IX,32H _ ...p-_ . 03- DO 540 J = 1 » N J L • " _ READ 550 9 J N I U F F ( 1 ) »FF ( 2 ) FORMAT ( I 1 0 » 2 F . 1 0 o 3 ) FORMATt 1X9 I 1 Q 9 2 F 1 0 . 3 ) PRINT 5419 J N L U F F d ) » F F ( 2 ) DO 540 L = l 9 2 I = ND(JNU.L) ; . ' I F ( I ) 5429542,543 IF(FF(L))544,540*544 PRINT _ _ L i L _ J j _ J _ , — : F O R M A T ( I X 9 33H S T R E S S FROM FORCE I 2 . 1 1 H AT J O I N T 115H NOT I N C L U D E D . ) . • . GO TO 540 . __ F K ( I ) = FK( I ) + F F ( L ) CONTINUE BACK S U B S T I T U T E ONE ; ' J 2 = (NB*NU)+1 . J3=(NB+1)*NU DO 311 I = 1 » N U ; ; . KK= J 2 + I - 1 S(KK)=FK(I) S ( J 2 ) =S ( J 2 )VS( 1 ) ; • DO 302 L = L l 9 L 2 9 N B Jl=J2+((L-l)/NB) J = J1-K3 : : : K1=L-NBS IF(K1) 30493049303 K1 = 1+ ( ( L - 1 ) / N B ) J = J2 K2=L-K3 DO 301 K = K19_K2_9j<3__ S(Jl)=S(Jl)-S(K)*S(J) J=J+1 : 550 541 542 544 545 543 540 C 310 : 311 304 303 : 301 3_o_2 C _sj. J 3 J L f _ _ J i ) /s( L > ( BACK SUBSTITUTE TWO ' _ ; ; ; 14. 309 308 312 C 400 401 C T 408 12 JL 409 ••10 9 3 7 6 5 4 3 410 500 411 L3=L2+1 S(J3) =S(J3)/S(L2) DO 308 L R = L 1 , L 2 » N B L= L 3 - L R J l = J2+( ( L - D / N B ) J =J1+ 1 K1=L+1 K2=L+K3 DO 309 K=K1,K2 S(J1)=S(J1)-S(K)*S(J) J = J+1 S(J1)=S(Jl)/S(L) DO 312 1 = 1, Nil KK=J2+I-1 D(I)=S(KK) BACK S U B S T I T U T E THREE PRINT 400 FORMAT(1H1,17HJ0INT DEFLECTIONS) PRINT 401 F O R M A T ( I X , 5 0 H J O I N T NUMBER X-DEFLECTION JOINT DEFLECTIONS DO 501 K=1,NRJ N2 = N D ( K » 1 ) N3 = N D ( K , 2 ) IF(N2) 408,409,410 N4 = -N2 DX = - D ( N 4 ) 60 TO 500 DX = 0.0 60 TO 500 DX = D ( N 2 ) IF(N3) 411,412,413 N5 = -N3 DY = - D ( N 5 ) 60 TO 501 o 1 Y-DEFLECTION) 412 413 501 414 415 416 417 C 424 423 425 421 426 428 418 DY =0.0 GO TO 501 DY = D ( N 3 ) PRINT 414,K,DX,DY FORMAT(IX,I7,9X,1PE14.7,4X,1PE14.7) PRINT 4 1 5 FORMAT(1H1 ,14HMEMBER S T R A I N S ) PRINT 4 1 6 STRAIN FORMAT(IX,83HMEMBER SHEAR STRAIN 1 ' STRAIN STRAIN) PRINT 417 AT BOTTOM FORMAT(IX,83HNUMBER STRAIN AT TOP 1 AT RIGHT AT L E F T ) MEMBER S T R A I N S DO 418 K=l,NRM DO 421 1 = 1 ,8 IF(NP(K,I ) )423,424,425 DM(I) = 0.0 GO TO 421 NPI = - N P ( K , I ) DM(I) = - D ( N P I ) GO TO 421 NPI = N P ( K , I) DM(I) = D(NPI ) CONTINUE GAM(K) = 0 . 5 * ( ( D M ( 1 ) + D M ( 4 ) - D M ( 5 ) - D M ( 8 ) ) / A S P ( K ) + D M ( 2 + D M ( 7 ) - D M ( 3 ) 1-DM(6))/HORL(K) EXT(K) = ( D M ( D - DM(4))/HORL(K) EXB(K) =(DM(8) - D M ( 5 ) ) / H O R L ( K ) E Y R ( K ) -=(DM(2) - D M ( 7 ) ) / ( A S P ( K ) # H O R L ( K ) ) E Y L ( K ) =(DM(3) - D M ( 6 ) ) / ( A S P ( K ) * H O R L ( K ) ) PRINT 428,K,GAM(K) » E X T ( K ) » E X B ( K ) »EYR ( K ) » E YL ( K )• FORMAT(IX,14 ,4X,1PE14.7,1P4E16.7) GO TO 418 CONTINUE • i 1 i) ,3 i i I • 440 1 I 1 441 • 442 C t 494 493 495 492 - 446 T 1 12 447 ..__!! 2 ' iO 9 3 3 7 4 6 5 4 5 3 452 448 453 PRINT 4 4 0 FORMAT(1H1,12HN0DAL FORCES) PRINT 4 4 1 0FORMAT(1X»126HMEMBER Fl F2 1 F4 F5 . F6 2 F8) PRINT 442 FORMAT ( I X , 6HNUMBER) MEMBER FORCES DO 4 4 8 K = l , NRM DO 492 1=1,8 IF(NP(K,I)1493,494,495 • D M ( I ) = OoO GO TO 4 9 2 NPI = - N P ( K , I ) DM(I) = -D(NPI ) GO TO 4 9 2 NPI = NP(K , I ) DM(T) = D ( N P I ) CONTINUE DO 4 4 6 J = l , 8 N=((J*(17-J))/2)-8 DO 4 4 6 I = J , 8 IN = I+N S A ( I , J ) =• SMM( K, IN ) DO 4 4 7 1=1,8 DO 4 4 7 J = l , 8 SA(I,J) = SA(J,I) DO 4 5 2 1=1,8 F ( I ) = 0.0 DO 4 5 2 J = l , 8 F(I) = F ( I )+ SA(I•J)*DM(J) PRINT 45 3 , K , F ( 1 ) , F ( 2 ) , F ( 3 ) , F ( 4 ) , F ( 5 ) , F ( 6 ) » F ( 7 ) , F ( 8 ) F O R M A T ( I X , 1 4 , 1 P E 1 5 . 7 , 1 P 7 E l 6.7 ) PRINT 4 5 5 F3 F7 ' i N) 1 o 455 456 .1 458 466 467 457 f FORMAT ( 1H 1 » 27HMEMBER S T R E S S E S AT C E N T R O I D ) PRINT 4 5 6 F O R M A T ( 1 X » 6 1 H M E M B E R NUMBER X-STRESS Y-STRESS IR S T R E S S ) PAR = l . - U X Y * U Y X DO 4 5 7 K = 1 » N R M S X Y ( K ) = GAM(K)*G SX(K) = E X * ( E X T ( K ) + E X B ( K ) + U X Y * ( E Y L ( K ) + E Y R ( K ) ) )/(2o*PAR ) SY(K) = EY*(EYL(K)+EYR(K)+UYX*(EXT(K)+EXB(K)))/(2o*PAR) PRINT 4 6 7 » K * S X ( K ) » S Y ( K ) ? S X Y ( K ) FORMAT( 1 X » I 7 » l P E 2 1 o 7 » l P 2 E 1 7 o 7 ) CONTINUE GO TO 9 9 9 STOP END i 12, 1 11 10 ' 2 9. 3 8 . 7 4 6. 5 4 5 3. 6 SHEA , ; vn —
- 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 |
Aggregated Source Repository | 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050614/manifest