UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Laterally loaded wood compression members : finite element and reliability analysis Koka, Exaud Noe 1987

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1987_A7 K64.pdf [ 4.74MB ]
Metadata
JSON: 831-1.0062671.json
JSON-LD: 831-1.0062671-ld.json
RDF/XML (Pretty): 831-1.0062671-rdf.xml
RDF/JSON: 831-1.0062671-rdf.json
Turtle: 831-1.0062671-turtle.txt
N-Triples: 831-1.0062671-rdf-ntriples.txt
Original Record: 831-1.0062671-source.json
Full Text
831-1.0062671-fulltext.txt
Citation
831-1.0062671.ris

Full Text

LATERALLY LOADED WOOD COMPRESSION MEMBERS: FINITE ELEMENT AND RELIABILITY ANALYSIS by EXAUD NOE KOKA B.Sc.Eng.(Hons).University Of DSM , TANZANIA,1984 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of C i v i l Engineering We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1987 © EXAUD NOE KOKA , 1987 In presenting t h i s thesis in p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at The University of B r i t i s h Columbia, I agree that the Library s h a l l make i t freel y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of C i v i l Engineering The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: October 1987 ABSTRACT This thesis consists of two parts. The f i r s t part describes the analysis and implementation of a f i n i t e element computer model for the general prediction of f a i l u r e of wood members in bending or in combined bending and a x i a l compression. Both i n s t a b i l i t y and material strength f a i l u r e s are included. The program i s v e r i f i e d using available a n a l y t i c a l and test r e s u l t s . A good agreement with the results predicted by th i s program i s observed. The second part describes a procedure for the structural r e l i a b i l i t y evaluation of a compression member assuming random loads and material variables. The program developed here for the r e l i a b i l i t y study l i n k s the f i n i t e element program and the Rackwitz-Fiessler algorithm for the ca l c u l a t i o n of the r e l i a b i l i t y index /3. The gradient of the f a i l u r e function, which i s a necessary input to the Rackwitz-Fiessler algorithm, i s computed numerically using the f i n i t e element routine. The results of the r e l i a b i l i t y study for a t y p i c a l column problem are compared against the available results obtained by following the code procedures [as outlined in CAN3-086.1-M84 (1984)] for di f f e r e n t slenderness r a t i o s . A performance factor <p = 0.75, for compression members of any length i s recommended in order to obtain a more accurate and consistent l e v e l of r e l i a b i l i t y in the design i i process. It i s estimated that i f t h i s factor <p^ = 0.75 adopted in the current design practices, a l e v e l r e l i a b i l i t y index of the order of 4.0 can be achieved. TABLE OF CONTENTS ABSTRACT i i LIST OF FIGURES v i LIST OF TABLES v i i i NOTATION ix Acknowledgment xi 1. INTRODUCTION 1 1.1. Background 1 1.2. Objectives 3 1.3. Thesis Organisation 3 2. CURRENT CODE REQUIREMENTS AND PREVIOUS RESEARCH WORK .. 5 2.1. Introduction 5 2.2. Current Code Requirements 5 2.2.1. Concentric Compression 5 2.2.2. Combined bending and Compression 7 2.3. Previous Research on Wood compression members . 8 3. FINITE ELEMENT ANALYSIS 10 3.1. Introduction 10 3.2. Assumptions 10 3.3. Stress s t r a i n r e l a t i o n s h i p 11 3.4. F i n i t e Element Approximation 17 3.4.1. Introduction 17 3.4.2. Kinematic Assumptions 17 3.4.3. Problem formulation 20 3.4.4. Interpolation Functions 22 3.4.5. Strain Displacement Relations 26 3.4.6. V i r t u a l Work Equations 27 3.4.7. Newton-Raphson Method 28 3.4.8. Computation procedure 31 3.4.9. Numerical Integration 33 3.5. Convergence c r i t e r i o n for solution vector .... 34 3.6. Obtaining the ultimate load Pmax 35 3.7. Failur e C r i t e r i o n and Size effects 39 3.7.1. Size e f f e c t in compression 40 3.7.2. Size e f f e c t in tension 41 3.8. Program Structure 44 3.9. Discussion 44 4. PROGRAM VERIFICATION 46 4.1. Introduction 46 4.2. Comparison of results 46 4.2.1. Presentation of Results 47 i v 4.3. Discussion 53 5. APPLICATIONS 55 5.1. Introduction 55 5.2. Numerical example 56 5.3. Observations 58 6. RELIABILITY ANALYSIS 59 6.1. Introduction 59 6.2. The 0 method for R e l i a b i l i t y analysis 62 6.2.1. Rackwitz-Fiessler Algorithm 64 6.3. Problem Formulation 65 6.3.1. Code Design Equation 67 6.4. The G function for the problem 67 6.5. The r e a l i a b i l i t y Program 69 6.6. R e l i a b i l i t y Results 71 6.6.1. Discussion of results 75 7. Conclusions and Recommendations 76 7.1. Conclusions 76 7.2. Recommendations 77 REFERENCES 79 APPENDIX A 81 v LIST OF FIGURES 1. Axial load-slenderness relationship for concentric loading 7 2. Stress s t r a i n assumptions for wood 13 3. B i l i n e a r stress s t r a i n r e l a t i o n s h i p for wood 14 4. Di s t r i b u t i o n of stresses and strains 15 5. Stress-strain relationship for various m 17 6. Large deformation of a beam element 18 7. Large deformation of axis of beam element 19 8. i f ck f i n i t e element in the x-coordinate system 21 9. Flow chart for obtaining the solution vector {X} 33 10. Estimating the i n i t i a l f a i l u r e load Pi 36 11. Iteration process for obtaining Pmax 38 12. Stress p r o f i l e across the section 42 13. Comparison of program (with m=-l) and a n a l y t i c a l results [3] 48 14. Axial load-Slenderness plots for the data of Table 2,(no size e f f e c t s included) 51 15(a). Axial load-Slenderness curve with size effect taken into account, e = 2mm 51 15(b). Axial load-slenderness curves with varying k c for a constant kfc 52 16. Some loading cases and support conditions 55 17(a). Ultimate strength interaction curves for simply supported columns subjected to uniformly di s t r i b u t e d load 57 v i 17(b). Non-dimensionalized ultimate strength interaction curves of Figure 17(a) 58 18. Typical problem for r e l i a b i l i t y evaluation 59 19. P r o b a b i l i t y density function for the variable G 61 20. D e f i n i t i o n of the r e l i a b i l i t y index 0 63 21. R e l i a b i l i t y results as a function of <j>^', size effects not included 74 22. R e l i a b i l i t y results with size e f f e c t included in the r e l i a b i l i t y program 75 v i i LIST OF TABLES 1. Maximum deflections of a fixed ended uniformly loaded beam. (E = 10000 Mpa, 2x4-in section, L = 2m) 47 2. Axial load-slenderness data for a pin-ended 2x4-in beam (size effect neglected). 49 v i i i N O T A T I O N The f o l l o w i n g symbols are f r e q u e n t l y used i n t h i s r e p o r t P = a x i a l compressive l o a d NDS = N a t i o n a l Design S p e c i f i c a t i o n s (USA). a = normal s t r e s s A = c r o s s s e c t i o n a l area E Q = mean modulus of e l a s t i c i t y ( i n Mpa) m = slope of the s t r e s s s t r a i n curve, f a l l i n g branch e = s t r a i n u = a x i a l deformation w = t r a n s v e r s e deformation B = width of c r o s s s e c t i o n H = depth of c r o s s s e c t i o n L = Length of member £ , TJ = normalized c o o r d i n a t e s {6} = nodal displacement v e c t o r M, M,, M 2 = Shape f u n c t i o n s f o r the w displacements N, N, = Shape f u n c t i o n s f o r the u displacements K T = g l o b a l tangent s t i f f n e s s matrix {6 Q} = i n i t i a l displacement v e c t o r {A5} = incremental displacement v e c t o r AP = l o a d increment {X Q} = i n i t i a l s o l u t i o n v e c t o r {X} = New s o l u t i o n v e c t o r k t = s i z e e f f e c t f a c t o r i n t e n s i o n ix k c = size effect factor in compression G = f a i l u r e function Pf = pr o b a b i l i t y of f a i l u r e /3 = r e l i a b i l i t y index = factored compressive resistance p a r a l l e l to grain 0p = resistance factor in pure comression K c = Slenderness factor d = dead load variable / = l i v e load variable Design Load Factors a D = dead load factor a L = l i v e load factor 7! = r a t i o of nominal dead load to nominal l i v e load Subscripts c = compression t = tension PDN = n o m i - n a l dead load PLN = n o m ^ n a l l i v e load f = f a i l u r e x ACKNOWLEDGMENT I wish to express my deepest gratitude to my supervisor, Professor R. 0 . Foschi for his invaluable advice and guidance throughout the research work and in the preparation of t h i s thesis. Thanks are also due to Professor B. Madsen for his advice and support. The f i n a n c i a l support from the Commonwealth Scholarship and Fellowship Committee and some additional funding from the department of C i v i l Engineering of the University of B r i t i s h Columbia i s g r a t e f u l l y acknowledged. Lastly, thanks are due to those who have rendered help in t h i s project, d i r e c t l y or i n d i r e c t l y . xi 1. I N T R O D U C T I O N 1.1. B A C K G R O U N D Wood compression members subjected to l a t e r a l loads occur very frequently, such as in building frames, bridge or roof trusses and other important engineering structures. They are usually proportioned to s a t i s f y some l i m i t i n g stress c r i t e r i o n set by design s p e c i f i c a t i o n s or codes. The stresses developed at any cross section in such members consist of: 1. the a x i a l stress caused by the compressive forces , 2 . the primary bending stress due to the l a t e r a l loads, and 3. the secondary bending stress resulting from the amplification of the deflections produced by the compressive forces . The secondary bending stress becomes p a r t i c u l a r l y important for members with a high slenderness r a t i o and large compressive forces. The procedures for computing the secondary stresses in e l a s t i c columns are described in the l i t e r a t u r e on s t a b i l i t y theory [1] . Although e l a s t i c analysis i s used extensively in design computations, i t does not give an accurate indication of the true load-carrying capacity, p a r t i c u l a r l y for columns which are not very slender. L a t e r a l l y loaded columns generally f a i l by excessive bending after the stresses in some 1 2 portions of the member exceed a maximum value. To determine the ultimate strength of such columns, i t i s necessary to perform a s t a b i l i t y analysis that considers the el a s t o - p l a s t i c behaviour of the material. Most available design codes and s p e c i f i c a t i o n s use the t r a d i t i o n a l approach, which consists of assuming a linear e l a s t i c material with a maximum normal stress f a i l u r e c r i t e r i o n . Previous a n a l y t i c a l and experimental studies on wood, as reported in the l i t e r a t u r e [5,6,7], have shown that : 1. wood has a non-linear stress s t r a i n r e l a t i o n s h i p in compression, e.g.bilinear e l a s t o - p l a s t i c r e l a t i o n s h i p , and 2. this material c h a r a c t e r i s t i c contributes s i g n i f i c a n t l y to the behaviour of the column, p a r t i c u l a r l y at small slenderness r a t i o s . Furthermore,there are s t i l l some problems which remain unsolved: 1. The codes do not give guidance for c a l c u l a t i n g moments result i n g from beam-column deflections. 2. An account for p o s s i b i l i t i e s of du c t i l e y i e l d i n g in the compression zone or tension f a i l u r e in the tension zone i s not given. 3 1.2. OBJECTIVES This study i s aimed at achieving three main objectives, namely : 1 . To develop a f i n i t e element analysis for the general prediction of the f a i l u r e of a compression member under transverse loads. The analysis w i l l take into account the n o n - l i n e a r i t i e s due to slenderness e f f e c t s (geometric),a non-linear s t r e s s - s t r a i n r e l a t i o n s h i p for the material, and estimation of f a i l u r e load controlled by either tension or compression. 2 . The analysis w i l l be implemented in a computer program. The computer program w i l l allow f l e x i b i l i t y in accomodating various support conditions and load conf igurations. 3 . To evaluate the r e l i a b i l i t y of wood compression members assuming random loads and material variables. 1.3. THESIS ORGANISATION Part 2 provides a summary of current design code recommendations and previous research on wood compression members. Part 3 describes a general formulation of the f i n i t e element analysis and the computer implementation. Part 4 provides a v e r i f i c a t i o n of the computer program developed in Part 3 , using experimental results as reported by previous researchers [ 6 , 7 ] , Part 5 presents the 4 application of the analysis to a l a t e r a l l y loaded compression member, where a x i a l load versus transverse load interaction diagrams are developed for d i f f e r e n t slenderness ratios using a 2x4-in SPF cross section. Part 6 discusses the concept of r e l i a b i l i t y evaluation. Here, a computer program for the evaluation of the r e l i a b i l i t y index 0 of a wood compression member i s constructed, using the program developed in Part 3 and the Rackwitz-Fiessler algorithm. A summary of the results obtained from th i s study for a s p e c i f i c problem i s given at the end of the chapter. And l a s t l y , Part 7 provides a general conclusion of the report and some recommendations for further research and study. 2. CURRENT CODE REQUIREMENTS AND PREVIOUS RESEARCH WORK 2.1. INTRODUCTION The f a i l u r e c h a r a c t e r i s t i c s of a compression member depend on i t s slenderness. The ultimate capacity of short compression members i s d i r e c t l y related to the strength of the material in compression. With an increase in the length of the member, a change to a buckling type of f a i l u r e i s observed. Thus, a l a t e r a l i n s t a b i l i t y f a i l u r e i s ch a r a c t e r i s t i c of slender compression members. For a member of intermediate length, there i s a t r a n s i t i o n between these two types of f a i l u r e regimes, in which case the load capacity depends on both the compression strength and the s t i f f n e s s of the material. 2.2. CURRENT CODE REQUIREMENTS 2.2.1. C o n c e n t r i c Compression Current design codes c l a s s i f y compression members into short, intermediate or slender members according to their slenderness r a t i o C c. For rectangular cross-sections, 6 where L = le n g t h of the member d = dimension of the c r o s s - s e c t i o n of the member i n the d i r e c t i o n of b u c k l i n g . Thus, Short members are c o n s i d e r e d to be those with slenderness r a t i o s of 10 or l e s s . They w i l l normally f a i l by c r u s h i n g p a r a l l e l to the g r a i n . T h e i r design a l l o w a b l e l o a d i s based on the s p e c i f i e d s t r e n g t h i n compression p a r a l l e l to g r a i n , F , t h e i r c r o s s s e c t i o n a l area A, and a performance f a c t o r <f> , Cr P — - F 0 (2) A c P Slender members g e n e r a l l y f o l l o w the E u l e r b u c k l i n g r e l a t i o n and have slenderness r a t i o s exceeding C^, a number dependent on the mean e l a s t i c modulus E Q and the s p e c i f i e d compression s t r e n g t h , F , of the column m a t e r i a l . The number C. i s given by : where E c i s c a l l e d the modulus of e l a s t i c i t y f o r compression members and i s equal to 0.74E Q. For lumber, can vary between 20 to 25, depending on the grade of the member under considerat ion. Intermediate members have slenderness rat i o s between 10 and C^. They are designed using a modified compression strength which empirically interpolates between slenderness rati o s of 10 and C^. The above c l a s s i f i c a t i o n i s i l l u s t r a t e d - — Crushing Euler Figure 1. Axial load-slenderness re l a t i o n s h i p for concentric loading. 2 . 2 . 2 . Combined bending and Compression A compression member i s often subjected to bending about either one or both axes, and the combined effect of the bending and a x i a l loading must be considered. For this type of loading, most current codes specify a simple f a i l u r e c r i t e r i o n based on a linear interaction between the a x i a l 8 load capacity of a concentrically loaded column and the moment capacity in bending alone. Therefore, t h i s approach may be applicable as long as the wood member remains in the linear e l a s t i c range. Very l i t t l e has been done so far to predict the behaviour beyond the linear range. This may be attributed, in part, to the u n c e r t a i n i t i e s about the precise form of the c u r v i l i n e a r s t r e s s - s t r a i n r e l a t i o n s h i p of wood in compression. 2.3. PREVIOUS RESEARCH ON WOOD COMPRESSION MEMBERS Most previous studies on wood columns and beam-columns (Newlin and Trayer 1925; Wood 1950) have considered wood to be a linear e l a s t i c material which f a i l s when a l i m i t i n g compression stress i s reached. Thus, Larsen and Theilgaard (1979) tested wood members with combined a x i a l and transverse loads to v e r i f y their theory for beam-column behaviour. They used a second order l i n e a r d i f f e r e n t i a l equation to predict the deflections of e l a s t i c beams and beam-columns. Bleau (1983) and Buchanan (1984), conducted an extensive jo i n t experimental study on e c c e n t r i c a l l y loaded columns to ca l i b r a t e and v e r i f y their strength models. Their models are able to predict the strength of f u l l size lumber, using results of a x i a l tension and compression tests on similar members. Buchanan used a mean modulus of e l a s t i c i t y , E Q , 9 e q u a l t o 10000 Mpa t o c a l i b r a t e h i s model. Zann (1985), used B l e a u ' s d a t a (1983), w i t h E Q e q u a l t o 10400 Mpa t o c a l i b r a t e h i s s t r e n g t h model. Zann's model (1985), i s based on the NDS (1975) d e s i g n recommendations, and t a k e s account of b i a x i a l a c count of b i a x i a l b e n d i n g . A l t h o u g h i n both c a s e s good agreement w i t h the t e s t r e s u l t s i s r e p o r t e d , a q u e s t i o n which remains unanswered i s t h e f a c t t h a t i n each case a d i f f e r e n t i s used, and t h i s E„ i s not the one o ' o c o r r e s p o n d i n g t o the mean t e s t r e s u l t s . B l e a u (1983), r e p o r t s an E Q of 9660 Mpa f o r t h e p o p u l a t i o n t e s t e d . The model d e v e l o p e d i n t h i s r e p o r t i n c o r p o r a t e s some of the i d e a s d i s c u s s e d by the p r e v i o u s r e s e a r c h e r s , and p r o v i d e s a more g e n e r a l s o l u t i o n t o the beam-column problem. The method of f o r m u l a t i o n and t h e c o r r e s p o n d i n g computer i m p l e m e n t a t i o n w i l l be d i s c u s s e d i n P a r t 3 of t h i s r e p o r t . 3. FINITE ELEMENT ANALYSIS 3.1. INTRODUCTION This chapter describes the formulation of a f i n i t e element analysis for predicting the f a i l u r e of a wood member under d i r e c t a x i a l compression and l a t e r a l loads. The theory and assumptions in t h i s chapter w i l l be described along with the basis of a computer program developed to implement the model. 3.2. ASSUMPTIONS The following assumptions are made: 1. plane sections remain plane. 2. the s t r e s s - s t r a i n law for the material i s known. 3 . material properties are constant along the length of the member. 4. bending in only one plane i s considered. 5. no torsional or out of plane deformations are considered. 6. duration of load effects are not considered. 7. shear deformations are small, hence neglected. 10 11 3.3. STRESS STRAIN RELATIONSHIP Various studies [4,5,8] have focussed on the st r e s s - s t r a i n behaviour of wood in compression p a r a l l e l to grain, with the aim of deriving a mathematical relationship to represent th i s behaviour. Recently, Malhotra and Mazur (1970), proposed a stress s t r a i n r e l a t i o n of the form : 1 a e = — [ c a - ( l - c ) f l n d - — )] (4) E c f o c where: e = s t r a i n . a = stress f c = maximum compression stress E Q = mean modulus of e l a s t i c i t y c = shape parameter. For c = 0.99,. the curve described by Equation (4) i s shown as A in Figure 2. A mathematical equation for the stress s t r a i n curve for clear dry wood in compression at various grain angles was also developed by O 1Halloran(1973). The equation takes the following form a = E Qe - hen (5) Where a, e, and E are defined above and A,n are equation 1 2 constants determined by f i t t i n g the equation to a given set of experimental data. A plot of t h i s equation i s shown as curve B in Figure 2. This equation cannot be used beyond maximum s t r a i n because i t may take on negative values very rapidly. A comprehensive study on the s t r e s s - s t r a i n r e l a t i o n s h i p of timber with defects, in compression p a r a l l e l to grain, has been done by Glos (1978), as reported by Buchanan (1984). Based on experimental data, the curve shown as C in Figure 2 was proposed. This curve i s characterised by a number of material parameters that depend on measurable wood properties, namely density, moisture content, knot r a t i o and the percentage of compression wood. Using t h i s curve for modelling purposes necessarily involves the c a l i b r a t i o n of these parameters. A simple b i l i n e a r proposal by Bazan (1980), as discussed by Buchanan (1984), appears to be the most recent one. In th i s proposal, i t i s assumed that the slope of the f a l l i n g branch i s a variable which can be a r b i t r a r i l y taken as that value which produces maximum bending moment for any neutral axis depth. A plot of t h i s curve i s shown as D in Figure 2. 13 Figure 2. Stress str a i n assumptions for wood The analysis in t h i s study uses the simple b i l i n e a r curve D, with the exception that the slope of the f a l l i n g branch of the s t r e s s - s t r a i n r e l a t i o n i s considered to be a material property, in agreement with Buchanan (1984). The curves in Figure 2 are characterised by a li n e a r e l a s t i c and a non-linear part. Therefore the stresses can generally be expressed as a = E Qe + F(e) (6) 14 The s t r e s s - s t r a i n r e l a t i o n s h i p adopted in t h i s study i s as shown in Figure 3 and includes linear e l a s t i c behaviour in tension, with a b i l i n e a r r e l a t i o n s h i p in compression and a f a l l i n g branch after maximum stress. Figure 3. Bi l i n e a r stress s t r a i n relationship for wood Using the above s t r e s s - s t r a i n relationship, the resu l t i n g d i s t r i b u t i o n of stresses and strains in a rectangular beam i s as shown in Figure 4. T r section strains stresses Figure 4. Di s t r i b u t i o n of stresses and strains The curve in Figure 3 can be mathematically represented by the following expressions : For segment 1-2 ; o = ~\£c\ ~ |fc|m - mEce (7) For segment 2-3 ; ° = V (v Or, in combination, a = E Qe - [E Qe + |fc|(1+m) + mE Qe](l - A(e + |e c|)) (9) 16 where A(e+|e j) i s the step function defined as follows : i f e > - |e c| ; A = 1 (10) i f e < - |e c| ; A = 0 Hence, for the case of e l a s t o - p e r f e c t l y p l a s t i c behaviour, m = 0; and Equation (9) reduces to a = Ee - {Ee + |f |} (1-A(e + |e I)) (11) For the e l a s t i c case, m = -1, and we have Equation (9) both for tension and compression. This explanation i s further i l l u s t r a t e d in Figure 5 below. m < o m = - l Figure 5. Str e s s - s t r a i n r e l a t i o n s h i p for various m. 1 7 3.4. FINITE ELEMENT APPROXIMATION 3.4.1. I n t r o d u c t i o n The f i n i t e element method i s a very powerful and v e r s a t i l e technique presently available for the numerical solution of problems of the type considered here. The advantages of the method have been recognised and i t s applications extensively demonstrated p a r t i c u l a r l y in steel and concrete structures, and for some wood structures such as wood f l o o r s , wood diaphragms and trusses. However, the application of the method to wood beam-column analysis has not been explored to an equivalent degree. 3.4.2. K i n e m a t i c Assumptions As the displacements become large, a geometric non-linearity is introduced in the deformation of a beam-column. Consider a beam element undergoing large deformations but small st r a i n s . For the geometry shown in Figure 6, u and w are, respectively, a x i a l and l a t e r a l displacements of the centreline of the beam. A and 0 are two points on the same plane such that O i s on the beam centerline (axis) and A i s at a distance z from O (positive z ) . Line OA represents conditions before deformation, while l i n e O'A1 represents conditions after deformation. Assuming that plane sections remain plane, the rotation 18 dw of the cross-section i s 6 = — . Figure 7 shows two points, dx A and B, at the same distance z from the centreline. After deformation, these points are at A* and B'. ^ beam Figure 6. Large deformation of a beam element 19 Figure 7. Large deformation of axis of beam element. From the geometry of Figure 7 i t follows that du. dw. ds 2 = dx 2 [(1 + — & ) 2 + ( — A ) 2 ] dx dx If the expression above i s expanded binomially, and i f the higher order terms are neglected, the following s i m p l i f i e d expression i s obtained. du, 1 dwx d s = d x d + — & + - ( — £ ) 2 + ) (13) dx 2 dx Therefore, the corresponding s t r a i n e. at a distance 2 is 20 du, 1 dw, A j . _ / A \ 2 A dx 2 dx dw But, from Figure 6, u. = u - z — (15) A dx thus, du, du d2w — A - - - - z — - (16) dx dx dx' Also, neglecting higher order terms, wA = w (17) Thus combining Equations (14),(16) and (17), we get the st r a i n at a height z as : du 1 dw d2w e= — + - ( — ) 2 - z (18) dx 2 dx dx 2 3 . 4 . 3 . Problem formulation A beam element with two end nodes i s used in the formulation. Let us choose a l o c a l coordinate li (-1 < ^ < 1) in each element such as the one shown in Figure 8 below. Thus, along the x axis coordinate system each element has two end nodes, i and j separated by a length 2A. 21 x=o £ = - 1 £=0 u du dx w dw dx u du. dx w dw dx Figure 8. i * " * 1 f i n i t e element in the x-coordinate system. Thus, the x-coordinate of any point within the element can be expressed as x = x c + A£. The elemental nodal degree of freedom vector is represented as in Equation (19). There are 4 degrees of freedom at each node, namely the two dispacements u and w and their respective f i r s t d e rivatives. { 5 } = u i du ( — - ) • dx 1 dw ( — )• dx 1 U j du ( — )• dx 3 w . D dw ( — )• dx ^ (19) 22 3.4.4. Interpolation Functions In t h i s study, complete cubic interpolations are used to approximate the displacements u and w within an element. It is important to note that in order to s a t i s f y compatibility conditions, the displacement u only requires a linear interpolation. However, a cubic interpolation was used to give an improved approximation of the a x i a l stress with fewer elements. A complete cubic interpolation requires 4 parameters to define the function. The displacements and the f i r s t derivatives at the two nodes provide s u f f i c i e n t parameters to f u l l y describe a cubic polynomial function. The displacements u and w are thus given as follows: du u( *> - ( ~ " - U 2 + - £ 3 ) u . + - ( 1 - * - £ 2 + S 3 ) ( _ ) U V < ' ~ O A A 1 8 1 dx (20) 1 3 1 1 du + ( - + - £ + - £ 3)u.+ - (-1-£+£ 2+£ 3) ( — )• 2 4 4 ^ 8 dx ^ / - \ 1 3 1 1 dw W U ) = ( _ - _ £+ _ £ 3 ) W . + - ( I " ? " £ 2 + £ 3 ) ( _ ) 2 4 4 1 8 dx (21) 1 3 1 1 dw + ( _ + _ € + _ £ 3)w.+ _ (-1-$+£ 2+£ 3) ( ) 2 4 4 J 8 dx J 2x - 2x where £ = — 2A 23 In vector matrix notation, we can write U = {N} T{5} ; — = { N , } T {5} (22) dx dw d2w w = {M}T{6} ; — = {M,} T {5} ; = = {M 2} T{6} dx dx 2 where N , N , , M, M, and M 2, are vector functions of £ given by the following expressions: r 1 3 1 - - - £ + - £ 3 2 4 4 A - d - £ - £ 2 + £ 3 ) 4 0 . 0 {N} = 0 . 0 1 3 1 2 4 4 - ( - 1 . 0 - £ + £ 2 + £ 3 ) (23) 0 . 0 0 . 0 24 {N,} = 1 3 3 _ (- _ + A 4 4 - ( - 1 . 0 - 2 £ + 3 £ 2 ) 4 0 . 0 0 . 0 1 3 3 - ( - - - £ 2 ) A 4 4 1 ( - 1 . 0 + 2 £ + 3 $ 2 ) 0 . 0 0 . 0 (24) 0 . 0 0 . 0 {M} = 1 3 1 - - - £ + - £ 3 2 4 4 A - ( 1 . 0 - * - $ * + $ 3 ) (25) 0 . 0 0 . 0 1 3 1 - + - £- - £ 3 2 4 4 (-1 . 0 - £ + £ 2 + £ 3 ) 0 . 0 25 0 . 0 1 3 3 - ( - - + - i ; 2 ) A 4 4 - ( - 1 . 0 - 2 £ + 3 £ 2 ) 4 0 . 0 0 . 0 1 (!-!«•> A 4 4 1 ( - 1 . 0 + 2 £ + 3 £ 2 ) 4 0 . 0 0 . 0 (26) 1 3 — ( - O A 2 2 — (-2.0 + 6U 4A 0.0 0 . 0 1 3 — ( - - £ ) A 2 2 1 — ( 2 . 0 + 6 £ ) 4A (27) 26 3.4.5. Strain Displacement Relations For a l a t e r a l l y loaded column problem with large deformations, the w displacements w i l l be much larger than the a x i a l displacements u. Thus the s t r a i n displacement terms considered are similar to those derived in Equation (18) above, where: du d2w 1 dw e = — - z — + - ( — ) 2 (28) dx dx 2 2 dx Substituting the displacement functions into Equation (28), the results in symbolic form are e = [B+B(6)]6 (29) where B, represents the linear s t r a i n displacement terms, while B(6), which i s a function of the 8 parameters, includes the contribution of the non-linear s t r a i n displacement terms. Thus h B = N,1" - zM 2 T, where z = - T? (30) ~ 2 1 B(5) = - {8} T M,^!T 27 3.4.6. V i r t u a l Work Equations The system of equations governing the problem is obtained v i a the p r i n c i p l e of v i r t u a l work. Defining # as a v i r t u a l dislacement of the nodal variables, the resulting v i r t u a l s t r a i n s , 7, are given by 2 = [B + C(6)]S (31) 3 where C(5) = — {B(6)} i s a linear function of 8. Thus, i f 95 -we neglect i n e r t i a forces, the v i r t u a l work equation reduces to / v£odv = (3 2) where P i s the consistent load vector, calculated using the shape functions as indicated in Equations (23), (24), (25), (26) and (27). V i s the volume of the member. Substitution of the res u l t i n g equation for 7_ into the Equation (32) leads to the system of governing equations for an element, that i s , J v t B + C(6)] Tadv = P (33) Assembling the element equations in the usual f i n i t e element manner leads to the global system of equations. In order to 28 f i n d the solution to the nonlinear system of Equations (33) , i t i s convenient to introduce the vector function of 8, $( 8), such that *(8) = J v[B + C ( 8 ) ] T adv - P (34) The solution 8 now has to s a t i s f y *(8) = 0. The zeros of $(8) may be found numerically via the Newton-Raphson procedure as outlined below. 3.4 .7. Newton-Raphson Method This i s a commonly used technique to solve non-linear equations. The method uses a f i r s t order approximation technique to solve non-linear equations through i t e r a t i o n . Thus, at 8+A8, the f i r s t order approximation for the function $ w i l l be d$(8) $(8+AS) = $(5) + [ ]A8 (35) 8{8} where $(8) i s a function of the displacement vector {8}. The 29 above equation can further be s i m p l i f i e d into (36) (37) where [K T] represents the tangent s t i f f n e s s matrix. D i f f e r e n t i a t i n g the right hand side of Equation (34) by parts, we obtain a s i m p l i f i e d expression for K T. Equation (36) permits an i t e r a t i v e procedure to determine the vector {8} s t a r t i n g from an i n i t i a l approximation {8^}. Thus, in general, A6 = - KT~1<I»(8) d{$(8)} . = [ ] 1 <HS) d{8} 6 i + 1 " 5 i " [ K T ] * ( V (38) * ( 6 i ) = -MBo + B < 6 i > ] a i " p (39) where the matrix K T i s obtained as shown in Equation (40) , which follows. 30 + EJ* vB TC ( 5)dV E / v C T ( 5 ) B d V E / V C T ( 6 ) C ( 5)dV E/ vB TBdV E o(l . 0 + m ) / v ( l . 0 E0(1 . 0 - m ) / v ( l . 0 E0(1 . 0 - m ) / v ( l .0 E o ( 1 - 0 " m ) / v ( l .0 / VM,M,adV Me + | e c | ) ) B T BdV (40) A ( 6 + \ec\))BT C ( 5)dV + I e c | ) ) C T ( 6 ) B d V Me + | e c | ) ) C T ( 5 ) C ( 5)dV 31 3.4.8. Computation procedure The Cholesky decomposition of the matrix K T i s u t i l i z e d to determine the vector {A8}from Equation (38). The boundary conditions are f i r s t applied to the s t i f f n e s s matrix K T and to the vector {$(6)}. The boundary condition codes for t h i s program are as follows 1 = u du 2 = — dx 3 = w dw 4 = — dx Thus, to enforce a boundary condition equal to zero, zeros are placed into the off diagonal locations for the row and column corresponding to the s p e c i f i e d degree of freedom in [ K T ] , while a zero i s placed for the same degree of freedom in the returned load vector {$(5)}. In addition to t h i s , a value 1 i s placed into the diagonal term of the sp e c i f i e d degree of freedom in the [K T] matrix. Then the matrix [K T] i s decomposed and f i n a l l y a solution {A5} i s obtained. Each element of the vector {AS} i s compared against an acceptable tolerance sp e c i f i e d by the user to determine whether a need to do more i t e r a t i o n s is necessary in order to refine the solution. A summary of the whole procedure given in the flow chart in Figure 9. 32 is GEOMETRY /LOADING AND MATERIAL PROPERTIES INITIALISE ARRAYS T INCREMENT P * COUNT ITERATIONS ESTABLISH GLOBAL VALUES {DELTA R), CK] V SOLVE SYSTEM OF EQUATIONS (KI{DELTA X} = {DELTA R} {X> = {Xo> + {DELTA X} i CONVERGENCE [ ACHIEVED ± SOLUTION =» {X} Figure 9. Flow chart for obtaining the solution vector {X}. 33 3.4.9. Numerical Integration Since the volume integrals in the expression for K T are complicated, i t i s d i f f i c u l t to obtain closed form solutions. Hence numerical integration i s used. Gaussian quadrature scheme has been applied due to i t s s u i t a b i l i t y in the l o c a l coordinate system varying from -1 to +1. According to Zienkiewicz (1979), the maximum order of the polynomial appearing in the integral determines the number of Gaussian points necessary to accurately integrate the function. Thus, the term {B} in Equation (40) contains a fourth order polynomial in £ and at the same time {B} i s squared in the expression for [ K T ] . Therefore, the highest order polynomial term in the integrals i s of order 8. Knowing that a k point Gaussian scheme w i l l integrate exactly a (2k-l) order polynomial, i t follows that a 5-point Gaussian scheme i s needed in the numerical integration. Thus, the integrals over the volume V become 2BHA r\ 2BHA N N lv = ~A— / d ^ d r i " = - 7 - £ £ K T { i ' 3 ) w i w - i ( 4 1 ) v 4 J_x J_l 4 I=I j=j J where N = 5 was chosen, and w^  and w_j are the corresponding Gaussian weights. 34 3.5. CONVERGENCE CRITERION FOR SOLUTION VECTOR The convergence of the solution vector at every load step is checked by the Euclidean norm c r i t e r i o n . If we l e t 3 Q be the previous solution vector and 3 be the present solution, then, as shown below, le t Ax = |d - d Q| represent the difference between the lengths of d and d . We can then write l 3 o l 2 = X o ( i ) | 3 | 2 = X 2 ( i ) where X Q ( i ) and X(i) are the components of d Q and d, respect i v e l y . 35 Then, Ax = (x(i) - x o ( i ) ) 2 (42) The convergence c r i t e r i o n based on the Euclidean norm i s defined as 3.6. OBTAINING THE ULTIMATE LOAD PMAX The f a i l u r e load Pmax i s obtained by an i t e r a t i v e procedure. For fast convergence to the solution Pmax, the following approach for estimating an i n i t i a l guess for the f a i l u r e load i s chosen. F i r s t of a l l , the crushing strength P c as well as the Euler buckling load Per of the member are computed. Regardless of the support conditions and member length, the ultimate load w i l l be less than the smallest value between P„ and Per and w i l l l i e within the shaded c region of Figure 10. Ax | 3 0 | < sp e c i f i e d tolerance (43) 36 Figure 10. Estimating the i n i t i a l f a i l u r e load P i . The minimum of P and Per i s then taken to be the c i n i t i a l f a i l u r e load P^. As a f i r s t step, we l e t the solution l i e between two load values, namely P1=0 and P2=P^. The average load P3=(P1+P2)/2 becomes the f i r s t t r i a l load. The f i n i t e element solution i s obtained for P=P3. If f a i l u r e occurs, i t means that the solution i s between the values P=P1 and P=P3. Therefore we set the minimum and the maximum loads for next i t e r a t i o n as P1=P1 and P2=P3. A new P3 (P1+P2)/2 i s calculated and the f i n i t e element program re-run. If the member survives, i t means that the solution 37 is now between the values P=P3 and P=P2. Therefore, we set the minimum and the maximum loads for next i t e r a t i o n as P1=P3 and P2=P2. This process i s repeated several times u n t i l an acceptable tolerance i s reached between two successive estimates of Pmax. If t h i s tolerance i s defined as TOLP, the it e r a t i o n s are stopped when P2-P3 TOP = < TOLP P3 The process i s summarized in the flow chart of Figure 11, where TOLP i s the allowable tolerance normally set by the user. 38 Figure 11. Iteration process for obtaining Pmax. 39 3.7. FAILURE CRITERION AND SIZE EFFECTS Due to the b r i t t l e fracture phenomenon which i s commonly observed in wood members, i t may be important that the associated size effects be incorporated in the analysis. In a b r i t t l e material, a decrease in member strength i s normally observed as a result of a corresponding increase in member si z e . If no size e f f e c t s are considered, the f a i l u r e c r i t e r i o n i s "max " F t < 4 4 ) where c m a x i s the maximum te n s i l e stress in the member. This c r i t e r i o n , although simple, does not resu l t in di f f e r e n t strengths between pure tension and pure bending. Such differences are accountable through the incorporation of size e f f e c t s . Weibull's theory of b r i t t l e fracture w i l l be applied to incorporate the size effect phenomenon. Thus, for a member of volume V, f a i l u r e i s related to the int e g r a l I = J ok dv (45) 40 with a corresponding f a i l u r e c r i t e r i o n given by I = ( a * ) k (46) where a = stresses in the member, * a = strength of a unit volume under uniform stress, k = size e f f e c t factor, V = volume of the stresses domain. 3 .7.1. Size e f f e c t in compression t The parameter |f | is the f a i l u r e stress in pure compression (buckling restrained). This may be considered 41 subject to size e f f e c t s , according to k ^ k f c CV = (F c ) C (47) * (48) where * F c = f a i l u r e stress in pure compression for a unit volume, k c = size effect parameter in compression, V = t o t a l volume of the domain under compression that i s , the entire member. 3.7.2. Size e f f e c t in tension Let F,p be the strength in pure tension. Then, from Equations (44) and (45) we have at any p r o b a b i l i t y l e v e l : (49) where kfc i s the size effect factor in tension, V i s the t o t a l volume and V T i s the domain of the t e n s i l e stresses. In the context of the analysis presented here, consider a f i n i t e element i and the l o c a l ^-coordinate system, as shown in Figure 12. The stresses within the element are assumed to follow the stress s t r a i n relationship as indicated Figure 4. L e t us introduce a l o c a l coordinate 77 (0 ^ 77 ^  Osuch that = rjh. The t e n s i l e stresses w i l l be linear in y, or a = r\a where a„ i s the maximum stress at the edge y = h. J H ^ ' 5 < ) : i y l 2A 1 Figure 12. Stress p r o f i l e across the section. 43 Equation (49) can then be expressed as 2^BHA / d£ / drj a u = F T UV (50) where N i s the number of elements. Since a = Oy(£)r), f'hU) k f ' k 1 A ( £ ) k / a T t({) d£ / r? fcdTj = / a T t(£)d£ f j i ; • H H J0 k t + 1 y _ i H then, Equation (50) becomes 1 v^/hU) k r k r — ; r/ J * T U)d£ = Fm 2N(k.+l)4—'I H T T t i=l or, f i n a l l y , kfc (a Tmax) k t^U/h(£) a T<«) kfc  T 2(k f c+l)N 4-^ 1 H aTmax '-I The location of the neutral axis, h(£), where the stresses a change sign, can be obtained by interpolation of the stress f i e l d . The implementation of the procedure in the f i n i t e element computer program follows the equations as derived above. A summary of the steps follow below. 1 . c T(£) is determined at a l l points £ and for a l l elements; 2. obtain the largest of the 0 T ( £ ) , aTmax to normalize the 44 stresses. 3. obtain h(£) for any cross section. A section f u l l y in compression w i l l result in h(£) = 0. 4. Integrate over each element and add, according to Equation (52). 5. Compare the aTmax with the maximum stress possible according to the f a i l u r e c r i t e r i o n of Equation (52). 3.8. PROGRAM STRUCTURE The computer program consists of a number of subroutines which read the structure's geometry and load data, carry out numerical integration, decompose matrices, solves system of equations and checks the convergence of the solution vector. The program enables the user to analyize beams, columns or beam-columns of various configurations. A time subroutine has been provided to give the amount of computer time used to solve each s p e c i f i c problem. This time i s calculated in cpu seconds. A l i s t i n g of the program has been provided in Appendix A. 3.9. DISCUSSION The analysis developed here o f f e r s numerous p o s s i b i l i t i e s . The material behaviour law can be modified to study d i f f e r e n t materials or the effect of several parameters in a single material. Also the dimensions of a 45 member cross-section, the e c c e n t r i c i t y of a x i a l load, l a t e r a l l y acting loads and support conditions can be varied. In the following chapter, the model w i l l be v e r i f i e d by considering some problems for which there are available experimental or theoreti c a l r e s u l t s . The computer program developed here does not take into account torsional or out of plane deformations. Also creep effects were not included in the analysis. It i s also anticipated that there could be a s i g n i f i c a n t variation of modulus of e l a s t i c i t y E Q along the length of the member. However, without loss of generality, and in the presence of r e l i a b l e experimental data, the program can be e a s i l y modified to accommodate such variations in E . The o approximation for the s t r e s s - s t r a i n relationship used i s suitable for small and intermediate levels of s t r a i n , but obviously can not be extrapolated to very large s t r a i n s . 4. PROGRAM VERIFICATION 4.1. INTRODUCTION In t h i s chapter, the f i n i t e element computer program developed in the previous chapter i s v e r i f i e d with reference to 1. t h e o r e t i c a l results from the theory of e l a s t i c beam-columns, and 2. the results of an extensive experimental program on a large number of timber members in s t r u c t u r a l sizes [as reported by Bleau (1983) and Buchanan (1984)]. The test material was SPF lumber, purchased in 16ft. (4.88m) lengths as 'Number 2 and Better' grade in Quebec, Canada. The program i s v e r i f i e d using the mean test r e s u l t s , namely modulus of e l a s t i c i t y E Q = 9660 Mpa, compressive strength f = 32.3 Mpa and t e n s i l e strength f f c = 30.35 Mpa. 4.2. COMPARISON OF RESULTS The f i r s t comparison considers a n a l y t i c a l results [3] and the computer program's e l a s t i c predictions using m = -1, where m is the slope of the f a l l i n g branch of the s t r e s s - s t r a i n curve in compression. The second comparison presents plots and tables of a x i a l load versus slenderness r a t i o to compare the mean maximum load from tests with what the present analysis predicts for several end e c c e n t r i c i t i e s 46 e. No size effects are considered. The t h i r d comparison i s similar to the second one, except that in th i s case the size e f f e c t phenomenon is taken into consideration, and the effe c t of varying k for a chosen k. i s evaluated. 4.2.1. P r e s e n t a t i o n of R e s u l t s Table 1 shows a comparison of linear and non-linear t h e o r e t i c a l results [3] and computer predictions for a uniformly loaded fixed ended beam. The data of Table 1 i s plotted in Figure 13. Qo Wraax [m] [kn/ml iTlmoshenko] [program] I l i n e a r ] 00.000 06.885 13.771 20.656 27.541 34.426 41.312 48.197 55.082 0.00000 0.01291 0.02447 0.03516 0.04406 0.05162 0.05874 0.06497 0.07076 0.00000 0.01266 0.02437 0.03469 0.04371 0.05159 0.05859 0.06485 0.07053 0.00000 0.01285 0.02570 0.03855 0.05140 0.06426 0.07711 0.08996 0.10281 Table 1. Maximum deflections of a fixed ended uniformly loaded beam. (E = 10000 Mpa, 2x4-in section, L = 2m) 48 0.11 L= 2m E Q= 10000 Mpa Q Q-kN/m o ^ Linear, — Timoshenko , x program 2x4-in SPF 40 60 Figure 13. Comparison of program (with m=-1) and a n a l y t i c a l results [3]. Table 2 shows a comparison of f a i l u r e loads as obtained by the computer program to test results for d i f f e r e n t e c c e n t r i c i t i e s e. A graphical plot of the data in th i s table i s shown in Figure 14, with no size effects considered. Similar results with size effects included are shown in Figures 15(a) and 15(b). 49 COMPUTER RESULTS TEST RESULTS e = 2mm e = 39mm e = 2mm e = 39mm c c Pmax [Kn] PmaxtKn] 3.37 100.953 41.327 5.10 100.111 40.066 104.35 48.21 6.74 98.008 38.593 8.99 95.064 36 . 490 11.24 90.437 34.177 14.61 80.131 30.601 69.02 32.68 16.85 70.022 28.175 19.10 60.125 25.676 20.22 55.170 24.442 48.75 24.71 21.35 50.897 23.376 24.72 40.024 20.356 25.80 36.934 19.410 34.98 — 28.10 31.793 17.759 32.60 24.220 14.829 35.96 20.054 13.194 20.52 14 .11 40.45 15.974 11.258 e -it-Table 2. Axial load-slenderness data for a pin-ended 2x4-in beam (size effect neglected). Data Input : 2x4-in SPF section mean E Q , f c , f f c Figure 14. Axial load-Slenderness plots for the data of Table 2,(no size e f f e c t ) . 110 0 10 20 30 40 C c 51 110 10 H 1 | | 1 1 ; | 1 1 5 15 25 35 45 Figure 15(a). Axial load-Slenderness curve with size effect taken into account, e = 2mm. Input Data: mean E^, f and f. o c t k = 20.0 and k. = 5.0 c t Figure 15(b). Axial load-slenderness curves with varyi ng for a constant kfc. Tests, l - k = 5 2 - k c = 10, 3 - k c = 20 4- k = 100, 5 - k = 150 c c 53 4.3. DISCUSSION The r e s u l t s as presented i n t h i s c hapter, show that there i s a r e l a t i v e l y good agreement between the t e s t and the p r e d i c t i o n s by the computer program, the agreement being a very good one f o r members with a high s l e n d e r n e s s r a t i o . For compression members i n the intermediate range, the program p r e d i c t i o n s are s l i g h t l y higher than the t e s t r e s u l t s . For very short members the program p r e d i c t i o n s are s l i g h t l y below. S e v e r a l e x p l a n a t i o n s can be put forward to e x p l a i n these d i s c r e p a n c i e s . The most obvious one i s that the s t r e s s - s t r a i n curve used f o r t h i s study may not be a tr u e r e p r e s e n t a t i o n of the a c t u a l behaviour. N e v e r t h e l e s s , s i n c e t h i s f e a t u r e may be changed i n the a n a l y s i s , the f i n i t e element technique developed here remains a powerful and general t o o l to study the behaviour and design c o n s i d e r a t i o n s of timber columns and beam-columns. A p p l i c a t i o n of t h i s computer program to wood beam-columns w i l l be d i s c c u s s e d i n the f o l l o w i n g chapter. As shown in F i g u r e 15(a), when k = 20.0 and k. = 5.0 are taken as input i n t o the program, the r e s u l t s are s l i g h t l y improved with r e s p e c t t o the ones where no s i z e e f f e c t was c o n s i d e r e d . A l s o , i t i s noted that s i z e e f f e c t s i n compression have l i t t l e s i g n i f i c a n c e f o r very slender members, while these s i z e e f f e c t s p lay a major r o l e i n very short and intermediate members. The reason f o r t h i s i s that 54 for short members, the volume subjected to tension i s small or non-existent. For slender members, the f a i l u r e i s controlled by the modulus of e l a s t i c i t y and member i n s t a b i l i t y . When k c i s very large, the results obtained are the same as the ones in Table 2, meaning that there i s no size effect for large k c« It appears from Figure 15(b) that k = 20.0 gives a best f i t to the test r e s u l t s . 5. APPLICATIONS 5.1. INTRODUCTION The application of the program to solve wood beam-columns w i l l be discussed in th i s chapter. This program can handle multiple spans with d i f f e r e n t load and support configurations. Among them are the ones shown in figure 16. Q p P Q J . (a) (b) Q o Q p P - U I I I I I I I I l i \ i i i i i r r (c) (d) Figure16. Some loading cases and support conditions. The members are assumed to be prismatic. The desired responses can usually be represented as load versus centre deflection curve or any other convenient way for a 55 56 p a r t i c u l a r l a t e r a l load Qo. Once the complete curves are obtained, the maximum loads can be ea s i l y determined from the peak of the curves. The computer program developed in th i s study provides an easier approach to the above process in that one gets the maximum load d i r e c t l y by supplying the program with the appropriate information. In a l l cases of Figure 16, the l a t e r a l loads Q or Qo cause bending moments about the major axis of the cross-section. It i s further assumed that weak axis buckling and l a t e r a l - tor s i o n a l buckling are e f f e c t i v e l y prevented so that f a i l u r e i s always caused by excessive bending in the plane of the applied l a t e r a l load . In performing the numerical procedure, i t i s assumed that the l a t e r a l load Qo i s applied f i r s t and maintained at a constant value as the a x i a l compressive load P increases or decreases. 5.2. NUMERICAL EXAMPLE As a numerical example, case (c) in Figure 16 has been considered in thi s study, using a 2x4-in SPF section and mean values for E Q , f and f f c . Also kfc = 5.0, m = 0.02 and k c = 10.0 has been used in obtaining the P versus Qo results as shown in Figures 17(a) and 17(b). Figure 17(a). Ultimate strength interaction curves for simply supported columns subjected to uniformly distributed load. Qo/Qou 0 0 2 0.4- 0 4 0.8 1 Figure 17(b). Non-dimensionalized ultimate strength interaction curves of Figure 17(a). 58 5 . 3 . OBSERVATIONS From Figures 17(a) and 17(b), i t can be noticed that P-Qo relationships predicted by the computer model are not a linear one as i t i s normally assumed in the current design practice for d i f f e r e n t slenderness r a t i o s . Additional research i s needed here in order to come up with a si m p l i f i e d design procedure for wood beam-columns. 6. RELIABILITY ANALYSIS 6.1. INTRODUCTION This chapter describes the procedure for the str u c t u r a l r e l i a b i l i t y analysis of a wood compression member. The problem to be studied is as shown in Figure 18; where P represents the applied a x i a l compressive load (for only dead and l i v e loads). L represents the length of the member while H and B represents the height and breadth of the cross section. yfS~r I Q u ^ L , B , H S P F Figure 18. Typical problem for r e l i a b i l i t y evaluation. The r e l i a b i l i t y of a member simply means the prob a b i l i t y that i t w i l l perform as intended in a prescribed s i t u a t i o n . It i s influenced by the demands on the structure and the capacity of the structure to respond to those demands. In general, one can define a performance or f a i l u r e function G to characterize the state of the structure in rel a t i o n to some performance c r i t e r i o n . This function G can be expressed 59 60 as G = C - D where C = structural capacity D = demands on the structure. The function G as defined above is posi t i v e whenever the capacity exceeds the demand, therefore the structure meets the performance c r i t e r i o n . On the other hand, the function G w i l l be negative whenever the demands exceed the capacity, resulting in the structure not meeting the required performance. When the function G i s exactly equal to zero, the structure i s on the threshold between meeting and f a i l i n g to meet the performance c r i t e r i o n , and such a state i s defined as " l i m i t state". e The p r o b a b i l i t y of f a i l u r e p^ i s the compliment of the r e l i a b i l i t y . Thus Pj = 1.0 - r e l i a b i l i t y According to the d e f i n i t i o n of G above, the pr o b a b i l i t y of f a i l u r e i s then given as 61 p f = Probability ( G < 0) Each design problem w i l l contain a set of intervening variables, and depending on the nature of the problem, some of the variables may be random, obeying some d i s t r i b u t i o n function. Thus, i f some of the basic variables are random, i t i s obvious that G w i l l be i t s e l f a random variable. The pro b a b i l i t y d i s t r i b u t i o n for G could be derived from a knowledge of the ind i v i d u a l p r o b a b i l i t y d i s t r i b u t i o n s for the basic variables, and the result would be as shown in Figure 19. G = 0 G G Figure 19. Probability density function for the variable G . The p r o b a b i l i t y of f a i l u r e p^ w i l l be the area under the curve to the l e f t of the o r i g i n G = 0. If this p r o b a b i l i t y 62 of f a i l u r e exceeds some desired value, one or more of the design variables would be changed and p^ recalculated u n t i l i t meets the required target. The p r o b a b i l i t y d i s t r i b u t i o n for G could be obtained by a n a l y t i c a l means using multiple integrations and the joint p r o b a b i l i t y d i s t r i b u t i o n s between the basic variables. This i s a very tedious and d i f f i c u l t approach. MonteCarlo simulation can be used to obtain the pr o b a b i l i t y of f a i l u r e in an approximate manner. In t h i s approach the value of G i s computed for a large number of combinations of the basic variables and p^ i s estimated from the proportion of times the G was negative. The selection of values for the basic variables must obey their joint p r obability d i s t r i b u t i o n s , and when more than two variables are involved, the procedure becomes d i f f i c u l t , tedious and expensive. In the following section, an approximate and fast procedure for estimating p^ w i l l be discussed. 6.2. THE 0 METHOD FOR RELIABILITY ANALYSIS In order to estimate the p r o b a b i l i t y of f a i l u r e p^ with s u f f i c i e n t accuracy but without resorting to complicated integrations or computer simulations, Hasofer and Lind [1974] introduced the concept of r e l i a b i l i t y index 0 using geometric approach. Thus, for a design problem containing N uncorrelated random variables X-, i = 1,..,N, with mean X. 63 and standard deviation a^, a set of "normalised " variables x^ is introduced. These variables have zero mean and standard deviation equal to 1 .0 , and are given as X . - X . x. = 3, ( 5 3) °i The f a i l u r e function can now be expressed in terms of the new, normalised variables x^ as shown schematicaly in the figure below, in which the horizontal plane represents the space of the variables x and the v e r t i c a l axis the function G. C(x) Figure 2 0 . D e f i n i t i o n of the r e l i a b i l i t y index j3. Hasofer and Lind showed that the r e l i a b i l i t y index 0 can be intepreted as the minimum distance between the o r i g i n 0 64 and the l i m i t state G Q. This i s a geometric problem which can be solved by successive i t e r a t i o n s using, for example, Hasofer and Lind's proposed algorithm. Knowing /3, the pr o b a b i l i t y of f a i l u r e i s obtained from p f = # ( - 0 ) (54) where $ i s the standardised normal p r o b a b i l i t y function. For p£ to be exact, we require that a l l the basic variables be normally d i s t r i b u t e d and G be linear in the basic variables. Figure 20 shows the case when the mean point belongs to the "safe domain" G > 0. The combinations of x^ which correspond to G = 0 (the l i m i t state) are represented by the curve G Q. 6.2.1. Rackwitz-Fiessler Algorithm This i s in actual fact the modification of the Hasofer and Lind Algorithm in order to improve the estimation of the pro b a b i l i t y of f a i l u r e . The modification refers to the case when the basic variables are non-normal. Rackwitz and Fi e s s l e r [1978], suggested a transformation of the o r i g i n a l random variables X^ into a set of normalised uncorrelated standard variables z^ using the following transformation z. = #~ 1[F(X.)] (55) 65 where $ i s the standard normal p r o b a b i l i t y d i s t r i b u t i o n function and F(X^) i s the cumulative d i s t r i b u t i o n function for the variable X^. The standard algorithm from Hasofer and Lind i s then used for the new variable z^. This modification improves the prediction of p^ because i t meets one of the two conditions mentioned e a r l i e r namely, that a l l variables be normally d i s t r i b u t e d . This algorithm i s presently the accepted norm for the evaluation of the r e l i a b i l i t y index 0. 6 . 3 . P R O B L E M F O R M U L A T I O N To i l l u s t r a t e the a p p l i c a b i l i t y of the theory discussed above to normal practice, l e t us consider the column problem of Figure 18. The cross sectional dimensions of the column are B for width and H for depth. The length of the column i s represented as L. It i s assumed simply supported under an a x i a l compressive load P (for both dead and l i v e loads). The demand on the structure i s the applied load P. Thus, D = P = P D + P L where D = demand P D = dead load P L = l i v e load If 66 where Then d = - ° - (56) P DN / = (57) PLN PLN = n o m i n a l (design) l i v e load PDN = n o m i n a l (design) dead load D = P L N [ 7 i d + /] (58) where d and / are considered to be random variables. The P factor 7, i s a constant defined as 7, = -SH or the r a t i o of PLN nominal dead load to nominal l i v e load. The capacity C i s the maximum load, Pmax, the member can carry; thus C = Pmax = P{E 0,fc,ft,B,H,L,m] (59) and the f a i l u r e function can be expressed as G = C - D G = Pmax - P (60) where f = strength in compression f f c = strength in tension. m = slope of st r e s s - s t r a i n curve in compression. The problem can now be studied using the Rackwitz-Fiessler algorithm and the f i n i t e element computer program developed 67 in part 3 of t h i s thesis. However, i t i s convenient for the purpose of future code development to modify Equation (60) above to bring in the design equation format adopted for the code. 6 . 3 . 1 . Code Design Equation For members subjected to pure a x i a l compression, the Canadian Code, CAN3-086.1-M84 (1984) sp e c i f i e s the following desing equation. A D P D N + A L P L N * *p A F c K C ( 6 1 ) where #p = performance factor in compression. A = cross sectional area of member. K = slenderness factor c F c = F i f t h percentile compression strength (a D,a L) = load factors (1.25 and 1.5 respectively). 6 . 4 . THE G FUNCTION FOR THE PROBLEM Considering the l i m i t i n g case of Equation (61), we obtain the following equation: 68 P L N [ a D ? i + a L ] = *p A F c K c ( 6 2 ) where <p AFcKc P L N = -B (63) L N V + a L Combining Equations (60), (62) and (63),we can express the f a i l u r e function as 0 AFcKc G = Pmax - -E [ 7 , d + / ] a D7,+a L or 0 A Fc Kc G = P{E ,f f ,B,H,L,m} - -E [ 7 l d + /] (64) a D7,+a L For the purpose of t h i s study, the following variables have been considered random modulus of e l a s t i c i t y E Q compressive strength f te n s i l e strength f f c dead load variable d l i v e load variable / and the following have been considered to be constants with mean average values height of cross section H 69 breadth of cross section B length of member L slope m. 6.5. T H E R E A L I A B I L I T Y P R O G R A M The computer program which implements the derivation above i s attached in Appendix A. As part of the input, the program requests the number of random variables (in thi s case 5 ) , the type of their d i s t r i b u t i o n (according to a d i s t r i b u t i o n code), and the relevant parameter information to characterize the d i s t r i b u t i o n s . The program can accept the following d i s t r i b u t i o n s Code D i s t r i b u t i o n 1 Normal 2 Lognormal 3 Weibull A Gumbel 5 Ranked Data The fixed parameters 7 , , <j>^, a D and a L are provided by the user for each p a r t i c u l a r problem. The subroutine GXPR computes the function G and i t s gradient by c a l l i n g the f i n i t e element subprogram. The GXPR routine returns the value of G and the gradient vector DELTA. For the column problem discussed in thi s thesis the elements of the gradient vector corresponding to the f i r s t 3 random 70 variables were obtained numerically, while the remaining two were obtained by d i f f e r e n t i a t i n g the f a i l u r e function e x p l i c i t l y . Thus, the t o t a l elements of the gradient vector considering only f i v e random variables are obtained as G(E *) - G(E ") Delta(l) = 2 2 — 2AE^ o G(f +) - G(f ~) Delta(2) = S $ — 2Af c G(f. +) - G ( f / ) DeltaO) = S £ — 2Af t <t> AFcKc Delta ( 4 ) = — E 7 l a D 7 i + a L <(> AFcKc Delta (5) = -a D 7 i + a L The fixed parameters are passed onto the routines GXPR and COLUMN through a COMMON block. 71 6.6. R E L I A B I L I T Y R E S U L T S Keeping the r a t i o 7 l = 1.0, m = 0.02 and using a 2x4-in SPF section, the factor <p^ was changed and the corresponding r e l i a b i l i t y index 0 was computed for columns of d i f f e r e n t slenderness r a t i o s . Figure 21 shows the results for the r e l i a b i l i t y index 0 as a function of the performance factor 0p for the case of no size e f f e c t considered in the program. Figure 22 shows r e l i a b i l i t y results with size e f f e c t s included in the computer program. In obtaining the results for the two cases studied, the following information has been used for the random variables. (1) E Q : 3-parameter Weibull. Scale = 6738.0 Mpa Location = 3514.0 Mpa Shape = 3.97 Mean = 9660.0 Mpa (2) f : 3-parameter Weibull Scale = 33.845 Mpa Location = 0.0 Shape = 7.8559 Mpa F i f t h percentile = 15.87 Mpa (3) f. : 3-parameter Weibull Scale = 29.861 Mpa Location = 4.03 Mpa Shape = 2.911 Mpa (4) d = dead load variable : Normal Mean = 1.0 Standard deviation = 0.15 (5) / = l i v e load variable : normal Mean = 0.75 Standard deviation = 0.15 F i g u r e 22. R e l i a b i l i t y r e s u l t s wi th s i z e e f f e c t i n c l u d e d i n the r e l i a b i l i t y program. 75 6.6.1. Discussion of results From the results of Figure 21 i t is noted that for a performance factor = 0.75, a r e l i a b i l i t y index 0 of the order j3 = 4.0 i s achieved for a l l slenderness ra t i o s considered. Figure 22 shows a small increase of the r e l i a b i l i t y index for the same performance factor <p . cr Considering the two cases, a performance factor 0 = 0.75 Cr can be taken as a reasonable value to be included in the current design practices for columns of any length. The procedure outlined in th i s chapter for the r e l i a b i l i t y analysis of columns does not take into account the duration of load e f f e c t over the length of the servive l i f e of the column. A r e l i a b i l i t y study for beams taking into account the duration of load effect i s currently been done in the Department of C i v i l Engineering of the University of B r i t i s h Columbia. It w i l l be of interest for further research, to integrate the model developed here to th i s study. 7. CONCLUSIONS AND RECOMMENDATIONS 7.1. CONCLUSIONS From the results of part one of t h i s study, i t i s seen that the f i n i t e element analysis, including large deformations and non-linear material properties, can model wood column behaviour s a t i s f a c t o r i l y . The model does require accurate and r e l i a b l e input information on modulus of e l a s t i c i t y , compressive and t e n s i l e strengths. The results including size e f f e c t s show that for the size e f f e c t parameters k c = 20.0 and kfc = 5.0, the computer predictions for the maximum load Pmax agree f a i r l y well with test r e s u l t s . However, k c does have significance influence for short and intermediate columns, and should be known with some accuracy. For the r e l i a b i l i t y r esults in part two of this study, i t i s observed that the current performance factors <l> , as hr given in CAN3-086.1-M84 (1984), are more conservative than what th i s model predicts. A value of 0p = 0.75 appears to be a reasonable one for a l l slenderness r a t i o s . It i s estimated that i f t h i s new value of <f>^ i s adopted in the current design practice, i t w i l l give r i s e to a r e l i a b i l i t y index /3 of the order of 4.0. If a lower 0 i s required, a d i f f e r e n t #p should be introduced for short and intermediate columns. This points to a deficiency in the "column formula" giving 76 77 the slenderness adjustment factor K c. Idealy, t h i s factor should r e f l e c t the changes due to slenderness in such a way that the same <j>^ - 0 r e l a t i o n s h i p be obtained for a l l column lengths. 7.2. RECOMMENDATIONS It i s recommended that a performance factor <f>^ = 0.75 be used in the current design practice for a l l slenderness r a t i o s . However, prior to adopting this recommendation, there i s need to do more research in thi s area. In p a r t i c u l a r , the research should cover duration of load e f f e c t s , and the case of correlated variables; neither of which has been included in the analysis. The application of the Rackwitz-Fiessler algorithm requires a l l the variables involved to be uncorrelated. However, in some p r a c t i c a l cases some or more of the intervening variables w i l l be correlated. For example, in the context of the problems discussed in thi s report, the strength of beams, columns or beam-columns under combined a x i a l and l a t e r a l loads w i l l depend on the modulus of e l a s t i c i t y E o, the compression strength f c and the te n s i l e strength f f c. For lumber, these variables are p a r t i a l l y correlated and thi s must be dealt with, using for example the procedures available in the l i t e r a t u r e [12], before using the Rackwitz-Fiesler algorithm. 78 There i s not enough data available at present on size e f fects in both tension and compression, hence further p r a c t i c a l as well as theoret i c a l study i s necessary in order to come up with a r e a l i s t i c design recommendation applicable to lumber of a l l grades and species. REFERENCES 1. Chen,W.F., and Astuta,T., 1976, Theory of Beam Columns,  Vol. II, Space Behaviour and Design, McGraw-Hill Book Company, New York, 732p. 2. Zienkiewicz,0., 1971, The F i n i t e Element method in  Engineering Science. McGraw-Hill Book Company, London, England. 3. Timoshenko S., and Woinowsky-Krieger S., 1959, Theory of  plates and s h e l l s . McGraw-Hill Book Company, New York., pp. 396-428. 4. Malhotra, S.K., and Mazur, S.J., 1970, Buckling Strength  of s o l i d timber Columns, Engrg.J. 13(A-4), I-VII. 5. Larsen,H.J., and Thielgaard, E. (1979)., L a t e r a l l y loaded  timber columns , J. Struct. Div., ASCE, 105(7) 1347-1363. 6. Bleau, R., (1983).,Etude de Resistance des Colonnes en  bois de Qualite Commerciale, "Thesis presented to the University of Shebrooke at Quebec, in p a r t i a l fulfilment of the requirements for the degree of Doctor of Philosophy". 7. Buchanan, A. H., Strength Model and design methods for 79 80 bending and a x i a l load interaction in timber members., "Thesis presented to the University of B r i t i s h Columbia, at Vancouver, B.C., in p a r t i a l fulfilment of the requirements for the degree of Doctor of Philosophy". 8. Zann, J . J . , (1985)., Strength of lumber under combined  bending and compression., "USDA Forest Service Research paper FPL 391., Washington, D.C". 9. Hasofer, A.M., and Lind, N.C., 1974., "Exact and Invariant Second Moment Coe Format", Journal of Engineering  Mechanics D i v i s i o n , ASCE, Vol. 100, No. EM1, pp. 111-121. 10. Rackwitz, R., and F i e s s l e r , B. 1978. "Structural R e l i a b i l i t y Under Combined Load Sequences", Computers and  Structures , Vol. 9, pp. 489-494. 11. Foschi, R.O., 1979., "A discussion on the application of the safety index concept to wood structures", Canadian  Journal of C i v i l Engineering, Vol. 6, No. 1, pp. 51-58. 12. Der Kiureghian, A., 1986., "Structural R e l i a b i l i t y Under Incomplete Probability Information", Journal of Engineering  Mechanics , ASCE, Vol. 112, No. 1, pp. 85-104. 81 A P P E N D I X A 1. Program COLUMN.FOR (source code) 8 2 2 . Sample Input/Output F i l e s 96 3. Program RBETA.FOR (source code) 99 4. Sample Input/Output F i l e 119 PROGRAM COLUMN.FOR 8 3 c c c c c c c c c c c c c c c c c c c c c c c c c c c c c ******************************************** c c c c c c c c c c c c c c c c c c c c * * * * * * * * * * * * * * * * * * * * * * * * * * COLUMN.FOR Ve r s i o n 2.0 2 August, 1987 A PROGRAM FOR THE CALCULATION OF THE ULTIMATE COLUMN (OR BEAM-COLUMN) LOAD ON A MATERIAL BEHAVIOUR IS ELASTIC IN TENSION WITH BRITTLE FRACTURE, AND ELASTIC IN COMPRESSION UP TO A LIMITING COMPRESSION STRESS, WITH A FALLING LINEAR BRANCH BEYOND THAT LIMIT. 1 SIZE EFFECTS ARE CONSIDERED BOTH IN TENSION AND COMPRESSION. END LOAD IS A COMPRESSION LOAD. END LOAD CAN BE APPLIED ECCENTRICALLY. LATERAL LOADS CAN BE DISTRIBUTED OR CONCENTRATED. THE PROGRAM FINDS THE ULTIMATE END LOAD CORRESPONDING TO A GIVEN END ECCENTRICITY AND GIVEN LATERAL LOADS. THE PROGRAM CAN ALSO FIND THE ULTIMATE LATERAL LOAD WHEN THE END LOAD IS SPECIFIED TO BE ZERO ( NP = 0 ). PROBLEM DATA IS READ FROM UNIT #1. OUTPUT IS STORED IN UNIT #2. PROGRAM WRITTEN BY E. KOKA AND R.O. FOSCHI, UBC. * * * * * * * * * * * * * * * * * * * * * * * * * * ************************************************************ IMPLICIT REAL*8(A-H,0-Y) CHARACTER*20 NAMED 1 ,NAMEA1 ,NANS DIMENSION IX(21,4),F(8),NBC(21),TKO(672),XE(8) 1, R(84),XO(84),X(84),B(84),B1(8,8),B2(8,8),B3(8,8),B4(8,8) 2, B5(8,8),B6(8,8),B7(8,8),B8(8,8),B9(8,8),Y(5),RE(8),XP(84) 3, Q(20),IQ(20), ESTR(7), FI(7) COMMON/C1/GAP(5),GAW{5),EN 1(8,5),EM1(8,5),EM2(8,5), NGAUSS COMMON/C2/DIFP, NINT COMMON/C3/DEFL,PDEFL ********************************************************** * DEFINE VARIABLES * ********************************************************** NELEM = NO OF ELEMENTS NJBC = NO OF JOINTS WITH B.C. NBC(I) = NO OF B.C. AT NODE I IX = B.C. CODE 1 = U 2 UX 3 W 4 WX EN 1 = INTERPOLATION FUNCTIONS FOR u EMI,EM2 = INTERPOLATION FUNCTIONS FOR w GAP = CORDINATE AT GAUSS POINT GAW = CORRESPONDING WEIGHT NGAUSS = NO OF GAUSS POINTS NITER = MAX. NO OF ITERATIONS TOP = TOLERANCE FOR LOAD STEP EPSLON = TOLERANCE FOR SOLUTION VECTOR FC = MATERIAL STRENGTH IN COMPRESSION * * * * * * * * * * * * * * * * * 84 c * FT MATERIAL STRENGTH IN TENSION * c * EO MOE OF THE MATRIAL (RANDOM) * c * EN = SLOPE OF THE STRESS-STRAIN CURVE IN COMPR. * c * SPAN = MEMBER LENGTH * c * W = WIDTH OF SECTION * c * H = DEPTH OF SECTION * c * E = ECCENTRICITY OF AXIAL LOAD * c * NEQ = NO OF EQUATIONS TO BE SOLVED * c * NJOINT = NO OF NODES * c * NDOF = NO OF VARIABLES PER NODE * c * NODEL = NO OF NODES PER ELEMENT * c * NDIMB = NO OF VARIABLES PER NODE * c * LBW,LHB = HALF BANDWIDTH INCLUD. THE DIAG. * c * NA = NO OF UNKNOWNS FOR TOTAL PROBLEM * c * RE = ELEMENTAL LOAD VECTOR * c * R = STRUCTURE LOAD VECTOR * c * B GLOBAL LOAD VECTOR RETURNED * c * TKO = GLOBAL TANGENT MATRIX * c * XE = ELEMENT DISPLACEMENT VECTOR * c * X = GLOBAL SOLUTION VECTOR * c * B1,..B9 = ARRAYS FOR TEMPORARY STORAGE * c ********************************************** ************ WRITE( * , 6 ) 6 FORMAT(' ENTER DATA FILE NAME '/) READ(* ,8) NAMED1 WRITE( *,7) 7 FORMAT(' ENTER OUTPUT FILE NAME '/) READ(* ,8) NAMEA1 8 FORMAT(A) OPEN (1, FILE = NAMED1,STATUS='OLD') OPEN (2,FILE = NAMEA1,STATUS='NEW) READ(1 ,*> NELEM, NGAUSS NDOF = 4 NJOINT = NELEM+1 NG1 = NGAUSS + 1 NG2 = NGAUSS + 2 READ (1,*) NP,NQ,Q0 IF (NQ.EQ.O) GO TO 44 DO 4 3 I = 1, NQ 43 READ (1,*) I Q ( I ) , Q (l) 44 E = 0.0D0 IF (NP.NE.O) READ(1,*) E DO 65 I = 1, NJOINT 65 NBC(I)=0 READ ( 1,*) NJBC DO 75 I = 1, NJBC READ (1,*) N, NBC(N) READ (1,*) (IX'(N,J) , J=1 ,NBC(N) ) 7 5 CONTINUE NEQ = NDOF*NJOINT NODEL = 2 C C * READS MATERIAL STRENGTH IN COMPRESSION (FC) AND TENSION (FT), C BOTH CORRESPONDING TO THE SPECIFIED CROSS-SECTION AND THE C REFERENCE SPAN SREF. XKC AND XKT ARE THE WEIBULL SIZE EFFECT C SHAPE PARAMETERS IN COMPRESSION AND TENSION RESPECTIVELY. C READ(1,*) FC,FT,SREF,XKC, XKT NDIMB = NODEL*NDOF LBW = NDIMB LHB = LBW NA = LBW*NEQ C READ MOE AND SLOPE m OF CURVE READ(1,*) EO,EN C READ PROBLEM SIZE L, B, H READ(I,*) SPAN,W,H AR = W*H XI = W*H**3/12.D0 DEL = SPAN/(2.D0*NELEM) SLAMDA = SPAN/H C * ADJUST STRENGTHS TO THE ACTUAL VOLUME FC = FC *(SREF/SPAN)**(1.0/XKC) FT = FT *(SREF/SPAN)**(1.0/XKT) C OBTAIN SHAPE FUNCTIONS AND DERIVATIVES : N1,M1,M2 CALL SHAPE(DEL) WRITE(*,79) 79 FORMAT(' TOLERANCE FOR PMAX ? '/) READ(*,*) TOP WRITE(*,790) 790 FORMAT(' TOLERANCE FOR CONVERGENCE? '/) READ(*,*) EPSLON WRITE(*,791) 791 FORMAT(' MAX. NUMBER OF ITERATIONS? '/) READ(*,*) NITER WRITE(*,799) 799 FORMATC WANT TO SEE INTERMEDIATE RESULTS? (Y/N)'/) READ(*,8) NANS NINT = 0 IF (NANS.EQ.'Y'.OR.NANS.EQ.'y') NINT = 1 IF (NP.EQ.O) GO TO 761 PC = AR*FC PCR = 3.14159D0**2*E0*XI/(SPAN**2) PI = PC IF(PCR .LE. PC) PI=PCR P2 = PI P1 = 0.0D0 P3 = (P1 + P2)/2.0D0 NFAIL = 0 SMAX1 = 0.0 GO TO 760 761 FQ1 = 0.0D0 FQ2 = 1.ODO FQ3 = FQ2 NFLAG = 0 760 DO 792 J = 1, NEQ 792 XP(J) = 0.0D0 C C START CALCULATIONS FOR TRIAL LOAD LEVELS CALL TIME(ZIM) ' ZIM0 = ZIM 3773 CONTINUE P = 0.OD0 FQ = 1.0D0 IF (NP.NE.0) P = P3 IF (NP.EQ.O) FQ = FQ3 IF (NINT.EQ.1.AND.NP.NE.0) WRITE(*,4000) P 86 IF (NINT.EQ.1.AND.NP.EQ.O) WRITE(*,4001) FQ 4000 FORMAT(//' SOLUTION FOR P =',E15.6,* :'/) 4001 FORMAT(//' SOLUTION FOR LATERAL LOAD FACTOR=',E15.6,*:'/) C INITIALISE ARRAYS DO 80 J = 1, NEQ XO(J) = XP(J) 80 R(J) = 0.DO C EXTERNAL LOAD VECTOR R IF (Q0.EQ.0.0D0) GO TO 87 DO 81 J = 1 , 8 RE(J)=0.D0 81 CONTINUE RE(3) = FQ*Q0*DEL RE(4) = FQ*Q0*DEL**2/3.D0 RE(7) = RE(3) RE(8) = -RE(4) DO 83 NE = 1, NELEM DO 82 J J = 1, 8 K = (NE-1)*NDOF + J J R(K) = R(K) + RE(JJ) 82 CONTINUE 83 CONTINUE 87 IF (NQ.EQ.0) GO TO 185 DO 180 J = 1, NQ JS = (IQ(J)-I)*NDOF + 3 180 R(JS) = R(JS) + Q(J)*FQ 185 EM = P*E J J = (NJOINT-1)*NDOF + 1 R(JJ) = R(JJ)-P R(1) = R(1) + P R(4)=R(4)-EM R(NEQ)=R(NEQ)+EM ITER=0 C C BEGIN ITERATIONS AT THE TRIAL LOAD LEVEL 777 CONTINUE DO 84 1 = 1 , NA 84 TKO(I) = O.0D0 DO 85 K = 1, NEQ 85 B(K) = -R(K) DO 645 IE = 1, NELEM C INITIALIZE ARRAYS DO 88 I = 1 , 8 F(I) = 0.0D0 DO 86 J = 1 , I B1(I,J) = 0.0D0 B2(I,J) = 0.0D0 B3(I,J) = 0.0D0 B4(I,J) = 0.0D0 B5(I,J) = 0.0D0 B6(I,J) = 0.0D0 B7(I,J) = 0.0D0 B8(I,J) = 0.ODO B9(I,J) = 0.ODO 86 CONTINUE 88 CONTINUE C PICK ELEMENT SOLUTION FROM GLOBAL VECTOR DO 90 J J = 1 , 8 87 K = (IE - 1)*NDOF + J J XE(JJ) = XO(K) .90 CONTINUE DO 101 K = 1, NGAUSS Y(K) = 0.D0 DO 91 1=1, 8 Y(K) = Y(K)+XE(I)*EM1(I,K) 91 CONTINUE C OBTAINING COMPONENTS OF EKT DO 93 I = 1, 8 DO 93 J = 1, I B1(I,J) = Bl(I,J)+E0*DEL*EN1(I,K)*Y(K)*AR* 1 EM1(J,K)*GAW(K) B2(I,J) = B2(I,J)+E0*DEL*EM1(I,K)*Y(K)*AR* 1 EN1(J,K)*GAW(K) B3(I,J) = B3(I ,J)+E0*DEL*EM1(I,K)*Y(K)*AR* 1 Y(K)*EM1(J,K)*GAW(K) B4(I,J) = B4(I,J)+(E0*AR*DEL*EN1(I,K)*EN1(J,K)+ 1 E0*XI*DEL*EM2(I,K)*EM2(J,K))*GAW(K) 93 CONTINUE DO 100 L = 1, NGAUSS C STRESSES AND STRAINS AT GAUSS POINT STR = 0.5D0*Y(K)**2 DO 96 MO = 1 , 8 STR = STR+(EN1(MO,K)-GAP(L)*H*0.5D0*EM2(MO,K))*XE(MO) 96 CONTINUE STRE = STR+FC/EO FAC = 1.0D0 IF(STRE.GE.O.DO) FAC=0.0D0 STRESS = E0*STR-((E0+EN*E0)*STR+FC*(1.DO+EN))*FAC DO 99 I = 1, 8 DO 98 J = 1 , I B5(I,J) = B5(I,J)+DEL*0.5D0*AR*(EN1(I,K)-GAP(L) * 1 H*0.5D0*EM2(I,K))*(E0+E0*EN)*FAC*(EN 1(J,K)-H*0.5D0* 2 GAP(L)*EM2(J,K))*GAW(K)*GAW(L) B6(I,J) = B6(I,J)+DEL*0.5D0*AR*(EN 1(I,K)-GAP(L)* 1 H*0.5D0*EM2(I,K))*(E0+EN*E0)*FAC*Y(K)*EM1(J,K)* 2 GAW(K)*GAW(L) B7(I,J) = B7(I,J)+DEL*0.5D0*EM1(I,K)*Y(K)*AR* 1 (E0+E0*EN)*FAC*(EN1(J,K)~H*0.5D0*GAP(L)*EM2(J,K))* 2 GAW(K)*GAW(L) B8(I,J) = B8(I,J)+DEL*0.5D0*EM1(I,K)*Y(K)*AR* 1 (E0+E0*EN)*FAC*Y(K)*EM1(J,K)*GAW(K)*GAW(L) B9(I,J) = B9(I,J)+AR*STRESS*EMI(I,K)*EM1(J,K)* 1 GAW(K)*GAW(L)*DEL*0.5D0 98 CONTINUE F(I) = F(I)+AR*DEL*0.5D0*STRESS*((EN1(I,K)-H*0.5D0* 1 GAP(L)*EM2(I,K))+Y(K)*EM1(I,K))*GAW(K)*GAW(L) 99 CONTINUE 100 CONTINUE 101 CONTINUE C OBTAIN ELEMENT TANGENT MATRIX C EKT IS THE (I,J) COMPONENT OF THE ELEMENT TANGENT MATRIX DO 105 I = 1, 8 II = (IE-1)*NDOF + I B(II) = B(II) + F(I ) . DO 102 J = 1 , I J J = (IE-1)*NDOF + J 8 8 EKT = BI(I,J)+B2(I,J)+B3(l,J)+B4(I,J)~ 1 B5(I,J)-B6(I,J)-B7(I,J)-B8(I,J)+B9(I,J) I J = (JJ-1)*(LBW-1) + I I TKO(IJ) = TK0(IJ)+EKT 102 CONTINUE 105 CONTINUE 645 CONTINUE C INTRODUCE BOUNDARY CONDITIONS DO 1 1 1 UO = 1 , NJOINT IF (NBC(IJO).EQ.O) GO TO 111 DO 110 J = 1, NBC(IJO) II = (UO -1)*ND0F + IX(IJO,J) LBW1 = LBW - 1 DO 108 K = 1, LBW1 J J = II - LBW + K IF (JJ.LE.O) GO TO 1080 I J = (JJ-1)*(LBW-1) + II TKO(IJ) = 0.0D0 1080 J J = II + K IF (JJ.GT.NEQ) GO TO 108 I J = (II-1)*(LBW-1) + J J TKO(IJ) = 0.0D0 108 CONTINUE I J = (II - 1)*(LBW-1) + II TKO(IJ) = 1.0D0 B(II) = 0.0D0 110 CONTINUE 111 CONTINUE C C SOLUTION OF THE SYSTEM C CALL DECOMP(NEQ,LBW,TKO,IERROR) IF(IERROR .EQ. 1) GO TO 377 4 CALL SOLVN(NEQ,LBW,TKO,B) DO 1121 = 1, NEQ X(I) = XO(I)-B(I) 112 CONTINUE CALL CONVRG(XO,X(IER,NEQ,EPSLON,ITER) ITER = ITER + 1 IF (ITER.EQ.NITER) GO TO 431 IF (IER.EQ.2) GO TO 430 IF(IER.EQ.0) GO TO 118; DO 115 I = 1, NE@ 115 XO(I) = X ( I ) GO TO 777 430 I ERROR = 1 GO TO 3774 431 WRITE(2,900) NITER, P 900 FORMATC NO CONVERGENCE IN',13,'ITERATIONS AT P=',E13.6/) GO TO 901 C C AFTER CONVERGENCE, OBTAIN STRESSES AND STRAINS AT C THE CURRENT LOAD LEVEL C 118 CONTINUE EMAXP = O.0DO EMAXN = 0.0D0 SUME = 0.0D0 89 DO 550 IE = 1, NELEM DO 500 J = 1, 8 K = (IE-1)*NDOF +J XE(J) = X(K) 500 CONTINUE DO 540 K = 1, NGAUSS FACTOR = 0.0 DO 501 I = 1, 8 501 FACTOR = FACTOR + XE(I)*EM1(I,K) EPLUS = 0.5D0 * FACTOR**2 EMINUS = EPLUS DO 505 I = 1 , 8 EPLUS = EPLUS + (EN 1(I,R)-H*0.5D0*EM2(I,K))*XE(I) EMINUS = EMINUS + (EN1(I,K)+H*0.5D0*EM2(I,K))*XE(I) 505 CONTINUE IF(EPLUS.GT.0.0D0 .AND. EMINUS.GT.0.0) GO TO 506 IF(EPLUS.GT.0.0D0 .AND. EMINUS.LE.0.0) GO TO 507 IF(EPLUS.LE.0.0D0 .AND. EMINUS.LE.0.0) GO TO 508 IF(EPLUS.LE.0.0D0 .AND. EMINUS.GT.0.0) GO TO 509 506 EPOS = EPLUS IF(EMINUS.GT.EPOS) EPOS=EMINUS ENEG = 0.0D0 GO TO 530 507 EPOS = EPLUS ENEG = EMINUS GO TO 510 508 EPOS = 0.0D0 ENEG = EPLUS IF (DABS(EMINUS).GT.DABS(ENEG)) ENEG = EMINUS GO TO 530 509 EPOS = EMINUS ENEG = EPLUS C C * FINDS THE POSITION OF THE NEUTRAL AXIS C 510 ESTR(1) = EMINUS F l ( 1 ) = -1.0D0 ESTR(NG2) = EPLUS Fl(NG2) = 1.0D0 DO 512 L = 1, NGAUSS SUM = 0.5*FACTOR**2 DO 511 I = 1,8 511 SUM = SUM + (EN1(I,K) - GAP(L)*H/2.0*EM2(I,K))*XE(I) ESTR(L+1) = SUM Fl(L+1) = GAP(L) 512 CONTINUE DO 515 I = 1, NG1 PROD = ESTR(I)*ESTR(I+1) IF (PROD.LE.0.0D0) GO TO 516 515 CONTINUE 516 XN = F I ( I ) - ESTR(I)*(FI(I+1)-FI(I))/(ESTR(I+1)-ESTR(I)) IF (ESTR(I).EQ.O.0DO) GO TO 518 IF (ESTR(I).LT.0.0D0) HN = (1.0D0 - XN)*H/2.0D0 IF (ESTR(I).GT.0.0D0) HN = (1.0D0 + XN)*H/2.0D0 GO TO 520 518 IP (ESTR(I+1).LT.0.0D0) HN = (1.0D0 + XN)*H/2.0D0 IF (ESTRd + 1 ) .GT.0.0D0) HN = ( 1 . 0D0 - XN)*H/2.0D0 520 SUME = SUME + (HN/H)*(E0*EPOS)**XKT*GAW(K) 90 530 IF(EPOS.LT.EMAXP) GO TO 538 EMAXP = EPOS 538 IF(DABS(ENEG).LT.DABS(EMAXN)) GO TO 540 EMAXN = ENEG 540 CONTINUE 550 CONTINUE SMAXP = E0*EMAXP SMAXN = E0*EMAXN IF (DABS(SMAXN).LE.FC) GO TO 560 SMAXN = SMAXN -((E0 + EN*E0)*EMAXN + FC*(1.0 + EN)) 560 IF (SUME.EQ.0.ODO.OR.SMAXP.EQ.0.ODO) GO TO 563 SUME = SUME/(2.0*NELEM*(XKT+1.0)*SMAXP**XKT) FTT = FT * SUME**(-1.0D0/XKT) GO TO 564 563 FTT = FT 564 IF (SMAXP.GE.FTT) GO TO 3774 DEFL = 0.0D0 DO 565 IE = 1, NELEM J = (IE-1)*NDOF + 3 IF (DABS(X(J)).GT.DABS(DEFL)) DEFL = X(J) 565 CONTINUE J = NEQ - 1 IF (DABS(X(J)).GT.DABS(DEFL)) DEFL = X(J) IF (NP.EQ.O) PDEFL = FQ3 IF (NP.NE.0) PDEFL = P3 377 4 CONTINUE IF (NINT.EQ.0) GO TO 8810 IF (IERROR.EQ.1) WRITE(*,8888) IF (IERROR.EQ.0.AND.SMAXP.LT.FTT) WRITE(*,8889) SMAXP IF (IERROR.EQ.0.AND.SMAXP.GE.FTT) WRITE(*,8890) SMAXP 8888 FORMAT(' IERROR=1,FAILS (DIVERGENCE OR SINGULAR MATRIX)'/) 8889 FORMAT(' IERROR=0 SMAXP = ',E15.6,' SURVIVES'/) 8890 FORMAT(' IERROR=0 SMAXP = ',E15.6,' FAILS'/) 8810 CONTINUE IF (NP.EQ.O) GO TO 4500 IF (IERROR.EQ.1) GO TO 7330 IF (SMAXP.GT.FTT) GO TO 7331 IF (SMAXP.EQ.FTT) GO TO 7337 PI = P3 IF (SUME.EQ.0.ODO.OR.SMAXP.EQ.0.ODO) GO TO 5650 SMAX1 = SMAXP*SUME**(1.0D0/XKT) GO TO 5655 5650 SMAX1 = SMAXP 5655 DO 833 J = 1, NEQ 833 XP(J) = X(J) GO TO 8334 7330 P2 = P3 GO TO 8334 7331 P2 = P3 NFAIL = 1 SMAX2 = SMAXP*SUME**(1.0D0/XKT) 8334 IF (P1.EQ.0.ODO) GO TO 8338 TOLP = (P2-P1)/P1 IF (TOLP.LE.TOP) GO TO 7338 GO TO 8336 8338 IF (P2.LE.0.1D0) GO TO 7338 8336 IF (NFAIL.EQ.1) GO TO 8340 P3 = (PI + P2)/2.0 91 GO TO 3773 8340 P3 = PI + (P2-P1)*(FT-SMAX1)/(SMAX2-SMAX1) GO TO 3773 7337 P = P3 PP = P3 PAV = P3 GO TO 7339 7338 IF (P1.EQ.0.ODO) P2 = 0.0D0 P = P2 PP = PI PAV = (P+PP)/2.0 7339 CALL TIME(ZIM) ZIM = ZIM - ZIMO WRITE(2,570) PP,P,PAV,SMAXP,SMAXN,DEFL,PDEFL,SLAMDA 570 FORMAT(2X,' FAILURE BETWEEN LOADS ',E15.6,2X,' AND', 12X,E15.6/* AVERAGE=',E15.6/' EDGE STRESS (+) =',E15.6/ 2' EDGE STRESS (-) =',E15.6/' MAX. DEFLECTION =',E15.6, 3 ' AT LOAD =',E15.6/' SLENDERNESS = ',F6.2/) WRITE(*,683) ZIM 683 FORMAT(' TIME =',F7.1,' SECS.'/) WRITE(*,570) PP,P,PAV,SMAXP,DEFL,PDEFL,SLAMDA GO TO 901 4500 IF (IERROR.EQ.1) GO TO 4330 IF (SMAXP.GT.FTT) GO TO 4330 IF (SMAXP.EQ.FTT) GO TO 4337 IF (NFLAG.EQ.1) GO TO 4331 FQ1 = FQ2 FQ2 = 2.0D0*FQ2 GO TO 4580 4331 FQ1 = FQ3 4580 DO 4833 J = 1,NEQ 4833 XP(J) = X(J) GO TO 4334 4330 NFLAG = 1 FQ2 = FQ3 4334 IF (FQ1.EQ.0.ODO) GO TO 5338 TOLP = (FQ2-FQ1)/FQ1 IF (TOLP.LE.TOP) GO TO 4338 5338 IF (NFLAG.EQ.0) FQ3 = FQ2 IF (NFLAG.EQ.1) FQ3 = (FQ1+FQ2)/2.ODO GO TO 3773 4337 P = FQ3 PP = FQ3 PAV = FQ3 GO TO 4339 4338 P = FQ2 PP = FQ1 PAV = (P+PP)/2.0 4339 CALL TIME(ZIM) ZIM = ZIM-ZIM0 WRI TE ( 2 , 67 0 ) PP , P , PAV, SMAXP, SMAXN ,'DEFL , PDEFL 670 FORMAT(2X,' FAILURE BETWEEN LOAD FACTORS ',E15.6,2X,' AND', 12X,E15.6/' AVERAGE =',E15.6/' EDGE STRESS (+) =',E15.6/ 2' EDGE STRESS (-) =',E15.6/' MAX. DEFLECTION = *,E15.6, 3' AT LOAD FACTOR =',E15.6/) WRITE(*,683) ZIM WRITE(*,6 7 0) PP,P,PAV,SMAXP,SMAXN,DEFL,PDEFL 901 CONTINUE 3131 C C c c* c 4 c 5 C 10 350 150 9 2 CONTINUE CLOSE (1,STATUS='KEEP') CLOSE (2,STATUS='KEEP') STOP END END OF MAIN PROGRAM SUBROUTINE SHAPE(DEL) THIS SUBROUTINE CALCULATES DERIVATIVES OF SHAPE FUNCTIONS IMPLICIT REAL*8(A-H,0-Y) COMMON/C1/GAP(5),GAW(5),EN1(8,5),EM1(8,5),EM2(8,5),NGAUSS IF (NGAUSS.EQ.5) GO TO 5 IF (NGAUSS.EQ.4) GO TO 4 *** 3 POINT GAUSSIAN INTEGRATION -0.774596669241483D0 0.0D0 -GAP(1) 0.555555555555556D0 0.888888888888889D0 GAW(1) GAP( 1 GAP (2 GAP (3 GAW( 1 GAW(2 GAW(3 GO TO *** 4 GAP( 1 GAP (2 GAP (3 GAP (4 GAW( 1 GAW(2 GAW(3 GAW(4 GO TO *** 5 GAP( 1 GAP (2 GAP (3 GAP (4 GAP (5 GAW GAW GAW GAW GAW 10 POINT = -0 GAUSSIAN INTEGRATION 861136311594053DO -0.339981043584856DO -GAP(2) -GAP(1) 0.347854845137454D0 0.652145154862546D0 GAW(2) GAW(1) 10 POINT = -0 1 ) 2) 3) 4) 5) INITIALISES GAUSSIAN INTEGRATION 906179845938664D0 -0.538469310105683DO 0.0D0 -GAP(2) -GAP(1) = 0.236926885056189D0 = 0.478628670499366D0 = 0.568888888888889D0 = GAW(2) = GAW(1) EN1,EM1,EM2 DO 150 IL = 1 DO 350 IK = 1 EN 1 (IL , IK) • EM1(IL,IK) • EM2(IL,IK) : CONTINUE CONTINUE DO 250 I = 1 , EN 1(1 ,1 ) = EN 1(2 ,1 ) = EN 1(5 ,1 ) = EN 1(6 ,1 ) = EMI (3 ,1) = EMI (4 ,1) = 8 . NGAUSS 0.0D0 0.0D0 0.0D0 NGAUSS (-0.75D0+0.75D0*GAP(I)**2)/DEL (-1.D0-2.D0*GAP(I)+3.D0*GAP(I)**2)*0.25D0 (0.75D0-0.75D0*GAP(I)**2)/DEL (-1.D0+2.D0*GAP(I)+3.D0*GAP(l)**2)*0.2 5D0 (-0.7 5DO+0.7 5DO*GAP(I)**2)/DEL (-1.D0-2.D0*GAP(I)+3.D0*GAP(l)**2)*0.25D0 93 EM1(7,1 EM1(8,I EM2(3,1 EM2(4,1 EM2(7,1 EM2(8,I 250 CONTINUE RETURN END (0.75D0-0.75D0*GAP(I)**2)/DEL' (-1.D0+2.D0*GAP(I)+3.D0*GAP(I)**2)*0.25D0 1.5D0*GAP(I)/(DEL**2) (-2.D0+6.D0*GAP(I))/(4.D0*DEL) -1.5D0*GAP(I)/(DEL**2) (2.D0+6.D0*GAP(I))/(4.D0*DEL)' SUBROUTINE DECOMP(NN,LHB,AA,IERROR) C * THIS SUBROUTINE DECOMPOSES A MATRIX USING CHOLESKY C METHOD FOR BANDED,SYMMETRIC,POS. DEFN. MATRICES IMPLICIT REAL*8(A-H,0-Y) DIMENSION AA(672) C TKO IS STORED COLUMNWISE. IERROR = 0 KB = LHB-1 C DECOMPOSITION IF(AA(1).LE.0.D0) IERROR=1 IF(IERROR.EQ.1) RETURN AA(1) = DSQRT(AA(1)) IF(NN.EQ.I) RETURN DO 551 I =2, LHB 551 AA(I) = AA(I)/AA(1) DO 590 J = 2, NN J1 = J-1 IJD = LHB*J-KB SUM = AA(IJD) KO = 1 IF(J.GT.LHB) KO=J-KB DO 555 K = KO, J l JK = KB*K+J-KB 555 SUM = SUM-AA(JK)*AA(JK) IF(SUM.LE.O.DO) IERROR=1 IF(I ERROR.EQ.1) RETURN AA(IJD) = DSQRT(SUM) DO 568 I = 1, KB II = J+I KO = 1 IF (II.GT.LHB) KO=II-KB " SUM = AA(IJD+I) IF(I.EQ.KB) GO TO 565 DO 540 K = KO, J1 JK = KB*K+J-KB IK = KB*K+II-KB 540 SUM = SUM-AA(IK)*AA(JK) 565 AA(IJD+I) = SUM/AA(IJD) 568 CONTINUE 590 CONTINUE RETURN END C SUBROUTINE SOLVN(NN,LHB,AA,S) C* THIS SUBROUTINE SOLVES THE SYSTEM OF EQUATIONS USING C THE DECOMPOSED MATRIX FROM THE PREVIOUS SUBROUTINE IMPLICIT REAL*8(A-H,0-Y) DIMENSION AA(672),S(84) 94 C FORWARD SUBSTITUTION KB = LHB-1 S(1 ) = S ( 1 )/AA( 1 ) IF(NN.EQ.1) GO TO 685 DO 680 I = 2, NN 11 = 1-1 KO = 1 IF(I.GT.LHB) KO=I-KB SUM = S(I) II = LHB*I-KB DO 67 5 K = KO, 11 IK = KB*K+I-KB 675 SUM = SUM-AA(IK)*S(K) S(I) = SUM/AA(II) 680 CONTINUE C BACKWARD SUBSTITUTION 685 Nl = NN-1 LB = LHB*NN-KB S(NN) = S(NN)/AA(LB) IF(NN.EQ.1) RETURN DO 699 I = 1, N1 I 1 = NN-I + 1 NI = NN-I KO = NN IF (I.GT.KB) KO=NI+KB SUM = S(NI) II = LHB*NI-KB DO 690 K = I 1, KO IK = KB*NI+K-KB 690 SUM = SUM-AA(IK)*S(K) S(NI) = SUM/AA(II) 699 CONTINUE RETURN END C SUBROUTINE CONVRG(XO,X,IER,NEQ,EPSLON,ITER) C* THIS SUBROUTINE CHECKS THE CONVERGENCE C OF SOLUTION VECTOR IMPLICIT REAL*8(A-H,0-Y) COMMON/C2/DIFP,NINT DIMENSION XO(84),X(84) IER = 0 PARXO = O.0D0 PARDIF = 0.0D0 PARX = 0.0D0 DO 602 1 = 1 , NEQ PARXO = PARXO + XO(I)**2 PARX = PARX + X(I)**2 602 PARDIF = PARDIF + (X(I)-XO(I))**2 IF (NINT.EQ.1) WRITE(*,1002) PARXO, PARX, PARDIF 1002 FORMATC NORMX0=',E13.6,'NORMX=',E13.6,'NORMDIF=',E13.6/) IF (ITER.EQ.0) GO TO 606 IF (PARDIF.GE.DIFP) GO TO 605 606 DIFP = PARDIF IF (PARXO.EQ.0.0DO) GO TO 603 DIF = DSQRT(PARDIF/PARX0) IF (DIF.LE.EPSLON) GO TO 604 603 IER = 1 95 RETURN 604 RETURN 605 IER = 2 RETURN END SUBROUTINE TIME(TIM) CALL GETTIM(IH,IM,IS,IHS) TIM = IH*3600 + IM*60 + IS + IHS/100.0 RETURN END SAMPLE INPUT/OUTPUT F I L E S . SAMPLE INPUT DATA FILE FOR AXIAL COMPRESSION 10 5 1 10 1 0 0.0 -0.001 2 1 2 1 3 1 1 1 3 32300.0 30350.0 2.0 10.0 5.0 9660000.0 -1.0 3.2 0.038 0.089 SAMPLE OUTPUT FILE FAILURE BETWEEN LOADS 0.l99730e+02 AND 0.201354e+02 AVERAGE = 0.200542e+02 EDGE STRESS (+) = 0.701410e+04 EDGE STRESS (-) = -0.188233e+05 MAX. DEFLECTION = 0.312672e-0l AT LOAD = 0.199730e+02 SLENDERNESS = 3 5.96 98 SAMPLE INPUT DATA FILE FOR PURE BENDING 8 5 10 1 0 1 0.0 5 1 .0 2 1 2 1 3 9 1 3 32300.0 30350.0 2.0 10.0 5.0 9660000.0 -1.0 3.2 0.038 0.089 SAMPLE OUTPUT FILE FAILURE BETWEEN LOAD FACTORS 0.406250e+0l AND 0.409375e+01 AVERAGE = 0.407813e+0l EDGE STRESS (+) = 0.646474e+05 EDGE STRESS (-) = -0.643671e+05 MAX. DEFLECTION = 0.128601e+00 AT LOAD FACTOR = 0.406250e+0l PROGRAM RBETA.FOR (SOURCE CODE) 1 0 0 c c c c c c c c c c c c c c c c c c c c c c c c c c c c ***************************************** * RBETA.FOR V e r s i o n 2.0 * (SHORTENED VERSION WITH SIZE EFFECTS CONSIDERED) * 15 August, 1987 * A PROGRAM FOR THE EVALUATION OF THE REIABILITY INDEX * BETA OF A COLUMN (OR BEAM-COLUMN) * MATERIAL BEHAVIOUR IS ELASTIC IN TENSION WITH BRITTLE FRACTURE, AND ELASTIC IN COMPRESSION UP TO A LIMITING COMPRESSION STRESS, WITH A FALLING LINEAR BRANCH BEYOND THAT LIMIT. END LOAD IS APPLIED CENTRALLY. LATERAL LOADS CAN BE DISTRIBUTED ALONG THE LENGTH OF THE MEMBER THE PROGRAM FINDS THE RELIABILITY INDEX BETA FOR A * BEAM-COLUMN TAKING INTO ACCOUNT 5 RANDOM VARIABLES * WHICH CONSTITUTE THE LOAD AND MATERIAL RESISTANCES * PROBLEM DATA IS READ FROM UNIT #1 * OUTPUT IS STORED IN UNIT #2. * * * ********************************************************* * MAXIMUM OF 10 VARIABLES * * MAXIMUM OF 20 ELEMENTS * IMPLICIT REAL*8(A-H,0-Z) REAL*8 INVNPR,NORMPR DIMENSION X(10),Y(10),U(10),DELTA(10),SIG(10) 1 ,AVER(10),STD(10),F1X(10),F2X(10) 2 ,SCALE(10),SHAPE(10),A(10),B(10),X0(10),XW(10) COMMON/CXI/GAP(5),GAW(5),EN1(8,5),EM1(8,5),EM2(8,5),NGAUSS COMMON/C2/F11,F21 COMMON/C4/W,H,SPAN,PLN,GAMA1,SREF,XKC,XKT COMMON/CX4/NELEM,NBC(21),IX(21,4) REAL*8 LOC(10),MU(10),NN(10), NNN(10) INTEGER*2 ICODEOO) INTEGER*2 MXC(10), MEX(10) OPEN(1,FILE='DET',STATUS='OLD') OPEN(2,FILE='OT',STATUS='NEW') PI 2 = DSQRT(8.0*DATAN(1.0D0)) CONST=1.0D0/PI2 C Q ********************************************************** C * DEFINE VARIABLES * Q ********************************************************** c * FCN COMPRESSIVE STRENGTH, FIFTH PERCENTILE * c * W WIDTH OF SECTION * c * H DEPTH OF SECTION * c * GAMA 1 RATIO OF NOMINAL DL TO LL * c * ALFD DEAD LOAD FACTOR * c * ALFL LIVE LOAD FACTOR * c * EMIN MODULUS OF ELASTICITY, MEAN VALUE * c * NELEM = NO OF ELEMENTS * c * NJBC NO OF JOINTS WITH B.C. * c * NGAUSS = NO OF INTEGRATION POINTS * 101 c * NBC(I) = NO OF B.C. AT NODE I * c * IX = B.C. CODE * c * 1 = U * c * 2 = UX * c * 3 = W * c * 4 = WX * c * N = NO OF RANDOM VARIABLES FOR TOTAL PROB. * c * EN1 ,EM1,EM2 = INTERPOLATION FUNCTIONS * c * GAP = CORDINATE AT GAUSS POINT * c * GAW = CORRESPONDING WEIGHT * c * NELEM = NO OF ELEMENTS * c * NGAUSS = NO OF GAUSS POINTS * c * NITER = MAX. NO OF ITERATIONS * c * TOP = TOLERANCE FOR LOAD • * c * EPSLON = TOLERANCE FOR SOLUTION VECTOR * c * FC = MATERIAL STRENGTH IN COMPRESSION * c * FT = MATERIAL STRENGTH IN TENSION * c * EO = MOE OF THE MATRIAL * c * EN = SLOPE OF THE STRESS-STRAIN CURVE * c * SPAN = MEMBER LENGTH * c * W = WIDTH OF SECTION * c * H = DEPTH OF SECTION * c * E = ECCENTRICITY OF AXIAL LOAD * c * NEQ = NO OF EQUATIONS TO BE SOLVED * c * NJOINT = NO OF NODES * c * NDOF = NO OF VARIABLES PER NODE * c * NODEL = NO OF NODES PER ELEMENT * c * SREF = REFERENCE SPAN * c * XKC = SIZE EFFECT SHAPE PARAMETER (COMP .) c * XKT = SIZE EFFECT SHAPE PARAMETER (TENS .) c * NDIMB = NO OF VARIABLES PER NODE * c * LBW,LHB = HALF BANDWIDTH INCLUD. THE DIAG. * c * NA = NO OF UNKNOWNS FOR TOTAL PROBLEM * Q ***************************** *•* *************************** C READ(1,*) FCN,W,H READ(1,*) GAMA1,ALFD,ALFL READ(1,*) EMIN,SREF,XKC,XKT READ(1,*) NELEM,NJBC,NGAUSS C C* READ BOUNDARY CONDITION CODES FOR THE PROBLEM DO 4455 1 = 1 , (NELEM+1) NBC(I) = 0 DO 4433 J = 1 , 4 IX(I,J) = 0 443 3 CONTINUE 4455 CONTINUE DO 9922 K = 1, NJBC READ(1,*) NJ,NBC(NJ) READ(1,*) (IX(NJ,JV),JV=1,NBC(NJ)) 992 2 CONTINUE C READ(1,*) N, (ICODE(I), I = 1,N) READ(1,*) (MXC(I), I = 1,N) DO 7779 I = 1,N MEX(I) = 0 IF (MXC(I).EQ.0) GO TO 7779 GO TO 7780 7779 CONTINUE GO TO 7782 7780 WRITE(*,7784) 7784 FORMAT(' ENTER EXPONENTS FOR DISTRIBUTION OF EXTREMES'/) READ(*,*) (MEX(I), 1=1,N) 7782 CONTINUE READ(1,*) TOLB READ(1,*) NITER C C ENTER THE CODES FOR EACH VARIABLE AND THEIR PARAMETERS C DO 9 IC = 1,N ICD = ICODE(IC) GO TO(11,12,13,14),ICD C C NORMAL ( CODE=1) C 11 READ(1,*) AVER(IC),STD(IC) GO TO 9 C C LOGNORMAL (CODE=2) C 12 READ(1,*) AVER(IC),STD(IC) GO TO 9 C C WEIBULL ( CODE=3 ) C 13 READ(1,*) LOC(IC),SCALE(IC),SHAPE(IC) GO TO 9 C C GUMBEL EXTREME TYPE I ( CODE=4 ) C 14 READ(1,*) B(IC),A(IC) C 9 CONTINUE C C ENTER INITIAL VECTOR X AND CHECK FOR CONSISTENCY IN THE CASE C OF THE WEIBULL DISTRIBUTION C DO 805 I = 1, N 1 51 READ(1,*) X(I) IF(ICODE(I).NE.3) GO TO 805 IF (X(I ) .GT.LOCd ) ) GO TO 805 WRITE(*,1270) 1270 FORMAT (' CHANGE INITIAL VALUE TO EXCEED THE' 1 ,/,' LOCATION PARAMETER FOR THE WEIBULL'/) GO TO 151 805 CONTINUE C DO 702 I = 1 , N 702 X0(I) = X(I) 155 NCOUNT=0 NBET = 0 I ERR 1 = 0 IERR = 0 READ(1,*) SPAN,R DELT = SPAN/(2.0D0*NELEM) 1 0 3 C C CALC VECTORS EN1, EM1, EM2 CALL SHAP(DELT) C C CORRECT FOR SLENDERNESS EFFECTS PC = W*H*FCN PCR = (3.14159D0**2)*EMIN*(W*H**3)/(12.D0*SPAN**2) CC = SPAN/H CA = DSQRT(0.9D0*0.74D0*EMIN/FCN) CK1 = 1.D0 - (1.D0/3.0D0)*((CC/CA)**4) CK2 = 3.14159D0**2*0.74D0*EMIN/(12.0*FCN*CC**2) IF (CC .GT. 10.0D0) GO TO 2080 CK = 1.0D0 GO TO 4080 2080 IF (CC .GT. CA) GO TO 3080 CK = CK1 GO TO 4080 3080 CK = CK2 4080 CONTINUE C C OBTAIN NOMINAL DESIGN LOAD PLN = R*W*H*FCN*CK/(ALFD*GAMA1+ALFL) C WRITE(2,1080)(ICODE(I), 1=1,N) 1080 FORMAT (' CODES : ',1015) C C START ITERATIONS: GIVEN THE VECTOR X ( l ) , COMPUTE C THE FAILURE FUNCTION GXP AND THE GRADIENT DELTA C USING THE SUBROUTINE GXPR, WHICH MUST BE PROVIDED C EXTERNALLY BY THE USER FOR EACH PARTICULAR CASE. C 2 CONTINUE DO 7722 J = 1, N 7722 XW(J) = X(J) CALL GXPR(XW,N,DELTA,GXP) C C CALC F1X(X), AND F2X(X) C CALL FFX(N,X,AVER,STD,F1X,F2X,ICODE,LOC,SCALE,SHAPE,A,B, 1 IERR, MXC, MEX) IF(IERR.EQ.1) GO TO 65 C C CALC Y-VALUES C DO 8 I = 1 , N Y(I) = INVNPR(F1X(I ) ) 8 CONTINUE C C CALC SIGMA AND MU VECTORS C DO 10 I = 1 , N IF (F2X(I).LE.0.0D0) GO TO 68 DSIG = DLOG(CONST) - Y(I)*Y(I)*0.5D0 - DLOG(F2X(I)) IF (DSIG.LT.-709.0D0) GO TO 865 SIG(I) = DEXP(DSIG) GO TO 87 865 SIG(I) = 0.0D0 87 MU(I)=-SIG(l)*Y(I)+X(l) 10 CONTINUE C C CALC NN SUM=0.0D0 DO 55 1=1,N 55 SUM = SUM + SIG(I)*SIG(I)*DELTA(I)*DELTA(I) SUM = DSQRT(SUM) DO 20 I = 1 ,N NN(I) = -SIG(I)*DELTA(I)/SUM 20 NNN(I) = DABS(NN(I)) C C CALC BETA SDMU=0.0D0 SDX = 0.0D0 DO 25 1=1,N SDMU = SDMU + DELTA(I)*MU(I) 25 SDX = SDX + DELTA(I)*X(I) BETA = (GXP + SDMU - SDX)/SUM DO 30 I = 1 ,N 30 U(I) = BETA * NN(I) NCOUNT = NCOUNT+1 IF (NCOUNT.GT.NITER) GO TO 66 IF(NCOUNT.EQ.1) GO TO 32 DIFFB = DABS(BETA - BETAP) BETAP = BETA NBET = 1 DO 80 I = 1, N TX = SIG(I)*U(I) + MU(I) 80 X(I) = TX CONFAC = (TOLB-DIFFB) IF (CONFAC.GT.0.0) GO TO 50 GO TO 2 32 BETAP = BETA DO 35 I = 1 ,N 35 X(I) = SIG(I)*U(I) + MU(I) GO TO 2 50 WRITE(2,51) BETA WRITE(2,710) NCOUNT 710 FORMAT(5X,'ITERATIONS =*,I5) WRITE(2,703) TOLB 703 FORMAT(5X,*TOLB =',F8.4) 705 WRITE(2,1280)(X0(I) ,1=1,N) 1280 FORMATC VECTOR XO ' , 1 0E1 3 . 5) WRITE(2,1300)(X(I) ,1=1,N) 1300 FORMAT(' VECTOR X *,10E13.5) WRITE(2,1320)(NNN(I),1=1,N) 1320 FORMAT(' SENSITIVITY COEFFS. ',10F8.4) WRITE(2,2088) SPAN,R 2088 FORMAT(' L=',E13.6,' fp = ' , E l 3 . 6 / ) 51 FORMAT(5X,'BETA = ',F10.3) GO TO 900 65 IF (NBET.EQ.1) GO TO 880 WRITE(2,1340) I ERR GO TO 900 880 WRITE( 2,1340 ) I ERR GO TO 900 68 IERR1 = 1 IF (NBET.EQ.1) GO TO 882 WRITE(2,1341) IERR1 GO TO 900 882 WRITE(2,1341 ) IERR 1 WRITE(2,1342) BETA GO TO 900 66 WRITE(2,1350)NITER 1350 FORMAT (' NO CONVERGENCE IN *,15,' ITERATIONS') 1340 FORMAT(' IERR =',12,' ERROR: NEGATIVE LOGNORMAL OR',/, 1 ' WEIBULL VARIABLE LESS THAN ITS',/, 2 ' LOCATION PARAMETER.',/, 3 ' TRY NEW INITIAL POINT') 1341 FORMAT(' I ERR 1 =',12,' NEGATIVE OR ZERO DENSITY F 2 X ( l ) ' , / , 1 ' TRY NEW INITIAL POINT'/) 1342 FORMAT(' LAST BETA WAS =', F10.3) 900 CONTINUE CLOSE (UNIT=1,STATUS='KEEP') CLOSE (UNIT=2,STATUS=*KEEP') STOP END C SUBROUTINE FFX(N,X,AVER,STD,F1X,F2X,ICODE,LOC,SC,SK,A,B, 1 IERR, MXC, MEX) IMPLICIT REAL*8(A-H,0-Z) REAL*8 NORMPR DIMENSION SC(N),SK(N),A(N),B(N),X(N),AVER(N) 1, STD(N),F1(N),F2X(N) COMMON/C2/F11,F21 INTEGER*2 ICODEOO) INTEGER*2 MXC(10), MEX(lO) REAL*8 LOC(N),MU PI2=(8.0*DATAN(1.0D0)) DO 20 I = 1 ,N IC = I CODEC I) GO TO(11,12,13,14) , IC C C NORMAL C 11 RATIO = (X(I) - AVER (I) ) /STD(I) F1X(I) = NORMPR(RATIO) F2X(I)=DEXP(-0.5D0*RATIO*RATIO)/(STD(l)*DSQRT(PI 2)) IF (MXC(I).EQ.0) GO TO 20 CALL EXTR (F 1X (I ) , F2X (I ) , MXC (I ) , MEX(D) GO TO 20 C C LOGNORMAL C 12 DLN = DLOGd.O + ( STD (I )/AVER (I ) ) ** 2 ) MU = DLOG(AVER(I)) - 0.5*DLN SDP = DSQRT(DLN) IF (X(I).LE.0.0) GO TO 99 PARAM = (DLOG(X(I)) - MU)/SDP F1X(I) = NORMPR(PARAM) POW = DEXP(-0.5D0*PARAM*PARAM) F2X(I) = POW/(SDP*X(I)*DSQRT(PI2)) IF (MXC(I).EQ.O) GO TO 2 0 CALL EXTR(FIXd)., F 2 X ( l ) , MXC (I ) , MEX(I)) GO TO 20 106 C WEIBULL C 13 IF (X(I).LE•LOC(I)) GO TO 99 POW = - ( ( X ( I ) - L O C ( l ) ) / S C ( I ) ) * * S K ( I ) POW = DEXP(POW) F1X(I) = 1.ODO - POW F2X(I) = ( S K ( I ) / S C ( I ) ) * ( ( X ( I ) - L O C ( l ) ) / S C ( l ) ) * * ( S K ( I ) 1- 1.0)*POW GO TO 20 . F2X(I) IF (MXC(I).EQ.0) CALL EXTR(F1X(I) GO TO 20 C C GUMBEL EXTREME TYPE I C 14 POW = - A ( I ) * ( X ( I ) -POW = DEXP(POW) F1X(I) = DEXP(-POW) F2X(I) = A(I)*POW * IF (MXC(I).EQ.0) CALL EXTR(F1X(I) 20 CONTINUE RETURN 99 IERR=1 RETURN END MXC(I), MEX(I)) B(I ) ) F1X(I) GO TO 20 F2X(I), MXC(I), MEX(I)) 10 SUBROUTINE EXTR(F1,F2,NC,M) IMPLICIT REAL*8(A-H,0-Z) INTEGER*2 NC, M IF (NC.EQ.2) GO TO 10 F2 = M*F2*F1**(M-1) F1 = F1**M RETURN F2 = M*F2*(1.0D0 - F1)**(M-1) F1 = 1.0D0 - (1.0D0 - F1)**M RETURN END E( 1 E(2 E(3 E(4 E(5 E(6 E(7 E(8 H(1 H(2 H(3 H(4 H(5 FUNCTION NORMPR(X) * NORMAL PROBABILITY INTEGRAL IMPLICIT REAL*8(A-H,0-Z) REAL*8 NORMPR DIMENSION E(16),H(16) PI = 2.ODO * DSQRT(DATAN(1.0D0)) IF (DABS(X).GT.5.0D0) GO TO 20 0.989400934991650EO 0.944575023073233EO 0.865631202387832E0 0.755404408355003E0 0.617876244402644E0 0.458016777657227E0 0.281603550779259E0 0.095012509837637E0 0.027152459411754E0 0.062253523938648E0 0.095158511682493E0 0.124628971255534E0 0. 1 49595988816577E0 107 H(6) = 0.169156519395003E0 H(7) = 0.182603415044924E0 H(8) = 0.189450610455068E0 DO I I = 1,8 E{ 1 7-1 ) = -E(I) 1 H(17-I) = H(I) Y = X/DSQRT(2.0D0) S = 0.0 DO 10 I = 1 , 16 Z = Y * E(I) Z = DEXP(-Z*Z) S = S + Z*H(I) 10 CONTINUE ERF = Y * S/PI NORMPR = (1.0D0 + ERF)/2.0D0 RETURN 20 IF (DABS(X).GT.37.5D0) GO TO 25 S = 1.0D0 - 1.0D0/(X**2) + 3.0D0/(X**4) - 15.0D0/(X**6) 1 + 105.0D0/(X**8) - 945.0D0/(X**10) + 10395.0D0/(X**12) S = S*DEXP(-X*X/2.0D0)/DABS(X) S = S*DSQRT(2.0D0)/PI IF (X.GT.0.0D0) NORMPR = 1.ODO - S/2.ODO IF (X.LT.0.0D0) NORMPR = S/2.ODO RETURN 25 IF (X.GT.0.0D0) NORMPR = 1.ODO IF (X.LT.0.0D0) NORMPR = 0.ODO RETURN END C C FUNCTION INVNPR(Y) C * INVERSE NORMAL PROBABILITY * C IMPLICIT REAL*8(A-H,0-Z) REAL*8 INVNPR REAL*8 NORMPR PI = DSQRT(8.0D0*DATAN(1.ODO)) TOL = 1.OE-8 IF (Y.EQ.0.50) GO TO 80 XO = -PI*(0.50D0 - Y) X1 = XO 5 S = NORMPR(X1) - Y S = S * DEXP(X1*X1/2.0D0) * PI X2 = X1 - S DIF = DABS(X2-X1) IF (DABS(DIF).LE.TOL) GO TO 20 XI = X2 GO TO 5 80 INVNPR = 0.0 RETURN 20 INVNPR=X2 RETURN END C SUBROUTINE COLUMN(XW,N,PAV) IMPLICIT REAL*8(A-H,0-Z) DIMENSION F(8),TK0(672),XE(8),XW(N) 1,R(84),X0(84),X(84),B(84),B1(8,8),B2(8,8),B3(8,8),B4(8,8) 108 2, B5(8,8),B6(8,8),B7(8,8),B8(8,8),B9(8,8),Y(5),RE(8),XP(84) 3, Q(20),IQ(20), ESTR(7), FI(7) C0MM0N/CX1/GAP(5),GAW(5),EN1(8,5),EM1(8,5),EM2(8,5),NGAUSS COMMON/CX2/DIFP,NINT COMMON/C3/DEFL,PDEFL COMMON/C 4/W,H,SPAN,PLN,GAMA1,SREF,XKC,XKT COMMON/CX4/NELEM,NBC(21),IX(21,4) c ********************************************** c * EN I,EM1 ,EM2 = INTERPOLATION FUNCTIONS * c * GAP = CORDINATE AT GAUSS POINT * c * GAW = CORRESPONDING WEIGHT * c * NELEM = NO OF ELEMENTS * c * NGAUSS = NO OF GAUSS POINTS * c * NITER = MAX. NO OF ITERATIONS * c * TOP = TOLERANCE FOR LOAD * c * EPSLON = TOLERANCE FOR SOLUTION VECTOR * c * FC = MATERIAL STRENGTH IN COMPRESSION * c * FT = MATERIAL STRENGTH IN TENSION * c * EO = MOE OF THE MATRIAL * c * EN = SLOPE OF THE STRESS-STRAIN CURVE * c * SPAN = MEMBER LENGTH * c * W = WIDTH OF SECTION * c * H = DEPTH OF SECTION * c * E = ECCENTRICITY OF AXIAL LOAD c * NEQ = NO OF EQUATIONS TO BE SOLVED * c * NJOINT = NO OF NODES * c * NDOF = NO OF VARIABLES PER NODE * c * NODEL = NO OF NODES PER ELEMENT * c * SREF = REFERENCE SPAN * c * XKC = SIZE EFFECT SHAPE PARAMETER (COMP. ) * c * XKT = SIZE EFFECT SHAPE PARAMETER (TENS. ) * c * NDIMB = NO OF VARIABLES PER NODE * c * LBW,LHB = HALF BANDWIDTH INCLUD. THE DIAG. * c * NA = NO OF UNKNOWNS FOR TOTAL PROBLEM * c ******************************************************** CONST = GAMA1 EO = XW(1)*1000.DO FC = XW(2)*1000.DO FT = XW(3)*1000.DO TOP = 0.01 DO EPSLON = 0.001 DO NITER = 10 EN = 0.02D0 NDOF = 4 NJOINT=NELEM+1 NG1 = NGAUSS + 1 NG2 = NGAUSS + 2 NP = 1 NQ = 0 Q0 = 0.0D0 IF (NQ.EQ.0) GO TO 44 DO 43 I = 1, NQ 43 READ(1,*) I Q ( I ) , Q (l) 44 ECEN = -0.002D0 IF (NP.NE.0) E=ECEN NEQ = NDOF*NJOINT NODEL = 2 NDIMB = NODEL*NDOF 109 LBW = NDIMB LHB = LBW NA = LBW*NEQ AR = W*H XI = W*H**3/12.D0 DEL == SPAN/(2.D0*NELEM) C C * ADJUST STRENGTHS TO THE ACTUAL VOLUME FC = FC *(SREF/SPAN)**(1.0/XKC) FT = FT *(SREF/SPAN)**(1.0/XKT) NINT = 0 IF (NP.BQ.O) GO TO 761 PC = AR*FC PCR = 3.14159D0**2*E0*XI/(SPAN**2) PI = PC IF(PCR .LE. PC) PI=PCR P2 = PI PI = 0.0D0 P3 = (PI + P2)/2.0D0 NFAIL = 0 SMAX1 = 0.0 GO TO 760 761 FQ1 = 0.0D0 FQ2 = 1.0D0 FQ3 = FQ2 NFLAG = 0 760 DO 792 J = 1, NEQ 792 XP(J) = O.ODO C C START CALCULATIONS FOR TRIAL LOAD LEVELS 3773 CONTINUE P = O.ODO FQ = 1.0D0 IF (NP.NE.O) P = P3 IF (NP.EQ.O) FQ = FQ3 IF (NINT.EQ.1.AND.NP.NE.O) WRITE(*,4000) P IF (NINT.EQ.1.AND.NP.EQ.O) WRITE(*,4001) FQ 4000 FORMAT(//' SOLUTION FOR P =',E15.6,' :'/) 4001 FORMAT(//* SOLUTION FOR LATERAL LOAD FACTOR= 1, E1 5. 6 , ':'/) C INITIALISE ARRAYS DO 80 J = 1, NEQ XO(J) = XP(J) 80 R(J) = 0.D0 C C EXTERNAL LOAD VECTOR R C RE = ELEMENT LOAD VECTOR IF (QO.EQ.O.ODO) GO TO 87 DO 81 J = 1 , 8 RE(J)=0.DO 81 CONTINUE RE(3) = FQ*Q0*DEL RE(4) = FQ*Q0*DEL**2/3.D0 RE(7) = RE(3) RE(8) = -RE(4) DO 83 NE=1, NELEM DO 82 J J = 1, 8 K = (NE-1)*NDOF + JJ R(K) = R(K) + RE(JJ) 110 82 CONTINUE 83 CONTINUE 87 IF (NQ.EQ.O) GO TO 185 DO 180 J = 1 , NQ JS = (IQ(J)-1)*NDOF + 3 180 R(JS) = R(JS) + Q(J)*FQ 185 EM = P*E ' J J = (NJOINT-1)*NDOF + 1 R(JJ) = R(JJ)-P R(1) = R(1) + P R(4) = R(4)-EM R(NEQ) = R(NEQ)+EM ITER = 0 C C BEGIN ITERATIONS AT THE TRIAL LOAD LEVEL 777 CONTINUE DO 84 I = 1 , NA 84 TKO(I) = O.ODO DO 85 K = 1, NEQ 85 B(K) = -R(K) DO 645 IE = 1, NELEM C INITIALIZE ARRAYS DO 88 I = 1 , 8 F(I) = O.ODO 86 88 DO 86 J = 1 , I B l ( I , J ) = = 0 .0D0 B2(I,J) = = 0 .0D0 B3(I,J) = = 0 .0D0 B4(I,J) = = 0 .0D0 B5(I,J) = = 0 .0D0 B6(I,J) = = 0 .0D0 B7(I,J) = = 0 .0D0 B8(I,J) = = 0 .0D0 B9(I,J) = = 0 .0D0 CONTINUE CONTINUE C PICK ELEMENT SOLUTION FROM GLOBAL VECTOR DO 90 J J = 1 , 8 K = (IE - 1)*NDOF + J J XE(JJ) = XO(K) 90 CONTINUE DO 101 K = 1, NGAUSS Y(K) = 0.D0 DO 91 1=1 , 8 Y(K) = Y(K) + XE(I)*EM1(I,K) 91 CONTINUE C OBTAINING COMPONENTS OF EKT DO 93 I = 1 , 8 DO 93 J = 1 , I B1(I,J) = B1(I,J)+E0*DEL*EN1(I,K)*Y(K)*AR* 1 EM1(J,K)*GAW(K) B2(I,J) = B2(I,J)+E0*DEL*EM1(I,K)*Y(K)*AR* 1 EN 1(J,K)*GAW(K) B3(I,J) = B3(I,J)+E0*DEL*EM1(I,K)*Y(K)*AR* 1 Y(K)*EM1(J,K)*GAW(K) B4(I,J) = B4(I,J)+(E0*AR*DEL*EN1(I,K)*EN1(J , K ) + 1 E0*XI*DEL*EM2(I,K)*EM2(J,K))*GAW(K) 93 CONTINUE 111 DO 100 L = 1 , NGAUSS C STRESSES AND STRAINS AT GAUSS POINT STR = 0.5D0*Y(K)**2 DO 96 MO = 1, 8 STR = STR+(EN1(MO,K)-GAP(L)*H*0.5D0*EM2(MO,K))*XE(MO) 96 CONTINUE 1 STRE = STR+FC/EO FAC = 1.ODO IF(STRE.GE.O.DO) FAC=0.0D0 STRESS = E0*STR-((E0+EN*E0)*STR+FC*(1.DO+EN))*FAC DO 99 I = 1, 8 DO 98 J = 1 , I B5(I,J) = B5(I,J)+DEL*0.5D0*AR*(EN 1(I,K)-GAP(L)* 1 H*0.5D0*EM2(I,K))*(E0+E0*EN)*FAC*(EN 1(J,K)-H*0.5D0* 2 GAP(L)*EM2(J,K))*GAW(K)*GAW(L) B6(I,J) = B6(I,J)+DEL*0.5D0*AR*(EN 1(I,K)-GAP(L)* 1 H*0.5D0*EM2(I,K))*(E0+EN*E0)*FAC*Y(K)*EM1(J,K)* 2 GAW(K)*GAW(L) B7(I,J) = B7(I,J)+DEL*0.5DO*EM1(I,K)*Y(K)*AR* 1 (E0+E0*EN)*FAC*(EN1(J,K)-H*0.5D0*GAP(L)*EM2(J,K))* 2 GAW(K)*GAW(L) B8(I,J) = B8(I,J)+DEL*0.5D0*EM1(I,K)*Y(K)*AR* 1 (E0+E0*EN)*FAC*Y(K)*EM1(J,K)*GAW(K)*GAW(L) B9(I,J) = B9(I,J)+AR*STRESS*EM1(I,K)*EM1(J,K)* 1 GAW(K)*GAW(L)*DEL*0.5D0 98 CONTINUE F( I ) = F(I)+AR*DEL*0.5D0*STRESS*((EN1(I,K)-H*0.5D0* 1 GAP(L)*EM2(I,K))+Y(K)*EM1(I,K))*GAW(K)*GAW(L) 99 CONTINUE 100 CONTINUE 101 CONTINUE C* OBTAIN ELEMENT TANGENT MATRIX C EKT IS THE (I,J) COMPONENT OF THE ELEMENT TANGENT MATRIX DO 105 I = 1, 8 II = (IE-1)*NDOF + I B(II) = B(II) + F ( I ) DO 102 J = 1, I J J = (IE-1)*NDOF + J EKT = B1(I,J)+B2(I,J)+B3(I,J)+B4(I,J)-1 B5(I,J)-B6(I,J)-B7(I,J)-B8(I,J)+B9(I,J) I J = (JJ-1)*(LBW-1) + II TKO(IJ) = TKO(IJ)+EKT 102 CONTINUE 105 CONTINUE 645 CONTINUE C INTRODUCE BOUNDARY CONDITIONS DO 111 IJO = 1, NJOINT IF (NBC(IJO).EQ.0) GO TO 111 DO 110 J = 1, NBC(IJO) II = ( U O -l)*NDOF + IX(UO,J) LBW1 = LBW - 1 DO 108 K = 1, LBW1 J J = II - LBW + K IF (JJ.LE.0) GO TO 1080 I J = (JJ-1)*(LBW-1) + II TKO(IJ) = 0.0D0 1080 J J = 11 + K IF (JJ.GT.NEQ) GO TO 108 I J = {11- 1)*(LBW-1) + J J TKO(IJ) = O.ODO 108 CONTINUE I J = (II - 1)*(LBW-1) + II TKO(IJ) = 1.0DO B(II) = O.ODO i 110 CONTINUE 111 CONTINUE C C SOLUTION OF THE SYSTEM CALL DECOMP(NEQ,LBW,TKO,IERROR) IF(IERROR .EQ. 1) GO TO 3774 CALL SOLVN(NEQ,LBW,TKO,B) DO 112 I = 1, NEQ X(I) = XOU)-B(I) 112 CONTINUE CALL CONVRG(XO,X,IER,NEQ,EPSLON,ITER) ITER = ITER + 1 IF (ITER.EQ.NITER) GO TO 431 IF (IER.EQ.2) GO TO 430 IF(IER.EQ.O) GO TO 118 DO 115 I = 1, NEQ 115 XO(I) = X(I) GO TO 777 430 IERROR = 1 GO TO 3774 431 WRITE(2,900) NITER, P 900 FORMAT(' NO CONVERGENCE IN',13,' ITERATIONS AT P=',E13.6/) GO TO 901 C* AFTER CONVERGENCE, OBTAIN STRESSES AND STRAINS C AT THE CURRENT LOAD LEVEL C 118 CONTINUE EMAXP = O.ODO EMAXN = O.ODO SUME = O.ODO DO 550 IE = 1, NELEM DO 500 J = 1, 8 K = (IE-1)*NDOF +J XE(J) = X(K) 500 CONTINUE DO 540 K = 1, NGAUSS FACTOR = 0.0 DO 501 I = 1, 8 501 FACTOR = FACTOR + XE(I)*EM1(I,K) EPLUS = 0.5D0 * FACTOR**2 EMINUS = EPLUS DO 505 I = 1, 8 EPLUS = EPLUS + (EN 1 (I,K)-H*0.5D0*EM2(I,K))*XE(I) EMINUS = EMINUS + (EN 1(I,K)+H*0.5D0*EM2(I,K))*XE(I) 505 CONTINUE IF(EPLUS.GT.O.ODO .AND. EMINUS.GT.0.0) GO TO 506 IF(EPLUS.GT.O.ODO .AND. EMINUS.LE.0.0) GO TO 507 IF(EPLUS.LE.0.0D0 .AND. EMINUS.LE.0.0) GO TO 508 IF(EPLUS.LE.O.ODO .AND. EMINUS.GT.0.0) GO TO 509 506 EPOS = EPLUS IF(EMINUS.GT.EPOS) EPOS=EMINUS ENEG = O.ODO 1 1 3 GO TO 530 507 EPOS = EPLUS ENEG = EMINUS GO TO 510 508 EPOS = O.ODO ENEG = EPLUS IF (DABS (EMINUS ) .GT.'DABS (ENEG) ) ENEG = EMINUS GO TO 530 509 EPOS = EMINUS ENEG = EPLUS C C * FINDS THE POSITION OF THE NEUTRAL AXIS 510 ESTR(1) = EMINUS F l ( 1 ) = -1.0D0 ESTR(NG2) = EPLUS Fl(NG2) = 1.0D0 DO 512 L = 1, NGAUSS SUM = 0.5*FACTOR**2 DO 511 I = 1,8 511 SUM = SUM + (EN1(I,K) - GAP(L)*H/2.0*EM2(I,K))*XE(I) ESTR(L+1) = SUM Fl(L+1) = GAP(L) 512 CONTINUE DO 515 I = 1, NG1 PROD = ESTR(I)*ESTR(I+1) IF (PROD.LE.O.ODO) GO TO 516 515 CONTINUE 516 XN = F I ( I ) - ESTR(I)*(FI(1+1)-FI(I))/(ESTR(I+1)-ESTR(l)) IF (ESTR(I).EQ.O.ODO) GO TO 518 IF (ESTR(I).LT.O.ODO) HN = (1.0D0 - XN)*H/2.0D0 IF (ESTR(I).GT.O.ODO) HN = (1.0D0 + XN)*H/2.0D0 GO TO 520 518 IF (ESTR(1+1).LT.O.ODO) HN = (1.0D0 + XN)*H/2.0D0 IF (ESTR(I+1).GT.O.ODO) HN = (1.0D0 - XN)*H/2.0D0 520 SUME = SUME + (HN/H)*(E0*EPOS)**XKT*GAW(K) 530 IF(EPOS.LT.EMAXP) GO TO 538 EMAXP = EPOS 538 IF(DABS(ENEG).LT.DABS(EMAXN)) GO TO 540 EMAXN = ENEG 540 CONTINUE 550 CONTINUE SMAXP = E0*EMAXP SMAXN = E0*EMAXN IF (DABS(SMAXN).LE.FC) GO TO 560 SMAXN = SMAXN -((E0 + EN*E0)*EMAXN + FC*(1.0 + EN)) 560 IF (SUME.EQ.O.ODO.OR.SMAXP.EQ.O.ODO) GO TO 563 SUME = SUME/(2.0*NELEM*(XKT+1.0)*SMAXP**XKT) FTT = FT * SUME**(-1.0D0/XKT) GO TO 564 563 FTT = FT 564 IF (SMAXP.GE.FTT) GO TO 3774 DEFL = O.ODO DO 565 IE = 1, NELEM J = (IE-1)*NDOF + 3 IF (DABS(X(J)).GT.DABS(DEFL)) DEFL = X(J) 565 CONTINUE J = NEQ - 1 IF (DABS(X(J)).GT.DABS(DEFL)) DEFL = X(J) 1 1 4 IF (NP.EQ.O) PDEFL = FQ3 IF (NP.NE.O) PDEFL = P3 3774 CONTINUE IF (NINT.EQ.0) GO TO 8810 IF (IERROR.EQ.1) WRITE(*,8888) IF (I ERROR. EQ.0.A.ND. SMAXP. LT. FTT) WRITE (*, 8889 ) SMAXP IF (IERROR.EQ.0.AND.SMAXP.GE.FTT) WRITE(*,8890) SMAXP 8888 FORMAT(' IERROR = 1,FAILS (DIVERGENCE OR SINGULAR MATRIX)'/] 8889 FORMAT(' IERROR = 0 SMAXP = *,E15.6,' SURVIVES'/) 8890 FORMAT(' IERROR = 0 SMAXP = ',E15.6,' FAILS'/) 8810 CONTINUE IF (NP.EQ.O) GO TO 4500 IF (IERROR.EQ.1) GO TO 7330 IF (SMAXP.GT.FTT) GO TO 7331 . IF (SMAXP.EQ.FTT) GO TO 7337 PI = P3 IF (SUME.EQ.0.ODO.OR.SMAXP.EQ.0.ODO) GO TO 5650 SMAX1 = SMAXP*SUME**(1.0D0/XKT) GO TO 5655 5650 SMAX1 = SMAXP 5655 DO 833 J = 1, NEQ 833 XP(J) = X(J) GO TO 8334 7330 P2 = P3 GO TO 8334 7331 P2 = P3 NFAIL = 1 SMAX2 = SMAXP*SUME**(1.0D0/XKT) 8334 IF (PI.EQ.0.ODO) GO TO 8338 TOLP = (P2-P1)/P1 IF (TOLP.LE.TOP) GO TO 7338 GO TO 8336 8338 IF (P2.LE.0.1D0) GO TO 7338 8336 IF (NFAIL.EQ.1) GO TO 8340 P3 = (PI + P2)/2.0 GO TO 3773 8340 P3 = P1 + (P2-P1)*(FT-SMAX1)/(SMAX2-SMAX1) GO TO 3773 7337 P = P3 PP = P3 PAV = P3 GO TO 7339 7338 IF (PI.EQ.0.ODO) P2 = 0.ODO P = P2 PP = PI PAV = (P+PP)/2.0 7339 CONTINUE GO TO 901 4500 IF (IERROR.EQ.1) GO TO 4330 IF (SMAXP.GT.FTT) GO TO 4330 IF (SMAXP.EQ.FTT) GO TO 4337 IF (NFLAG.EQ.1) GO TO 4331 FQ1 = FQ2 FQ2 = 2.0D0*FQ2 GO TO 4580 4331 FQ1 = FQ3 4580 DO 4833 J = 1,NEQ 4833 XP(J) = X(J) GO TO 4334 4330 NFLAG = 1 FQ2 = FQ3 4334 IF (FQ1.EQ.O.ODO) GO TO 5338 TOLP = (FQ2-FQ1)/FQ1 IF (TOLP.LE.TOP) GO TO 4338 5338 IF (NFLAG.EQ.O) FQ3 = FQ2 IF (NFLAG.EQ.1) FQ3 = (FQ1+FQ2)/2.0D0 GO TO 3773 , 4337 P = FQ3 PP = FQ3 PAV = FQ3 GO TO 4339 4338 P = FQ2 PP = FQ1 PAV = (P+PP)/2.0 4339 CONTINUE 901 RETURN END C SUBROUTINE SHAP(DELT) C* THIS SUBROUTINE CALCULATES DERIVATIVES OF SHAPE FUNCTIONS IMPLICIT REAL*8(A-H,0 -Z) C0MM0N/CX1/GAP(5),GAW(5),EN1(8,5),EM1(8,5),EM2(8,5),NGAUSS IF (NGAUSS.EQ.5) GO TO 5 IF (NGAUSS.EQ.4) GO TO 4 C *** 3 POINT GAUSSIAN INTEGRATION GAP(1) = -0.774596669241483D0 GAP(2) = O.ODO GAP(3) = -GAP(1) GAW(1) = 0.555555555555556D0 GAW(2) = 0.888888888888889D0 GAW(3) = GAW(1) GO TO 10 C *** 4 POINT GAUSSIAN INTEGRATION 4 GAP(1) = -0.861136311594053D0 GAP(2) = -0.339981043584856D0 GAP(3) = -GAP(2) GAP(4) = -GAP(1) GAW(1) = 0.347854845137454D0 GAW(2) = 0.652145154862546D0 GAW(3) = GAW(2) GAW(4) = GAW(1) GO TO 10 C *** 5 POINT GAUSSIAN INTEGRATION 5 GAP(1) = -0.906179845938664DO GAP(2) = -0.538469310105683DO GAP(3) = O.ODO GAP(4) = -GAP(2) GAP(5) = -GAP(1) GAW(1) = 0.236926885056189D0 GAW(2) = 0.478628670499366D0 GAW(3) = 0.568888888888889D0* GAW(4) = GAW(2) GAW(5) = GAW(1) C INITIALISES EN 1 ,EM1 ,EM2 10 DO 150 IL = 1, 8 DO 350 IK = 1, NGAUSS 116 EN 1(IL,IK) = O.ODO EM1(IL,IK) = O.ODO EM2(IL,IK) = O.ODO 350 CONTINUE 150 CONTINUE DO 250 I = 1, NGAUSS ENI(1,I) = (-0.75DO+0.75DO*GAP(I)**2)/DELT EN1(2,I) = (-1.D0-2.D0*GAP(I)+3.D0*GAP(I)**2)*0.25D0 EN1(5,I) = (0.75D0-0.75D0*GAP(I)**2)/DELT EN1(6,I) = (-1.D0+2.D0*GAP(I)+3.D0*GAP(I)**2)*0.25D0 EMI(3,1) = (-0.75D0+0.75D0*GAP(I)**2)/DELT EM1(4,I) = (-1.D0-2.D0*GAP(I)+3.D0*GAP(l)**2)*0.25D0 EM1(7,I) = (0.75D0-0.75D0*GAP(I)**2)/DELT EM1(8,I) = (-1.D0+2.D0*GAP(I)+3.D0*GAP(l)**2)*0.25D0 EM2(3,I) = 1.5D0*GAP(I)/(DELT**2) EM2(4,I) = (~2.D0+6.D0*GAP(I))/(4.D0*DELT) EM2(7,I) = -1.5D0*GAP(I)/(DELT**2) EM2(8,I) = (2.D0+6.D0*GAP(I))/(4.D0*DELT) 250 CONTINUE RETURN END C SUBROUTINE DECOMP(NN,LHB,AA,IERROR) C* THIS SUBROUTINE DECOMPOSES A MATRIX USING CHOLESKY C METHOD FOR BANDED,SYMMETRIC,POS. DEFN. MATRIX IMPLICIT REAL*8(A-H,0-Z) DIMENSION AA(672) C TKO IS STORED COLUMN - WISE. IERROR = 0 KB = LHB-1 C DECOMPOSITION IF(AA(1).LE.0.D0) IERROR=1 IF(IERROR.EQ.1) RETURN AA(1) = DSQRT(AA(1)) IF(NN.EQ.I) RETURN DO 551 I = 2, LHB 551 AA(I) = AA(I)/AA(1) DO 590 J = 2, NN J1 = J-1 IJD = LHB*J-KB SUM = AA(IJD) KO = 1 IF(J.GT.LHB) KO=J-KB DO 555 K = KO, J1 JK = KB*K+J-KB 555 SUM = SUM-AA(JK)*AA(JK) IF(SUM.LE.O.DO) IERROR=1 IF(I ERROR.EQ.1) RETURN AA(IJD) = DSQRT(SUM) DO 568 I = 1, KB II = J + I KO = 1 IF (II.GT.LHB) KO=II-KB SUM = AA(IJD+I) IF(I.EQ.KB) GO TO 565 DO 540 K = KO, J1 JK = KB*K+J-KB IK = KB*K+II-KB 540 SUM = SUM-AA(IK)*AA(JK) 565 AA(IJD+I) = SUM/AA(IJD) 568 CONTINUE 590 CONTINUE RETURN END C SUBROUTINE SOLVN(NN,LHB,AA,S) C* THIS SUROUTINE SOLVES CALLS A MATRIX SOLVER TO THE SYSTEM IMPLICIT REAL*8(A-H,0-Z) DIMENSION AA(672),S(84) C FORWARD SUBSTITUTION KB = LHB-1 S(1) = S(1)/AA(1) IF(NN.EQ.1) GO TO 685 DO 680 I = 2, NN 11 = 1-1 KO = 1 IF(I.GT.LHB) KO=I-KB SUM = S(I) II = LHB*I-KB DO 675 K = KO, I 1 IK = KB*K+I-KB 675 SUM = SUM-AA(IK)*S(K) S(I) = SUM/AA(II) 680 CONTINUE C BACKWARD SUBSTITUTION 685 N1 = NN-1 LB = LHB*NN-KB S(NN) = S(NN)/AA(LB) IF(NN.EQ.I) RETURN DO 699 I = 1, N1 I 1 = NN-I + 1 NI = NN-I KO = NN IF (I.GT.KB) KO=NI+KB SUM = S(NI) II = LHB*NI-KB DO 690 K = I 1, KO IK = KB*NI+K-KB 690 SUM = SUM-AA(IK)*S(K) S(NI) = SUM/AA(II) 699 CONTINUE RETURN END C SUBROUTINE CONVRG(XO,X,IER,NEQ,EPSLON,ITER) C* THIS SUBROUTINE CHECKS THE CONVERGENCE OF SOLUTION VECTOR IMPLICIT REAL*8(A-H,0-Z) COMMON/CX2/DIFP,NINT DIMENSION XO(84),X(84) IER = 0 PARXO = O.ODO PARDIF = O.ODO PARX = O.ODO DO 602 I = 1, NEQ PARXO = PARXO + XO(I)**2 PARX = PARX + X(I)**2 1 18 602 PARDIF = PARDIF + (X(I)-XO(I))**2 IF (NINT.EQ.1) WRITE(*,1002) PARXO, PARX, PARDIF 1002 FORMAT(' NORMX0=',E13.6,'NORMX=',E13.6,'NORMDIF=',E13.6/) IF (ITER.EQ.0) GO TO 606 IF (PARDIF.GE.DIFP) GO TO 605 606 DIFP = PARDIF 1 IF (PARX0.EQ.0.ODO) GO TO 603 DIF = DSQRT(PARDIF/PARX0) IF (DIF.LE.EPSLON) GO TO 604 603 IER = 1 RETURN 604 RETURN 605 IER = 2 RETURN END C SUBROUTINE GXPR(XW,N,DELTA,GXP) IMPLICIT REAL*8(A-H,0-Z) DIMENSION XW(N),DELTA(N) COMMON/CX1/GAP(5),GAW(5),EN1(8,5),EM1(8,5),EM2(8,5),NGAUSS COMMON/C4/W,H,SPAN,PLN,GAMA1,SREF,XKC,XKT COMMON/CX4/NELEM,NBC(21),IX(21,4) CALL COLUMN(XW,N,PU) GXP = PU - PLN*(GAMA1*XW(N-1) + XW(N)) 1 = 0 6644 1 = 1 + 1 XW(I) = XW(I)*1.01D0 CALL COLUMN(XW,N,PU1) XW(I) = XW(I)*0.99D0/1.01D0 CALL COLUMN(XW,N,PU2) XW(I) = XW(I)/0.99D0 DELTA(I) = (PU1 - PU2)/(0.02D0*XW(I)) IF (l.GE.(N-2)) GO TO 1202 GO TO 6644 1202 DELTA(N-1) = -PLN*GAMA1 DELTA(N) = -PLN RETURN END SAMPLE INPUT/OUTPUT FILE 1 20 SAMPLE INPUT DATA FILE FOR RELIABILITY ANALYSIS 15870.0 0.038 0.089 1.0 1 .25 U 5 9660000.0 2.0 10.0 5.0 4 2 3 1 2 1 3 5 1 3 5 3 3 3 1 1 0 0 0 0 0 0.01 10 3514.0 6738.0 3.97 0.0 33.845 7.8559 4.03 29.861 2.9111 1.0 0.15 0.75 0.15 4538.4 7.036 8.358 1 .025 0 .881 3.2 0.6 SAMPLE OUTPUT FILE CODES : 3 3 BETA = 5.136 ITERATIONS = 4 TOLB = 0.0100 VECTOR XO VECTOR X SENSITIVITY COEFFS. L = 3.2 </>p = 0. 6 4693.6 7.053 8.538 1.025 0.881 3878.8 32.302 30.358 1.3053 1.0553 0.8282 0.0000 0.0000 0.3963 0.3963 EXPLANATIONS Vector Xo : I n i t i a l ( t r i a l ) value f o r the v a r i a b l e s . V e ctor X : Coordinates of the most l i k e l y f a i l u r e point (design point) S e n s i t i v i t y c o e f f i c i e n t s S e n s i t i v i t y of (3 to each of the v a r i a b l e s . In t h i s case j3 i s most s e n s i t i v e to X ( 1 ) , X(3, and X(5). I t i s not s e n s i t i v e to X(2) and X(3). S o l u t i o n corresponding to L = 3. 2m, 4>p = 0.6 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0062671/manifest

Comment

Related Items