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: F I N I T E ELEMENT AND R E L I A B I L I T Y 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 E n g i n e e r i n g  We accept t h i s t h e s i s as conforming to the r e q u i r e d standard  THE UNIVERSITY OF BRITISH COLUMBIA October  1987  © EXAUD NOE KOKA , 1987  In p r e s e n t i n g this thesis in p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at The U n i v e r s i t y of B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I further agree that permission f o r e x t e n s i v e copying of t h i s t h e s i s f o r scholarly purposes may be granted by the Head of my Department or by h i s or her r e p r e s e n t a t i v e s . It is understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain shall not be allowed without my written permission.  Department of C i v i l  Engineering  The U n i v e r s i t y of B r i t i s h Columbia 2075 Wesbrook P l a c e Vancouver, Canada V6T 1W5 Date: October  1987  ABSTRACT  This  thesis  describes  the  consists analysis  of two and  parts.  The  implementation  first  of  a  part finite  element computer model f o r the general p r e d i c t i o n of f a i l u r e of wood members i n bending or i n combined bending and compression. Both i n s t a b i l i t y and m a t e r i a l s t r e n g t h are  included.  a n a l y t i c a l and  The  program  test  is verified  results.  A good  using  axial  failures available  agreement  with  the  r e s u l t s p r e d i c t e d by t h i s program i s observed. The  second part d e s c r i b e s a procedure f o r the s t r u c t u r a l  reliability  evaluation  of  a  compression  member  random loads and m a t e r i a l  v a r i a b l e s . The program  here f o r the r e l i a b i l i t y  study  program  and  the  Rackwitz-Fiessler  c a l c u l a t i o n of the r e l i a b i l i t y failure  function,  Rackwitz-Fiessler  l i n k s the  which  the f i n i t e element r o u t i n e .  developed  finite  algorithm  element for  index /3. The g r a d i e n t of  is a  algorithm,  assuming  necessary  input  to  using  reliability  study f o r a t y p i c a l column problem are compared a g a i n s t a v a i l a b l e r e s u l t s obtained [as  outlined  slenderness  in  by f o l l o w i n g the code  CAN3-086.1-M84  (1984)]  accurate  the  procedures  for  different  ratios.  A performance f a c t o r <p = 0.75, f o r compression of any  the the  i s computed n u m e r i c a l l y The r e s u l t s of the  the  length i s  recommended  and c o n s i s t e n t l e v e l of ii  i n order  to obtain  reliability  i n the  members a  more design  p r o c e s s . I t i s estimated adopted  in  reliability  the  current  that i f t h i s f a c t o r design  practices,  index of the order of 4.0  can be  <p^ = a  0.75  level  achieved.  TABLE OF CONTENTS ABSTRACT  i i  LIST OF FIGURES  vi  LIST OF TABLES  viii  NOTATION  ix  Acknowledgment  xi  1. INTRODUCTION 1.1. Background 1.2. O b j e c t i v e s 1.3. T h e s i s O r g a n i s a t i o n  1 1 3 3  2. CURRENT CODE REQUIREMENTS AND PREVIOUS RESEARCH WORK .. 5 2.1. I n t r o d u c t i o n 5 2.2. Current Code Requirements 5 2.2.1. C o n c e n t r i c Compression 5 2.2.2. Combined bending and Compression 7 2.3. Previous Research on Wood compression members . 8 3. FINITE 3.1. 3.2. 3.3. 3.4.  3.5. 3.6. 3.7. 3.8. 3.9.  ELEMENT ANALYSIS Introduction Assumptions Stress strain relationship F i n i t e Element Approximation 3.4.1. I n t r o d u c t i o n 3.4.2. Kinematic Assumptions 3.4.3. Problem f o r m u l a t i o n 3.4.4. I n t e r p o l a t i o n Functions 3.4.5. S t r a i n Displacement R e l a t i o n s 3.4.6. V i r t u a l Work Equations 3.4.7. Newton-Raphson Method 3.4.8. Computation procedure 3.4.9. Numerical I n t e g r a t i o n Convergence c r i t e r i o n f o r s o l u t i o n v e c t o r O b t a i n i n g the u l t i m a t e load Pmax F a i l u r e C r i t e r i o n and S i z e e f f e c t s 3.7.1. S i z e e f f e c t i n compression 3.7.2. S i z e e f f e c t i n t e n s i o n Program S t r u c t u r e Discussion  4. PROGRAM VERIFICATION 4.1. I n t r o d u c t i o n 4.2. Comparison of r e s u l t s 4.2.1. P r e s e n t a t i o n of R e s u l t s iv  ....  10 10 10 11 17 17 17 20 22 26 27 28 31 33 34 35 39 40 41 44 44 46 46 46 47  4.3. D i s c u s s i o n  53  5. APPLICATIONS 5.1. I n t r o d u c t i o n 5.2. Numerical example 5.3. Observations  55 55 56 58  6. RELIABILITY ANALYSIS 6.1. I n t r o d u c t i o n 6.2. The 0 method f o r R e l i a b i l i t y a n a l y s i s 6.2.1. R a c k w i t z - F i e s s l e r A l g o r i t h m 6.3. Problem Formulation 6.3.1. Code Design Equation 6.4. The G f u n c t i o n f o r the problem 6.5. The r e a l i a b i l i t y Program 6.6. R e l i a b i l i t y R e s u l t s 6.6.1. D i s c u s s i o n of r e s u l t s  59 59 62 64 65 67 67 69 71 75  7. C o n c l u s i o n s and Recommendations 7.1. C o n c l u s i o n s 7.2. Recommendations  76 76 77  REFERENCES  79  APPENDIX A  81  v  LIST OF FIGURES 1. A x i a l  load-slenderness  relationship  for concentric  loading  2. S t r e s s s t r a i n assumptions f o r wood 3. B i l i n e a r s t r e s s s t r a i n  relationship  4. D i s t r i b u t i o n of s t r e s s e s 5. S t r e s s - s t r a i n  13 f o r wood  and s t r a i n s  relationship  7  f o r various m  6. Large deformation of a beam element  14 15 17 18  7. Large deformation of a x i s of beam element 8. i k f i n i t e element i n the x-coordinate system  21  9. Flow chart f o r o b t a i n i n g the s o l u t i o n vector {X}  33  19  fc  10.  E s t i m a t i n g the i n i t i a l f a i l u r e load  Pi  36  11.  I t e r a t i o n process f o r o b t a i n i n g  12.  S t r e s s p r o f i l e across the s e c t i o n  13.  Comparison of program (with m=-l) and a n a l y t i c a l r e s u l t s [3] A x i a l load-Slenderness p l o t s f o r the data of Table 2,(no s i z e e f f e c t s i n c l u d e d )  14.  Pmax  15(a). A x i a l load-Slenderness curve with s i z e e f f e c t taken i n t o account, e = 2mm 15(b). A x i a l l o a d - s l e n d e r n e s s c u r v e s with v a r y i n g k f o r a constant k c  16.  Some l o a d i n g  fc  cases and support c o n d i t i o n s  17(a). U l t i m a t e s t r e n g t h i n t e r a c t i o n curves f o r simply supported columns s u b j e c t e d to uniformly d i s t r i b u t e d load  vi  38 42 48 51 51 52 55  57  17(b). Non-dimensionalized u l t i m a t e i n t e r a c t i o n curves of F i g u r e 18.  Typical  19.  Probability for  the  17(a)  problem f o r r e l i a b i l i t y density  variable  strength  evaluation  59  function  G  61  20.  D e f i n i t i o n of the  21.  R e l i a b i l i t y r e s u l t s as a f u n c t i o n of <j>^', s i z e e f f e c t s not i n c l u d e d R e l i a b i l i t y r e s u l t s with 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  22.  58  reliability  vii  index 0  63 74 75  L I S T 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  viii  NOTATION  The f o l l o w i n g  symbols a r e f r e q u e n t l y u s e d  P = axial  compressive  NDS = N a t i o n a l D e s i g n a = normal  load Specifications  m = s l o p e of the s t r e s s e = u  strain  falling  branch  deformation  B = width  of c r o s s  section  H = depth  of c r o s s  section  = Length  o f member  TJ = n o r m a l i z e d c o o r d i n a t e s  {6}  = nodal displacement  M, M,, M  2  = Shape  N, N, = Shape K  curve,  deformation  w = transverse  £,  ( i n Mpa)  strain  = axial  L  (USA).  area  = mean modulus o f e l a s t i c i t y  Q  report  stress  A = cross sectional E  in this  = global  T  vector  f u n c t i o n s f o r the w  functions f o r the u displacements  tangent  stiffness  = initial  {A5}  = incremental displacement  AP = l o a d  solution  {X} = New  solution  k  effect  t  = size  vector vector  increment  {X } = i n i t i a l Q  displacement  matrix  {6 } Q  displacements  vector  vector  factor  i n tension  ix  k  c  = s i z e e f f e c t f a c t o r i n compression  G = failure  function  Pf = p r 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  = f a c t o r e d compressive r e s i s t a n c e p a r a l l e l to g r a i n 0p = r e s i s t a n c e K  c  f a c t o r i n pure comression  = Slenderness f a c t o r  d = dead load v a r i a b l e / = l i v e load  variable  Design Load F a c t o r s a  D  = dead load f a c t o r  a  L  = l i v e load f a c t o r  7! = r a t i o of nominal dead load to nominal l i v e  Subscripts c = compression t = tension P  DN  =  n  P  LN  =  n  o  i-  m  o  m  ^  n  n  l  dead  load  l  live  load  a  a  f = failure  x  load  ACKNOWLEDGMENT I wish to express my Professor  0.  R.  deepest g r a t i t u d e to my  Foschi  for  his  invaluable  guidance throughout the research work and of t h i s  are a l s o  and  due  support.  to P r o f e s s o r  The  Commonwealth S c h o l a r s h i p and a d d i t i o n a l funding of  the  advice  and  preparation  thesis.  Thanks advice  i n the  supervisor,  University  financial  B. Madsen support  his  from  the  F e l l o w s h i p Committee and  some  from the department of C i v i l of  for  British  Columbia  is  Engineering gratefully  acknowledged. Lastly,  thanks are due  i n t h i s p r o j e c t , d i r e c t l y or  to those who indirectly.  xi  have rendered  help  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 s u b j e c t e d  to  lateral  loads  occur very f r e q u e n t l y , such as i n b u i l d i n g frames, b r i d g e or roof t r u s s e s They are  and  usually  stress c r i t e r i o n  other important proportioned  engineering  to  satisfy  structures.  some  limiting  set by design s p e c i f i c a t i o n s or codes.  s t r e s s e s developed  a t any  cross  section in  such  The  members  consist of: 1.  the a x i a l s t r e s s caused by the compressive  2.  the primary bending s t r e s s due t o the l a t e r a l l o a d s , and  3.  the  secondary  amplification compressive The for  bending of  the  stress  forces ,  resulting  deflections  produced  compressive  the  by  the  forces .  secondary bending s t r e s s becomes p a r t i c u l a r l y members  from  with  a  high  forces.  The  slenderness procedures  secondary s t r e s s e s i n e l a s t i c columns l i t e r a t u r e on s t a b i l i t y  ratio  important and  large  f o r computing  the  are d e s c r i b e d i n the  [1].  theory  Although e l a s t i c a n a l y s i s i s used e x t e n s i v e l y i n design computations,  i t does not g i v e an a c c u r a t e i n d i c a t i o n of the  true l o a d - c a r r y i n g c a p a c i t y , p a r t i c u l a r l y  f o r columns  are not  columns  fail  by  very s l e n d e r . excessive  L a t e r a l l y loaded  bending  after  1  the  stresses  which  generally in  some  2 p o r t i o n s of the member exceed a maximum v a l u e . To the  ultimate strength  perform  a  of such columns,  stability  analysis  e l a s t o - p l a s t i c behaviour design  codes  approach,  and  which  of  consists  i t i s necessary  that  considers  the m a t e r i a l .  specifications of  use  assuming  m a t e r i a l with a maximum normal s t r e s s  determine  Most the  a  to the  available traditional  linear  elastic  failure criterion .  Previous a n a l y t i c a l and experimental s t u d i e s on wood, as r e p o r t e d i n the l i t e r a t u r e [5,6,7], have shown that 1.  wood has  a  non-linear stress  compression, e . g . b i l i n e a r  strain  :  relationship  in  elasto-plastic  relationship,  characteristic contributes  significantly  and 2.  this material  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  still  some  problems  which  remain  unsolved: 1.  The codes do not resulting  2.  give guidance f o r c a l c u l a t i n g  from beam-column  moments  deflections.  An account f o r p o s s i b i l i t i e s of d u c t i l e y i e l d i n g  i n the  compression zone or t e n s i o n f a i l u r e i n the t e n s i o n i s not g i v e n .  zone  3  1.2.  OBJECTIVES T h i s study  i s aimed at a c h i e v i n g three main  objectives,  namely : 1.  To develop  a finite  element a n a l y s i s  f o r the  p r e d i c t i o n of the f a i l u r e of a compression t r a n s v e r s e l o a d s . The the  analysis will  non-linearities  due  to  general  member  under  take i n t o  account  slenderness  effects  (geometric),a n o n - 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 the m a t e r i a l , and e s t i m a t i o n of f a i l u r e l o a d by e i t h e r t e n s i o n or 2.  controlled  compression.  The a n a l y s i s w i l l be implemented i n a computer The  computer  accomodating  program various  for  will  allow  support  program.  flexibility  conditions  and  in load  conf i g u r a t i o n s . 3.  To e v a l u a t e the r e l i a b i l i t y of wood compression  members  assuming random loads and m a t e r i a l v a r i a b l e s .  1.3. THESIS ORGANISATION 2  Part  provides  recommendations and members. Part  3  f i n i t e element Part 4  by  summary  of  current  previous research  describes a  analysis  provides  developed  a  on wood  general  and the  a verification  the  i n Part 3 , using experimental  previous  researchers  [6,7],  Part  code  compression  formulation  computer of  design  of  the  implementation. computer  r e s u l t s as 5  program reported  presents  the  4 application  of  the  analysis  to  compression member, where a x i a l  a  laterally  load versus  transverse  i n t e r a c t i o n diagrams a r e developed f o r d i f f e r e n t r a t i o s using a 2x4-in  a  computer  reliability  index  program 0  of  f o r the a  wood  algorithm.  obtained  from t h i s study  the end  of  the chapter.  general c o n c l u s i o n of for  slenderness  A  evaluation.  evaluation  compression  c o n s t r u c t e d , using the program developed Rackwitz-Fiessler  load  SPF c r o s s s e c t i o n .  Part 6 d i s c u s s e s the concept of r e l i a b i l i t y Here,  loaded  of  the  member  i n Part 3 and the  summary  of  the  results  f o r a s p e c i f i c problem i s given And  l a s t l y , Part  the r e p o r t  f u r t h e r research and study.  is  and some  7  provides  at a  recommendations  2. CURRENT CODE REQUIREMENTS AND PREVIOUS RESEARCH WORK  2.1.  INTRODUCTION The  failure  depend on i t s  c h a r a c t e r i s t i c s of  a  s l e n d e r n e s s . The u l t i m a t e  compression members i s d i r e c t l y  compression c a p a c i t y of  observed.  Thus,  characteristic of intermediate two  a change t o  types  of  a  instability  l e n g t h , there i s a t r a n s i t i o n failure  regimes,  in  which  failure  both the compression  is is  member  between case  of  length  of f a i l u r e  of slender compression members. For a  c a p a c i t y depends on stiffness  a b u c k l i n g type  lateral  short  r e l a t e d t o the s t r e n g t h  the m a t e r i a l i n compression. With an i n c r e a s e i n the of the member,  member  the  these load  s t r e n g t h and the  of the m a t e r i a l .  2.2. CURRENT CODE REQUIREMENTS  2.2.1. C o n c e n t r i c C o m p r e s s i o n Current  design codes  s h o r t , intermediate slenderness  or  c l a s s i f y compression members slender members  according  to  r a t i o C . For r e c t a n g u l a r c r o s s - s e c t i o n s , c  into their  6 where L = l e n g t h o f t h e member d = dimension in  the d i r e c t i o n  Thus, Short slenderness crushing is to  o f t h e c r o s s - s e c t i o n o f t h e member  members a r e  ratios  parallel  F ,  performance  considered  o f 10 o r l e s s .  t o be  They w i l l  their  factor  <f>  strength  cross  those  normally  to the g r a i n . Their design  b a s e d on t h e s p e c i f i e d grain,  of b u c k l i n g .  fail  allowable  i n compression  sectional  with  area  by  load  parallel A,  and  a  ,  Cr  P  — - F A c  Slender relation  members  0 P  generally  a n d have s l e n d e r n e s s  dependent  (2)  on t h e mean e l a s t i c  follow ratios  the  Euler  exceeding  modulus E  Q  buckling  C^, a  and t h e  number  specified  compression  s t r e n g t h , F , o f t h e column m a t e r i a l . The number  C.  by :  i s given  where E  c  i s called  members a n d between  t h e modulus o f e l a s t i c i t y  i s equal  t o 0.74E . Q  f o r compression  F o r lumber,  can  20 t o 25, d e p e n d i n g on t h e g r a d e o f t h e member  vary under  considerat ion. Intermediate members have s l e n d e r n e s s r a t i o s between and C^.  They  are designed  s t r e n g t h which e m p i r i c a l l y  using  a  10  modified  compression  i n t e r p o l a t e s between  slenderness  r a t i o s of 10 and C^. The above c l a s s i f i c a t i o n  i s illustrated  - — Crushing  Euler  Figure  1. A x i a l l o a d - s l e n d e r n e s s r e l a t i o n s h i p f o r c o n c e n t r i c loading.  2.2.2.  Combined  b e n d i n g and  Compression  A compression member i s o f t e n s u b j e c t e d to bending about e i t h e r one  or both  axes, and  the combined  e f f e c t of the  bending and a x i a l l o a d i n g must be c o n s i d e r e d . For t h i s of l o a d i n g ,  most c u r r e n t  c r i t e r i o n based on  codes  a linear  specify a  interaction  simple  type  failure  between the  axial  8 load c a p a c i t y  of a  moment c a p a c i t y may  concentrically  in bending  be a p p l i c a b l e  the  attributed,  as long as the  in  approach  l i t t l e has  been  the  done so f a r  to  l i n e a r range. T h i s may  to the u n c e r t a i n i t i e s about the  form of the c u r v i l i n e a r  the  wood member remains i n  behaviour beyond the in part,  and  alone. Therefore, t h i s  l i n e a r e l a s t i c range. Very predict  loaded column  be  precise  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  compression.  2.3.  PREVIOUS RESEARCH ON WOOD COMPRESSION MEMBERS Most p r e v i o u s s t u d i e s  (Newlin and  Trayer  be a l i n e a r  tested  Wood 1950)  behaviour. They  to used  predict  when a  with  verify their a second  Theilgaard  combined theory  order  d e f l e c t i o n s of  axial  for  linear  to  limiting  Thus, Larsen and  members  the  beam-columns  have c o n s i d e r e d wood  which f a i l s  i s reached.  wood  t r a n s v e r s e loads  equation to  wood columns and  e l a s t i c material  compression s t r e s s (1979)  1925;  on  and  beam-column differential  e l a s t i c beams  and  beam-columns. Bleau  (1983) and  Buchanan (1984), conducted an  extensive  j o i n t experimental study on e c c e n t r i c a l l y loaded columns c a l i b r a t e and able to  v e r i f y t h e i r strength  predict  r e s u l t s of a x i a l members. Buchanan  the  strength  tension used a  and  of  models. T h e i r full  size  compression  mean modulus  to  models are  lumber,  t e s t s on  using similar  of e l a s t i c i t y ,  E , Q  9 equal  to  10000 Mpa  Bleau's data his  t o c a l i b r a t e h i s m o d e l . Zann ( 1 9 8 5 ) ,  (1983),  with E  equal  Q  s t r e n g t h m o d e l . Zann's m o d e l  (1975) d e s i g n account of  the  test  results i s the  different  and  corresponding r e p o r t s an E The the  The  is  used, '  to  the  mean  o f 9660 Mpa  discussed  a more g e n e r a l  method  of  implementation  test  be  the  solution  formulation will  the  a c c o u n t of  biaxial  i n both  fact that in this  E„ o  and  discussed  calibrate  i s b a s e d on  is  results.  f o r the p o p u l a t i o n  by  to  cases  i s reported,  model d e v e l o p e d i n t h i s r e p o r t  ideas  provides  Q  takes  Although  which remains unanswered o  10400 Mpa  (1985),  r e c o m m e n d a t i o n s , and  b i a x i a l bending.  agreement w i t h  to  used  a  NDS  good  question  each case not  the  Bleau  a one  (1983),  tested.  i n c o r p o r a t e s some  previous  researchers,  of and  t o t h e beam-column p r o b l e m . the  corresponding  i n Part  3 of t h i s  computer report.  3. F I N I T E ELEMENT ANALYSIS  3.1.  INTRODUCTION This  chapter  describes  the formulation  of  a  finite  element a n a l y s i s f o r p r e d i c t i n g 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 l o a d s . The theory and assumptions i n t h i s chapter the b a s i s of a computer  w i l l be d e s c r i b e d along  program developed  to implement  with the  model.  3.2. ASSUMPTIONS The  f o l l o w i n g assumptions a r e made:  1.  plane s e c t i o n s remain  2.  the s t r e s s - s t r a i n law f o r the m a t e r i a l i s known.  3.  m a t e r i a l p r o p e r t i e s are constant  plane.  along the l e n g t h of the  member. 4.  bending i n only one plane  i s considered.  5.  no  of  torsional  or  out  plane  deformations  considered. 6.  d u r a t i o n of l o a d e f f e c t s are not c o n s i d e r e d .  7.  shear deformations  are s m a l l , hence n e g l e c t e d .  10  are  11 3.3. STRESS STRAIN Various  RELATIONSHIP  studies  stress-strain  [4,5,8]  have  behaviour of wood  focussed  i n compression  g r a i n , with the aim of d e r i v i n g a mathematical to represent t h i s  behaviour. R e c e n t l y ,  (1970), proposed a s t r e s s s t r a i n  c  the  parallel  to  relationship  Malhotra and  relation  1 e = — [ca-(l-c)f E o  on  Mazur  of the form :  a lnd- — )] f c  (4)  where: e = strain. a = stress f  c  = maximum compression  E  Q  = mean modulus of e l a s t i c i t y  c = shape For c = 0.99,. the  stress  parameter.  curve d e s c r i b e d by Equation (4) i s shown  as A i n F i g u r e 2. A mathematical c l e a r dry wood  equation f o r the s t r e s s s t r a i n curve  i n compression  at  v a r i o u s g r a i n angles  was  The equation takes  the  a l s o developed by O H a l l o r a n ( 1 9 7 3 ) . 1  following  for  form  a = E e - he  (5)  n  Q  Where a, e, and  E  are d e f i n e d above  and A,n a r e  equation  12 constants determined  by f i t t i n g the equation t o a given  of experimental data. A curve B i n  F i g u r e 2.  p l o t of t h i s  equation i s shown  T h i s equation cannot  maximum s t r a i n because i t may  set  be used  as  beyond  take on negative values  very  rapidly. A comprehensive study on the s t r e s s - s t r a i n of timber with  d e f e c t s , i n compression  has been  by  done  Glos  (1978),  as  relationship  p a r a l l e l to  reported  by  grain, Buchanan  (1984). Based on experimental data, the curve shown as C Figure 2  was  proposed.  T h i s curve  number of m a t e r i a l parameters  i s characterised  in  by  a  that depend on measurable wood  p r o p e r t i e s , namely d e n s i t y , moisture content, knot r a t i o  and  the percentage  for  of compression  modelling purposes these  necessarily  wood. Using  t h i s curve  i n v o l v e s the c a l i b r a t i o n  of  parameters.  A simple b i l i n e a r p r o p o s a l by Bazan (1980), as d i s c u s s e d by Buchanan (1984), appears t h i s p r o p o s a l , i t i s assumed branch  to be the  most recent one.  that the slope of the  i s a v a r i a b l e which can be a r b i t r a r i l y  falling  taken as  value which produces maximum bending moment f o r any  In  that  neutral  a x i s depth. A p l o t of t h i s curve i s shown as D i n F i g u r e 2.  13  F i g u r e 2. S t r e s s s t r a i n assumptions f o r wood The a n a l y s i s  in this  curve D, with the  study  uses the  exception that the  branch of the s t r e s s - s t r a i n  relation  simple  slope of the  e l a s t i c and a  in Figure  2 are c h a r a c t e r i s e d  non-linear part. Therefore  g e n e r a l l y be expressed  falling  i s c o n s i d e r e d to be  m a t e r i a l p r o p e r t y , i n agreement with Buchanan The curves  bilinear  a  (1984). by a  linear  the s t r e s s e s  can  as  a = E e Q  + F(e)  (6)  14 The is  as  stress-strain relationship  shown  behaviour  in  in  Figure  tension,  3  and  with  a  adopted i n t h i s  includes bilinear  linear  study elastic  relationship  in  compression and a f a l l i n g branch a f t e r maximum s t r e s s .  F i g u r e 3. B i l i n e a r s t r e s s  Using the  strain relationship  above s t r e s s - s t r a i n  d i s t r i b u t i o n of s t r e s s e s i s as shown i n F i g u r e 4.  relationship,  and s t r a i n s  f o r wood  the  resulting  in a rectangular  beam  T  r  section  strains  stresses  F i g u r e 4. D i s t r i b u t i o n of s t r e s s e s and s t r a i n s The curve i n F i g u r e 3  can be mathematically represented  by  the f o l l o w i n g e x p r e s s i o n s :  For segment  1-2 ;  o = ~\£ \ c  ~ |f |m - mE e c  (7)  c  For segment 2-3 ;  (v  ° = V  Or,  i n combination,  a = E e - [ E e + |f |(1+m) + m E e ] ( l Q  Q  c  Q  - A(e + | e | ) ) c  (9)  16 where A(e+|e  j) i s the step f u n c t i o n d e f i n e d as f o l l o w s :  if  e > - |e | ; A = 1  if  e < - |e | ; A = 0  (10)  c  c  Hence, f o r 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 = 0; and Equation  a = Ee  For  the e l a s t i c  for  tension  m  <  m  (9) reduces to  - {Ee  +  |f |}  (1-A(e  +  |e  I))  (11)  case, m = -1, and we have Equation  and compression.  illustrated  behaviour,  This explanation  (9)  is  both  further  i n F i g u r e 5 below.  o m = - l  F i g u r e 5. S t r e s s - s t r a i n r e l a t i o n s h i p f o r v a r i o u s m.  17  3.4. F I N I T E ELEMENT APPROXIMATION  3.4.1.  Introduction  The  finite  element  v e r s a t i l e technique s o l u t i o n of  the  is a  very  presently available  problems  advantages of  method  of  the  method  type  have  powerful  f o r the  considered  been  a p p l i c a t i o n s e x t e n s i v e l y demonstrated  and  numerical here.  recognised  and i t s in  steel  and c o n c r e t e s t r u c t u r e s , and f o r some wood s t r u c t u r e s  such  as wood f l o o r s ,  wood diaphragms and  a p p l i c a t i o n of the method to  particularly  The  t r u s s e s . However, the  wood beam-column a n a l y s i s  has  not been explored to an e q u i v a l e n t degree.  3.4.2. K i n e m a t i c A s s u m p t i o n s As  the  displacements  non-linearity beam-column.  is  introduced  Consider  deformations but small F i g u r e 6,  u  become  and w  a  in  beam  large,  a  geometric  the  deformation  element  undergoing  s t r a i n s . For the  of  large  geometry shown  are, respectively, a x i a l  a  and  in  lateral  displacements of the c e n t r e l i n e of the beam. A and 0 a r e two p o i n t s on  the  same  plane  such  that O  i s on  centerline  (axis) and A i s at a d i s t a n c e z from O  the  (positive  z ) . L i n e OA represents c o n d i t i o n s before deformation, l i n e O'A  1  represents conditions a f t e r  beam  while  deformation.  Assuming that plane s e c t i o n s remain plane, the  rotation  18  dw of the c r o s s - s e c t i o n i s 6 = — . Figure dx A and B, at the same  distance  7 shows two  z from the c e n t r e l i n e .  deformation, these p o i n t s a r e at A* and B'. ^  Figure  6.  beam  Large deformation of a beam element  points, After  19  F i g u r e 7. Large deformation  of a x i s of beam element.  From the geometry of F i g u r e 7 i t f o l l o w s that  ds  2  = dx  2  [(1  +  du.  —&  dx  If the expression above i s higher order expression  terms are  2  +  dw.  (  — A  dx  )2]  expanded b i n o m i a l l y , and  n e g l e c t e d , the f o l l o w i n g  if  the  simplified  i s obtained. du, ds  = dxd  + —&  1 + -  dx Therefore,  )  the corresponding  2  dw  x  ( —£  )  2  +  )  (13)  dx  strain  e.  at a d i s t a n c e 2 i s  20 du, A dx  A  j .  1 dw, _ / A \ 2 dx  2  dw But, from F i g u r e 6, u. = u - z — dx  (15)  A  thus, du, —  du  dw 2  (16)  - - - - z — -  A  dx  dx  dx'  A l s o , n e g l e c t i n g higher order terms,  w  Thus combining  A  = w  (17)  Equations (14),(16)  and  (17), we  get the  s t r a i n at a height z as :  e=  3.4.3. A  Problem beam  du 1 dw — + - ( — ) dx 2 dx  dw 2  2  - z  dx  (18)  2  formulation  element with  two  end  nodes i s used  i n the  f o r m u l a t i o n . Let us choose a l o c a l c o o r d i n a t e li (-1 < ^ < 1) in each element  such as the  one shown i n  Thus, along the  x axis coordinate  Figure 8  below.  system each element  two end nodes, i and j separated by a l e n g t h 2A.  has  21  £=0  £=-1  x=o  u du dx w dw dx  F i g u r e 8. i * " *  1  u du. dx w dw dx  f i n i t e element  Thus, the x-coordinate of can be expressed as x = x of freedom v e c t o r are 4  degrees  i n the x - c o o r d i n a t e system.  any p o i n t w i t h i n the  element  + A£. The elemental nodal  c  degree  i s represented as in Equation (19).  of freedom  at  each node,  namely  dispacements u and w and t h e i r r e s p e c t i v e f i r s t  the  There two  derivatives.  i  u  du (—  -  )  dx  •  1  dw (  )•  —  dx {5}  1  =  (19) U  j  du (  )• 3  —  dx w.  D  (  dw —  dx  )•  ^  22 3.4.4. I n t e r p o l a t i o n F u n c t i o n s In t h i s study, complete approximate is  c u b i c i n t e r p o l a t i o n s are used to  the displacements u and w w i t h i n an element.  important  to note that i n order to s a t i s f y  conditions,  the  displacement  i n t e r p o l a t i o n . However, a give an fewer  improved  only  define  requires  cubic i n t e r p o l a t i o n  approximation  of  the  fully  was  the a x i a l  function.  The  describe  displacements  linear used  stress  u(*> < '  -  ( ~ "  -  ~  O  A  and  the  nodes p r o v i d e s u f f i c i e n t a  cubic  polynomial  displacements u and w are thus g i v e n as  V  a  c u b i c i n t e r p o l a t i o n r e q u i r e s 4 parameters  d e r i v a t i v e s at the two  U  compatibility  to with  elements.  A complete  to  u  It  U +  -  2  £ )u.+ 3  1  A  -  8  to first  parameters  function.  The  follows:  (1-*-  du  £ + S )( 2  _  3  )  dx  1  (20)  1 3 1 1 du + ( - + - £ + - £ )u.+ - (-1-£+£ +£ ) ( — )• 2 4 4 ^ 8 dx ^ 3  /- \ W  U )  1 =  (  _  3 -  _  2  +  1 3 ( _ + _ 2 4  4  1 £+  3  1  _  £ 3 )  W  . +  -  dw (I"?"  £ +  £ 3 ) (  2  _  (21) )  4 8 dx 1 1 dw + _ £ )w.+ _ (-1-$+£ +£ ) ( ) 4 J 8 dx J 1  3  €  2  2  3  23 where  In v e c t o r matrix  £ =  2x - 2x 2A  —  n o t a t i o n , we can w r i t e  U = {N} {5} ; — = { N , } dx T  w = {M} {6} T  dw ; — = {M,} dx  where N , N , , M, M, and  T  (22)  {5} dw 2  {5} ;  T  M ,  dx  2  = =  {M } {6} T  2  are v e c t o r f u n c t i o n s of £  2  given  by the f o l l o w i n g e x p r e s s i o n s :  r  1  3  1  - - - £ + 2  A  4  4  £  4  d - £ - £  2  3  +£ ) 3  0.0  {N}  0.0  =  -  (23)  1  3  1  2  4  4  (-1.0-£+£ +£ ) 2  0.0 0.0  3  24 1  3  _  3  (- _ +  A  4  4  4  (-1.0-2£+3£ ) 2  0.0  0.0 1  {N,} =  3  -  3  ( - - -  A  4  1  £ ) 2  4  (24)  (-1.0+2£+3$ ) 2  0.0  0.0  0.0  0.0 1  3  1  - - - £ + 2  {M} =  4  £  3  4  (25)  A (1.0-*-$*+$3)  -  0.0  0.0 1  2  3 +  1  - £- - £ 4  3  4  (-1 . 0 - £ + £ + £ 2  3)  25  0.0 0.0 1 3 3 - ( - + A 4 4 4  i; ) 2  (-1.0-2£+3£ ) 2  0.0 (26)  0.0  1 (!-!«•> A  4  1 4  4  (-1.0+2£+3£ ) 2  0.0 0.0 1  — A  3  ( -  2  2  (-2.0  — 4A  O  + 6U  0.0 0.0 1 3 — ( - - £ ) A 2 1 — (2.0+6£) 4A 2  (27)  26 3.4.5. S t r a i n Displacement For  a  laterally  Relations  loaded  column  problem  with  deformations, the w displacements  w i l l be much l a r g e r  the a x i a l  the  displacements  terms c o n s i d e r e d are  u. Thus  s i m i l a r t o those  strain  large than  displacement  derived i n  Equation  (18) above, where: du dw e = — - z — dx dx  1 dw + - ( — ) 2 dx  2  2  S u b s t i t u t i n g the displacement the r e s u l t s i n symbolic  (28)  2  f u n c t i o n s i n t o Equation  form a r e  e = [B+B(6)]6  where B, r e p r e s e n t s while B(6), includes  which  the  displacement  (28),  (29)  the l i n e a r is a  s t r a i n displacement  function  contribution  of  of the  the  8  terms,  parameters,  non-linear  strain  terms. Thus h B = N, " - z M , where z = - T? ~ 2 1  T  2  1  B(5) = - {8}  T  M,^!  T  (30)  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. D e f i n i n g # as v i r t u a l dislacement of  the nodal  v a r i a b l e s , the  a  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 l i n e a r 95 -  we n e g l e c t i n e r t i a  f u n c t i o n of 8. Thus, i f  f o r c e s , the v i r t u a l work equation  reduces  to / £odv =  (3 2)  v  where P i s the c o n s i s t e n t  load v e c t o r , c a l c u l a t e d using  shape f u n c t i o n s as i n d i c a t e d i n Equations  (23), (24), (25),  (26) and (27). V i s the volume of the member.  Substitution  of the r e s u l t i n g equation f o r 7_ i n t o the Equation to the system  of governing equations  the  (32) leads  f o r an element,  that  is,  J t B + C(6)] adv T  v  Assembling  = P  (33)  the element equations i n the usual f i n i t e element  manner l e a d s t o the g l o b a l system of e q u a t i o n s . In order  to  28  f i n d the s o l u t i o n t o the n o n l i n e a r system of Equations it  i s convenient t o  $( 8), such  solution  $(8)  may  of  8,  that  *(8)  The  i n t r o d u c e the v e c t o r f u n c t i o n  (33),  8  = J [B + C(8)]  now has t o  be  found  T  v  s a t i s f y *(8)  numerically  (34)  adv - P  =  0. The zeros  v i a the  of  Newton-Raphson  procedure as o u t l i n e d below.  3 . 4 . 7 . Newton-Raphson Method T h i s i s a commonly used e q u a t i o n s . The  method  technique t o s o l v e Thus,  at  8+A8,  uses  technique to s o l v e a  first  order  approximation  n o n - l i n e a r equations through the  first  order  non-linear  approximation  iteration. f o r the  f u n c t i o n $ w i l l be d$(8) $(8+AS)  where $ ( 8 )  = $(5)  + [  8{8}  ]A8  (35)  i s a f u n c t i o n of the displacement v e c t o r {8}. The  29  above equation can f u r t h e r be s i m p l i f i e d  into  A6 = - K ~ <I»(8)  (36)  1  T  d{$(8)} =  .  [  ]  1  (37)  <HS)  d{8}  where  [K ]  represents  T  Differentiating  the  p a r t s , we o b t a i n (36) {8}  the  tangent  r i g h t hand  side  from  an i n i t i a l  matrix. (34)  of Equation  a s i m p l i f i e d expression  permits an i t e r a t i v e procedure starting  stiffness  for K .  Equation  T  to determine the  approximation  {8^}.  by  vector  Thus,  in  general,  6  *  i (  " i 5  +1  6  i  )  where the matrix K which f o l l o w s .  =  T  " [K  -M o B  is  +  * V  (38)  (  ] T  B  < i> 6  ] a  i  "  p  obtained as shown i n Equation  (39)  (40),  30  EJ* B C(5)dV  +  T  v  E/ C (5)BdV T  v  E/ C (6) T  V  C(5)dV  E/ B BdV T  v  E (l . 0 + m ) / ( l . 0 o  v  Me  +  |e |))B  BdV  T  c  (40) E (1 . 0 - m ) / ( l . 0 0  v  A(  E (1 . 0 - m ) / ( l .0 0  E  o(  - 0 " m ) / ( l .0  / M,M,adV V  v  +  \e \))B  C(5)dV  T  c  +  v  1  6  I e |))C (6)BdV T  c  Me  + |e |))C (5) T  c  C(5)dV  31 3.4.8. Computation The  Cholesky  procedure  decomposition  of the matrix K  to determine the v e c t o r {A8}from Equation  is  utilized  (38). The  boundary  T  c o n d i t i o n s are f i r s t a p p l i e d to the s t i f f n e s s matrix K to the v e c t o r {$(6)}. The  and  T  boundary c o n d i t i o n codes f o r  this  program are as f o l l o w s  1 = u  2 =  du  —  dx  3 = w  4 =  dw  —  dx  Thus, to enforce a boundary are p l a c e d i n t o the o f f column corresponding  c o n d i t i o n equal to zero,  d i a g o n a l l o c a t i o n s f o r the row  to the s p e c i f i e d  1 i s placed  i n t o the d i a g o n a l  degree of freedom i n the i s decomposed and  {$(5)}.  Each element of the  In a d d i t i o n to t h i s , term of the  [K ]  T  v e c t o r {AS}  {A5}  T  i s obtained.  i s compared a g a i n s t  by the  a  specified  [ K ] matrix. Then the matrix  f i n a l l y a solution  in  freedom  T  i n the returned l o a d v e c t o r  and  degree of freedom  [ K ] , while a zero i s p l a c e d f o r the same degree of  value  zeros  acceptable tolerance  specified  user  to  whether a need to do  more i t e r a t i o n s i s necessary  an  determine in  order  32 to r e f i n e the s o l u t i o n . A summary of the whole procedure given  i n the flow chart  is  i n F i g u r e 9.  GEOMETRY /LOADING AND PROPERTIES INITIALISE  MATERIAL  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}  F i g u r e 9. Flow c h a r t f o r o b t a i n i n g the s o l u t i o n v e c t o r  {X}.  33  3.4.9. Numerical  Integration  Since the volume i n t e g r a l s i n the e x p r e s s i o n f o r K complicated,  it  is  s o l u t i o n s . Hence quadrature  difficult  numerical  to  obtain  integration  scheme has been a p p l i e d due  is  the polynomial  Zienkiewicz appearing  number of Gaussian  (1979),  i n the  the f u n c t i o n . Thus, the term {B} f o u r t h order polynomial squared  i n the e x p r e s s i o n  order polynomial Knowing that  term  a  k  in  point  scheme  is  the  of  determines  the  time  T  needed i n  integrate  (40) c o n t a i n s a  f o r [ K ] . T h e r e f o r e , the the i n t e g r a l s  in  +1.  to a c c u r a t e l y  Gaussian  Gaussian  maximum order  at the same  e x a c t l y a (2k-l) order polynomial, Gaussian  to  i n Equation  i n £ and  form  to i t s s u i t a b i l i t y  integral  p o i n t s necessary  closed used.  the l o c a l c o o r d i n a t e system v a r y i n g from -1 According to  are  T  is  scheme  of  {B}  is  highest order  will  8.  integrate  i t f o l l o w s t h a t a 5-point  the  numerical  integration.  Thus, the i n t e g r a l s over the volume V become  l  v  =  2BHA r\ ~— / ^ " 4 J_ J_ d  d r i  A  v  x  where N = 5 was Gaussian  =  l  chosen, and w^  weights.  2BHA N N -7-££ T '3 4 I=I j=j K  { i  ) w  i -i  and w_j are the  w  ( 4 1 )  J  corresponding  34  3.5. CONVERGENCE CRITERION FOR SOLUTION VECTOR The  convergence  of the  step i s checked by the 3  Q  be  the previous  s o l u t i o n vector  at every  E u c l i d e a n norm c r i t e r i o n .  solution vector  and 3  I f we  be the  load let  present  s o l u t i o n , then, as shown below,  l e t Ax  = |d  -  d | represent Q  lengths of d and d . We can then  l ol 3  |3|  where X ( i ) Q  respect i v e l y .  and  2  2  =  X  o  (  i  the d i f f e r e n c e  between  the  write  )  = X (i) 2  X ( i ) are  the components  of  d  Q  and  d,  35  Then, Ax =  The convergence  (x(i) - x ( i ) )  (42)  2  o  c r i t e r i o n based  on the E u c l i d e a n  norm i s  d e f i n e d as Ax  < specified  |3 |  (43)  tolerance  0  3.6. OBTAINING THE ULTIMATE LOAD PMAX The  failure  load  procedure. For f a s t  Pmax i s obtained  convergence t o the  by  an  iterative  s o l u t i o n Pmax, the  f o l l o w i n g approach f o r e s t i m a t i n g an  i n i t i a l guess f o r the  f a i l u r e l o a d i s chosen. F i r s t of a l l ,  the c r u s h i n g  P  c  as w e l l as the E u l e r b u c k l i n g  computed. Regardless of length, the u l t i m a t e value between  P„ and c region of F i g u r e 10.  load Per of the member a r e  the support  load w i l l Per and  strength  c o n d i t i o n s and  be l e s s will  than the  member smallest  l i e w i t h i n the  shaded  36  Figure  10. E s t i m a t i n g the i n i t i a l  The minimum initial solution  failure  of P c load  and P^.  Per i s As  f a i l u r e load P i .  then  a first  taken to step,  the  l e t the  l i e between two load v a l u e s , namely P1=0 and P2=P^.  The  average load P3=(P1+P2)/2 becomes the f i r s t  The  f i n i t e element s o l u t i o n i s o b t a i n e d  occurs,  we  be  i t means  that the  P=P1 and P=P3. Therefore loads f o r next  iteration  (P1+P2)/2 i s c a l c u l a t e d  trial  f o r P=P3. I f f a i l u r e  s o l u t i o n i s between the  we set the minimum and the as  P1=P1 and  and  r e - r u n . I f the member s u r v i v e s ,  the  load.  P2=P3. A  finite  values maximum  new P3  element  program  i t means that the  solution  37 i s now between the v a l u e s the minimum  and the  P1=P3 and  P2=P2. T h i s  until  acceptable  an  P=P3 and P=P2. T h e r e f o r e , we set  maximum loads process tolerance  s u c c e s s i v e estimates of Pmax. I f as TOLP, the i t e r a t i o n s  TOP =  P2-P3 P3  user.  i s repeated i s reached  iteration several between  as  times two  this tolerance i s defined  a r e stopped when  <  The process i s summarized where TOLP i s  f o r next  TOLP  i n the flow  the a l l o w a b l e t o l e r a n c e  c h a r t of F i g u r e  11,  normally set by the  38  Figure  11.  I t e r a t i o n process f o r o b t a i n i n g Pmax.  39  3.7.  FAILURE CRITERION AND SIZE EFFECTS Due  t o the b r i t t l e  observed i n wood  fracture  phenomenon which i s commonly  members, i t may  be important  that  the  a s s o c i a t e d s i z e e f f e c t s be i n c o r p o r a t e d i n the a n a l y s i s .  In  a  is  brittle  material,  a  decrease  i n member  strength  normally observed as a r e s u l t of a c o r r e s p o n d i n g i n c r e a s e i n member s i z e . I f no s i z e e f f e c t s a r e c o n s i d e r e d , the  failure  criterion is  "max " t  <  F  where c  m a x  i s the maximum t e n s i l e s t r e s s  c r i t e r i o n , although strengths differences size  between  simple, does pure  tension  are accountable  i n the member. T h i s  not r e s u l t and  4 4 )  pure  through the  in  different  bending.  Such  incorporation  of  w i l l be a p p l i e d  to  effects. W e i b u l l ' s theory of b r i t t l e  i n c o r p o r a t e the s i z e e f f e c t  phenomenon. Thus, f o r a  of volume V, f a i l u r e i s r e l a t e d  I = J o dv k  fracture  member  t o the i n t e g r a l  (45)  40  with a c o r r e s p o n d i n g f a i l u r e c r i t e r i o n given  I = (a*)  by  (46)  k  where a = stresses * a  = strength  uniform  i n the  member,  of a u n i t volume under  stress,  k = size effect factor, V = volume of the  stresses  domain.  3.7.1. S i z e e f f e c t in compression  t  The  parameter  compression  (buckling  |f |  is  the  failure  restrained). This  stress may  be  in  pure  considered  41 subject to s i z e e f f e c t s , a c c o r d i n g  f  k C c  V  = (F  to  ^ k )  (47)  C  c  * (48)  where * F  = f a i l u r e stress  c  for a u n i t k  c  in pure compression  volume,  = s i z e e f f e c t parameter i n compression,  V = t o t a l volume of the compression that 3.7.2. S i z e Let  domain under  i s , the  e n t i r e member.  e f f e c t in tension  F,p be  Equations (44)  the and  strength in (45)  we  pure t e n s i o n .  have at any  Then,  from  probability level:  (49)  where k  fc  is  the  size  t o t a l volume and In the context  V  T  follow  the  is  of the  f i n i t e element i and in F i g u r e 12.  effect factor  The  in tension,  the domain of the analysis  the  tensile  V is  the  stresses.  presented here, c o n s i d e r  a  l o c a l ^-coordinate system, as shown  stresses  within  the  element are  s t r e s s s t r a i n r e l a t i o n s h i p as  indicated  assumed to Figure  4.  Let  us introduce a l o c a l c o o r d i n a t e  = rjh. The t e n s i l e  77 (0 ^ 77 ^ O s u c h that  s t r e s s e s w i l l be l i n e a r  in y, or a = r\a  where a„ i s the maximum s t r e s s at the edge y = h.  <)  5  ^ '  J  H  l  :  Figure  i  y  2A  1  12. S t r e s s p r o f i l e a c r o s s the s e c t i o n .  43 Equation  (49) can then be expressed as  2^BHA / d£ / drj a  =F V  u  (50)  U  T  where N i s the number of elements. Since a =  f'hU) / •H  k a  t T  f'k d£ / r? J  ({)  H  A(£)  1 fc  / t y_i  dTj = k  a  + 1  0  Oy(£)r),  k (£)d£  t T  f j i ;  H  then, Equation (50) becomes  —;  1  v^/hU) r r/ J * U)d£  = Fm  (a max) ^U/h(£)  a <«)  2(k +l)N 4-^1  a max  k  2N(k.+l)4—'I H t i=l  or,  T  T  k  r  T  finally, k  kt  fc  T  T  fc  T  H  k  fc  T  '-I  The l o c a t i o n of the n e u t r a l a x i s , h ( £ ) , where the s t r e s s e s a change s i g n , can be obtained by i n t e r p o l a t i o n of the  stress  field. The element  implementation computer program  of  the  procedure  f o l l o w s the  in  the  equations as  finite derived  above. A summary of the steps f o l l o w below. 1.  c (£) T  is  determined  at  a l l points  £  and  for a l l  elements; 2.  o b t a i n the l a r g e s t of the  0  T  ( £ ) ,  a max to normalize T  the  44  stresses. 3.  o b t a i n h(£) f o r  any c r o s s s e c t i o n .  compression w i l l r e s u l t 4.  Integrate  over  each  A section f u l l y  in  add,  to  i n h(£) = 0. element  and  according  Equation (52). 5.  Compare the according  a max T  with  the  maximum  stress  possible  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 c o n s i s t s of a number of  subroutines  which read the s t r u c t u r e ' s geometry and load data, c a r r y out numerical  i n t e g r a t i o n , decompose m a t r i c e s ,  equations  and checks the convergence of the s o l u t i o n v e c t o r .  The  program enables  beam-columns of has been p r o v i d e d to cpu  the user  A time  or  subroutine  the amount of computer time  s o l v e each s p e c i f i c problem. seconds. A l i s t i n g of  of  t o a n a l y i z e beams, columns  various configurations. to give  s o l v e s system  used  T h i s time i s c a l c u l a t e d  the program has been p r o v i d e d  in in  Appendix A.  3.9. DISCUSSION The  analysis  developed  here  offers  numerous  p o s s i b i l i t i e s . The m a t e r i a l behaviour law can be m o d i f i e d to study  different  parameters i n a  materials  or  the  s i n g l e m a t e r i a l . Also  effect  of  several  the dimensions of  a  45  member  cross-section,  l a t e r a l l y acting In the  the  eccentricity  chapter,  the model  c o n s i d e r i n g some  problems  f o r which  experimental or t h e o r e t i c a l  account t o r s i o n a l or out  anticipated  not that  experimental  modified  to  Q  along  data,  accommodate  approximation f o r the f o r small  are  by  available  into  analysis.  It  creep  is  be a s i g n i f i c a n t v a r i a t i o n the length  However, without l o s s of g e n e r a l i t y , reliable  there  verified  here does not take  i n the  there c o u l d E  w i l l be  of plane deformations. A l s o  included  modulus of e l a s t i c i t y  suitable  load,  results.  computer program developed  e f f e c t s were  axial  loads and support c o n d i t i o n s can be v a r i e d .  following  The  of  the  such  of the  can  variations  be  and i n t e r m e d i a t e l e v e l s  used  of s t r a i n ,  o b v i o u s l y can not be e x t r a p o l a t e d to very l a r g e  of  easily  in E . o  stress-strain relationship  of  member.  and i n the presence program  also  strains.  The is but  4. PROGRAM V E R I F I C A T I O N  4.1. INTRODUCTION In  t h i s chapter,  developed  the f i n i t e  element computer  i n the p r e v i o u s chapter  program  i s v e r i f i e d with r e f e r e n c e  to 1.  theoretical  results  from  the  theory  of  elastic  beam-columns, and 2.  the r e s u l t s of  an e x t e n s i v e experimental  l a r g e number of timber  The  members i n s t r u c t u r a l s i z e s  r e p o r t e d by Bleau  (1983) and Buchanan  t e s t m a t e r i a l was  SPF lumber, purchased  l e n g t h s as  'Number  program on  2 and B e t t e r '  = 32.3  Mpa  Q  i n 16ft. (4.88m)  grade i n Quebec,  = 9660 Mpa,  and t e n s i l e s t r e n g t h f  fc  compressive  = 30.35  [as  (1984)].  Canada.  The program i s v e r i f i e d using the mean t e s t r e s u l t s , modulus of e l a s t i c i t y E  a  namely  strength f  Mpa.  4.2. COMPARISON OF RESULTS The and  first  comparison c o n s i d e r s  analytical results  [3]  the computer program's e l a s t i c p r e d i c t i o n s using m = -1,  where  m  is  the  s t r e s s - s t r a i n curve p r e s e n t s p l o t s and  slope  of  the  falling  i n compression. tables  of a x i a l  The  branch second  load versus  of  comparison slenderness  r a t i o to compare the mean maximum l o a d from t e s t s with the present a n a l y s i s p r e d i c t s f o r s e v e r a l end 46  the  what  eccentricities  e. No s i z e e f f e c t s are  c o n s i d e r e d . The t h i r d comparison  s i m i l a r t o the second one, except that e f f e c t phenomenon  is  e f f e c t of v a r y i n g k  taken  into  is  i n t h i s case the s i z e  consideration,  and  the  f o r a chosen k. i s e v a l u a t e d .  4.2.1. P r e s e n t a t i o n o f R e s u l t s Table  1 shows  a comparison  theoretical results  [3]  and  u n i f o r m l y loaded f i x e d ended plotted  computer  and  non-linear  predictions  beam. The data  for a  of Table 1 i s  i n F i g u r e 13.  Qo [kn/ml  Wraax [m] iTlmoshenko]  00.000 06.885 13.771 20.656 27.541 34.426 41.312 48.197 55.082  Table  of l i n e a r  1.  0.00000 0.01291 0.02447 0.03516 0.04406 0.05162 0.05874 0.06497 0.07076  Maximum d e f l e c t i o n s  loaded beam. (E  [program]  Ilinear ]  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  of  a fixed  ended  uniformly  = 10000 Mpa, 2x4-in s e c t i o n , L = 2m)  48  0.11  L= 2m E = 10000 Mpa Q  Q-kN/m o  Q  60  40  ^ Linear, — Timoshenko , x program  2x4-in SPF Figure  13. Comparison of program  (with m=-1)  and  analytical  results [3].  Table the  2 shows a comparison computer  program  of f a i l u r e loads as obtained  to  test  results  for  e c c e n t r i c i t i e s e. A g r a p h i c a l p l o t of the data i s shown  in Figure  14, with  Similar  results  with s i z e  Figures  15(a) and 15(b).  no s i z e  effects  different  in this  effects  i n c l u d e d are  by  table  considered. shown  in  49  COMPUTER RESULTS  e = 2mm c  c  3.37 5.10 6.74 8.99 11.24 14.61 16.85 19.10 20.22 21.35 24.72 25.80 28.10 32.60 35.96 40.45  TEST  e = 39mm  e = 2mm  Pmax [Kn] 100.953 100.111 98.008 95.064 90.437 80.131 70.022 60.125 55.170 50.897 40.024 36.934 31.793 24.220 20.054 15.974  RESULTS e = 39mm PmaxtKn]  41.327 40.066 38.593 36 . 490 34.177 30.601 28.175 25.676 24.442 23.376 20.356 19.410 17.759 14.829 13.194 11.258  104.35  48.21  69.02  32.68  48.75 34.98 20.52  24.71  —  14 .11  e -it-  Table  2. Axial load-slenderness data for a pin-ended beam (size effect neglected).  Data Input : 2x4-in SPF section mean E , f , f Q  c  fc  2x4-in  Figure  14.  A x i a l load-Slenderness Table 2,(no  plots  f o r the  data  of  size e f f e c t ) .  110  0  10  20  30  40  C c  51  110  10  H  1  5  Figure  | 15  |  1 25  1  ;  | 35  1  15(a). A x i a l load-Slenderness curve with s i z e taken i n t o account, e =  Input Data: mean E^, f and f. o c t k = 20.0 and k. = c t  5.0  2mm.  1 45  effect  Figure  15(b). A x i a l l o a d - s l e n d e r n e s s curves with f o r a constant  Tests, 2-  k  4-  k  c  = 10, = 100,  c  l - k 35  k .  = 5 k  = 20  c  - k  = c  150  fc  v a r y i ng  53 4.3.  DISCUSSION The  there the  is a  as  presented  relatively  in  good one  for  compression  this  good agreement  p r e d i c t i o n s by t h e computer  a very For  results  in  the  results.  For very  slightly  below. S e v e r a l e x p l a n a t i o n s  the  discrepancies.  this  finite and  feature  general  Application will  of  be d i s c c u s s e d  taken  slightly effect  as  be  and  developed study  test  study  one may  i n the here  the  columns  is  to that  n o t be  a  Nevertheless, analysis,  remains a  behaviour  program  are  the  powerful  and  design  and  beam-columns.  t o wood  beam-columns  i n the f o l l o w i n g chapter.  input  15(a), into  have these  when k  the  respect  considered. Also,  members, w h i l e  the  the  be p u t f o r w a r d  obvious  this  changed  computer  improved w i t h  was  for  timber  in Figure  in compression  short  to  of t h i s  As shown  ratio.  range,  than  the a c t u a l behaviour.  may  tool  can  The most  used  element t e c h n i q u e  considerations  are  curve  r e p r e s e n t a t i o n of  since  being  slenderness  higher  and  s h o r t members t h e p r o g r a m p r e d i c t i o n s  stress-strain  true  slightly  that  the t e s t  intermediate  predictions  these  show  t h e agreement  a high  program  explain  are  between  program,  members w i t h  members  chapter,  little size  =  program, to  it  the  t h e ones  i s noted  significance  effects  20.0  reason  k.  =  5.0  results  are  where no  size  that size  effects  for  slender  p l a y a major  i n t e r m e d i a t e members. The  and  very role  in  very  for this  is  that  54 for short members, the volume s u b j e c t e d t o t e n s i o n i s or  non-existent.  controlled  by  For  the  i n s t a b i l i t y . When k  c  modulus  k  members, of  the  elasticity  i s very l a r g e , the r e s u l t s  the same as the ones size effect  slender  i n Table  2, meaning  small  failure and  member  obtained are  that there i s no  f o r l a r g e k « I t appears from F i g u r e 15(b)  = 20.0 g i v e s a best  c  f i t to the t e s t  is  results.  that  5. APPLICATIONS  5.1.  INTRODUCTION  The  application  of  the  beam-columns w i l l be d i s c u s s e d can  handle m u l t i p l e  program  to  solve  wood  i n t h i s c h a p t e r . T h i s program  spans with d i f f e r e n t  l o a d and  support  c o n f i g u r a t i o n s . Among them are the ones shown i n f i g u r e 16.  Q  Q  p  J.  P  (b)  (a)  Q o P  -U  Q p  I I I I I I I I l i  \  (c)  i i i i i  (d)  Figure16. Some loading cases and support  The  members  are  assumed  responses can u s u a l l y be deflection  curve  r r  or  any  to  be  conditions.  prismatic.  The  r e p r e s e n t e d as load versus other 55  convenient  way  desired centre for  a  56 particular  lateral  load Qo.  obtained, the maximum the peak of  Once the  loads can be  the c u r v e s . The  complete  curves  e a s i l y determined  from  computer program developed  t h i s study p r o v i d e s an e a s i e r approach  t o the above  a p p r o p r i a t e i n f o r m a t i o n . In  F i g u r e 16, the l a t e r a l about  the major  assumed that  loads  a x i s of  weak a x i s  a l l cases  Q or Qo cause bending  the c r o s s - s e c t i o n . buckling  -  b u c k l i n g a r e e f f e c t i v e l y prevented so that f a i l u r e caused by lateral  e x c e s s i v e bending  load  assumed that  i n the  plane of  lateral  load Qo  i s applied  of  further  torsional i s always  the  . In performing the numerical procedure, the  the  moments  It i s  and l a t e r a l  in  process  in that one gets the maximum load d i r e c t l y by s u p p l y i n g program with the  are  applied i t is  first  and  maintained at a constant value as the a x i a l compressive  load  P i n c r e a s e s or decreases.  5.2.  NUMERICAL EXAMPLE As a numerical example,  considered i n t h i s  study, using  mean v a l u e s f o r E , f Q  k  c  case  = 10.0 has been used  and f . f c  (c) i n F i g u r e 16 has a 2x4-in  Also k  fc  been  SPF s e c t i o n  and  = 5.0, m = 0.02  and  i n o b t a i n i n g the P versus Qo r e s u l t s  as shown i n F i g u r e s 17(a) and 17(b).  Figure  17(a).  Ultimate  strength  interaction  simply supported columns subjected to uniformly  curves  for  distributed  load.  0  Figure  02  17(b).  0.4-  04  Non-dimensionalized  0.8  ultimate  i n t e r a c t i o n curves of Figure 17(a).  1  Qo/Qou  strength  58  5 . 3 . OBSERVATIONS From F i g u r e s P-Qo  17(a) and  17(b), i t can  be n o t i c e d  that  r e l a t i o n s h i p s p r e d i c t e d by the computer model are not a  l i n e a r one as i t i s practice  for  research  is  normally assumed i n the c u r r e n t  different needed  here  s i m p l i f i e d design procedure  slenderness in  order  f o r wood  ratios. to  come  design  Additional up  beam-columns.  with  a  6. R E L I A B I L I T Y ANALYSIS  6.1.  INTRODUCTION T h i s chapter d e s c r i b e s  reliability  analysis  problem to be represents and  of  studied  the procedure f o r the a  wood  i s as shown  B represents  the  member.  in Figure  the a p p l i e d a x i a l compressive load  l i v e l o a d s ) . L represents  H and  compression  structural The  18; where  P  ( f o r only dead  the length of the member while  height  and  breadth of  the  cross  section.  yfS~r  Figure  The that It  ^  18. T y p i c a l problem f o r r e l i a b i l i t y  reliability  i twill  of a member simply  perform as intended  i s influenced  c a p a c i t y of the general,  IL Q , B ,u H SPF  by the  evaluation.  means the p r o b a b i l i t y  in a prescribed  demands on the  s t r u c t u r e t o respond  s t r u c t u r e and the  to those demands.  one can d e f i n e a performance or f a i l u r e  to c h a r a c t e r i z e the  s t a t e of the  some performance c r i t e r i o n . T h i s  59  situation.  function  structure in relation  In G to  f u n c t i o n G can be expressed  60  as G = C - D  where C = structural capacity D = demands on the s t r u c t u r e .  The f u n c t i o n G as d e f i n e d above i s p o s i t i v e whenever the c a p a c i t y exceeds the demand,  therefore  the performance c r i t e r i o n . On the other w i l l be negative resulting  in  whenever  the  structure  structure  failing  is  on  the  not  meets  hand, the f u n c t i o n G  the demands exceed the  performance. When the f u n c t i o n G the  the s t r u c t u r e  meeting  capacity,  the  required  i s e x a c t l y equal to  threshold  between  zero,  meeting  to meet the performance c r i t e r i o n , and such a  i s d e f i n e d as " l i m i t  and state  state". e  The p r o b a b i l i t y of f a i l u r e reliability  failure  of  the  . Thus Pj = 1.0 -  According  p^ i s the compliment  reliability  to the d e f i n i t i o n of  i s then given  as  G above, the p r o b a b i l i t y  of  61 p  Each design problem variables,  it  w i l l contain a  may be random,  Thus, i f some of  i s obvious that G w i l l  probability distribution knowledge of the  set of  intervening  and depending on the nature of the problem,  of the v a r i a b l e s function.  ( G < 0)  = Probability  f  basic  the basic  variables  distribution are  random,  be i t s e l f a random v a r i a b l e . for G  the i n d i v i d u a l  variables,  obeying some  c o u l d be  derived  r e s u l t would be  The  from  probability distributions  and the  some  as shown  a  for in  F i g u r e 19.  G  Figure  =  0  19. P r o b a b i l i t y  G  density  G  function  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 curve to the l e f t of the  o r i g i n G = 0.  f o r the v a r i a b l e G .  be the area under the If t h i s  probability  62  of f a i l u r e design  exceeds some d e s i r e d value,  v a r i a b l e s would be changed and  i t meets the r e q u i r e d for  t a r g e t . The  G c o u l d be obtained  i n t e g r a t i o n s and  one  or more of  p^ r e c a l c u l a t e d  probability  the  until  distribution  by a n a l y t i c a l means using  multiple  the j o i n t 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 b a s i c v a r i a b l e s .  This i s a  very  tedious and  difficult  approach. MonteCarlo p r o b a b i l i t y of  simulation failure  approach the value of  in  can  be  used  to  an approximate  G i s computed f o r  values  for  the  basic  p r o b a b i l i t y d i s t r i b u t i o n s , and are  i n v o l v e d , the procedure  manner. In  this  i s estimated  negative.  variables  the  a l a r g e number  combinations of the b a s i c v a r i a b l e s and p^ the p r o p o r t i o n of times the G was  obtain  must  The  obey  from  s e l e c t i o n of their  when more than two  becomes d i f f i c u l t ,  of  joint  variables  tedious  and  expensive. In the f o l l o w i n g s e c t i o n , an approximate and  fast  procedure f o r e s t i m a t i n g p^ w i l l  6.2.  THE  0 METHOD FOR  In order  integrations  or  [1974] introduced  discussed.  R E L I A B I L I T Y ANALYSIS  to estimate  s u f f i c i e n t accuracy  be  the p r o b a b i l i t y of f a i l u r e p^  but  computer  without r e s o r t i n g simulations,  the concept of  uncorrelated  random v a r i a b l e s X-,  Hasofer  reliability  geometric approach. Thus, f o r a design i  to  with  complicated and  index 0  Lind using  problem c o n t a i n i n g  = 1,..,N, with mean  N X.  63  and standard d e v i a t i o n x^  is  introduced.  a^, a set of "normalised "  These  variables  standard d e v i a t i o n equal to 1 . 0 ,  x. =  X.  -  have  variables  zero  and are given  mean  and  as  X.  3,  ( ) 5 3  °i The f a i l u r e new,  f u n c t i o n can  now  normalised v a r i a b l e s  f i g u r e below, i n which  be expressed  x^ as shown  i n terms of  the  schematicaly in  the  the h o r i z o n t a l plane represents  space of the v a r i a b l e s x and the v e r t i c a l a x i s the  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 L i n d showed that the r e l i a b i l i t y be i n t e p r e t e d as the minimum  distance  between  index 0 can  the o r i g i n  0  64  and  the l i m i t  can  be s o l v e d by  Hasofer  and  state  G . T h i s i s a geometric problem  successive  Lind's  i t e r a t i o n s u s i n g , f o r example,  proposed  algorithm.  p r o b a b i l i t y of f a i l u r e i s obtained  p  normally  /3, the  Knowing  from  (54)  = #(-0)  f  where $ i s the s t a n d a r d i s e d p£ t o be exact,  which  Q  we  normal p r o b a b i l i t y f u n c t i o n . For  r e q u i r e that a l l the b a s i c v a r i a b l e s  d i s t r i b u t e d and G be l i n e a r  be  i n the b a s i c v a r i a b l e s .  F i g u r e 20 shows the case when the mean p o i n t belongs to the "safe domain" G > 0. The combinations of x^ which correspond to G = 0 (the l i m i t  s t a t e ) are represented  6.2.1. R a c k w i t z - F i e s s l e r This i s in actual and  L i n d Algorithm  by the curve G . Q  Algorithm  f a c t the m o d i f i c a t i o n of the  i n order  Hasofer  t o improve the e s t i m a t i o n  of the  p r o b a b i l i t y of f a i l u r e . The m o d i f i c a t i o n r e f e r s t o the when  the  Fiessler  basic  variables  non-normal.  [1978], suggested a t r a n s f o r m a t i o n  random v a r i a b l e s X^ standard  are  of the  i n t o a set of normalised  v a r i a b l e s z^ using the f o l l o w i n g  z. = # ~ [ F ( X . ) ] 1  Rackwitz  case and  original  uncorrelated  transformation  (55)  65 where $  is  the standard  f u n c t i o n and  F(X^)  is  for the v a r i a b l e X^.  normal  The  standard  improves the p r e d i c t i o n of  algorithm  v a r i a b l e z^. T h i s  p^ because i t  d i s t r i b u t e d . This  accepted norm f o r the e v a l u a t i o n  6.3.  PROBLEM  To  algorithm  is  of  the  variables  presently  of the r e l i a b i l i t y  index  the 0.  FORMULATION  above to normal p r a c t i c e , l e t us c o n s i d e r 18. The  cross  are B f o r width and represented  and  modification  meets one  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  of F i g u r e  function  from Hasofer  c o n d i t i o n s mentioned e a r l i e r namely, that a l l  be normally  as L.  theory  s e c t i o n a l dimensions of the  H f o r depth. The I t i s assumed  D +  where D = demand P  D  = dead l o a d  P  L  = live  load  P  L  column  l e n g t h of the column i s  simply  supported under live  demand on the s t r u c t u r e i s the a p p l i e d l o a d P.  D = P = P  discussed  the column problem  a x i a l compressive l o a d P ( f o r both dead and  If  distribution  the cumulative d i s t r i b u t i o n  L i n d i s then used f o r the new  two  probability  loads).  Thus,  an The  66  d = -°P DN  (56)  /  (57)  = P  where  LN  P  LN  =  n  o  m  i  n  a  l  (design)  live  load  P  DN  =  n  o  m  i  n  a  l  (design)  dead l o a d  Then D = P  where d and  L N  [7id  (58)  + /]  / are considered  f a c t o r 7, i s a constant  to be  random v a r i a b l e s . The P d e f i n e d as 7, = -SH or the r a t i o of LN P  nominal dead load to  nominal l i v e l o a d .  The c a p a c i t y C i s  the maximum l o a d , Pmax, the member can c a r r y ; thus C = Pmax = P{E ,fc,ft,B,H,L,m] 0  and  the f a i l u r e  (59)  f u n c t i o n can be expressed as  G = C - D (60)  G = Pmax - P  where f f  fc  = strength  i n compression  = strength  in tension.  m = slope of s t r e s s - s t r a i n curve i n compression. The  problem can now be  s t u d i e d using the  Rackwitz-Fiessler  a l g o r i t h m and the f i n i t e element computer program  developed  67 in part 3 of t h i s t h e s i s . However, i t i s convenient f o r the purpose of f u t u r e code development  t o modify Equation  (60)  above to b r i n g i n the design equation format adopted f o r the code.  6.3.1.  For  Code Design Equation members s u b j e c t e d  to pure  a x i a l compression,  the  Canadian Code, CAN3-086.1-M84 (1984) s p e c i f i e s the f o l l o w i n g desing e q u a t i o n .  D DN  A  P  +  A  L LN P  * *p  A  F  c  K  C  (  6  1  )  where #p = performance  factor  i n compression.  A = c r o s s s e c t i o n a l area of member. K F  = slenderness factor  c  = F i f t h p e r c e n t i l e compression  c  ( a , a ) = load f a c t o r s D  L  strength  (1.25 and  1.5 r e s p e c t i v e l y ) .  6.4.  THE G FUNCTION FOR THE PROBLEM Considering  the  limiting  o b t a i n the f o l l o w i n g e q u a t i o n :  case of  Equation  (61), we  68  P  LN  [ a  D?i  +  a  L  *p  ] =  A  F  c  c  K  (  6  2  )  where P L  L N N  <p AFcKc = -B V L +  (63)  a  Combining Equations (60), (62) and (63),we can the  express  f a i l u r e f u n c t i o n as 0 AFcKc G = Pmax - -E [ ,d + / ] a 7,+a 7  D  L  or  G = P{E , f  0 A Fc Kc f ,B,H,L,m} - -E [ a 7,+a D  For  the purpose of  7 l  t h i s study, the f o l l o w i n g  have been c o n s i d e r e d  d + /]  (64)  L  variables  random  modulus of e l a s t i c i t y E  Q  compressive s t r e n g t h f t e n s i l e strength f  f c  dead l o a d v a r i a b l e d l i v e load v a r i a b l e /  and  the f o l l o w i n g have been  considered  with mean average v a l u e s height  of c r o s s  section H  to be  constants  69 breadth  of c r o s s s e c t i o n B  l e n g t h 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  above i s a t t a c h e d  i n Appendix  program requests  the number  case 5 ) ,  the type  of t h e i r  d i s t r i b u t i o n code), to c h a r a c t e r i z e the the  and  the d e r i v a t i o n  A. As part of the input, of random  the  variables (in this  distribution  (according  the r e l e v a n t parameter  d i s t r i b u t i o n s . The  to a  information  program can  accept  following distributions Code  Distribution  1  Normal  2  Lognormal  3  Weibull  A  Gumbel  5  Ranked Data  The f i x e d parameters 7 , ,  <j>^, a  D  and a  the user  f o r each p a r t i c u l a r problem.  computes  the  function G  f i n i t e element value of G problem gradient  discussed vector  The subroutine  The GXPR  routine  gradient v e c t o r DELTA. in  are provided  and i t s g r a d i e n t by  subprogram.  and the  L  this  thesis  corresponding  to  by GXPR  calling  the  returns  the  For the  column  the  elements  of the  the  first  3  random  70  v a r i a b l e s were obtained n u m e r i c a l l y , while the remaining two were  obtained  by  differentiating  the  failure  function  e x p l i c i t l y . Thus, the t o t a l elements of the g r a d i e n t  vector  c o n s i d e r i n g o n l y f i v e random v a r i a b l e s a r e o b t a i n e d as  Delta(l) =  G(E  *) - G(E ") — 2AE^o  2  2  G(f ) - G(f ~) = S $— +  Delta(2)  2Af  G(f. ) S  DeltaO) =  c  - G(f/) £ —  +  2Af Delta(4)  t  <t> AFcKc = — E a 7i + a D  7  l  L  <(> AFcKc Delta(5)  = a 7i a +  D  The  fixed  L  parameters are passed  COLUMN through  a COMMON b l o c k .  onto the r o u t i n e s GXPR  and  71 6.6.  RELIABILITY  RESULTS  Keeping the r a t i o  7 l  = 1.0, m = 0.02 and using a  2x4-in  SPF s e c t i o n , the f a c t o r <p^ was changed and the corresponding 0  reliability  index  slenderness  ratios.  reliability  index  was computed F i g u r e 21  f o r columns of  shows  the r e s u l t s  0 as a f u n c t i o n of the performance  0p f o r the case of no s i z e e f f e c t c o n s i d e r e d Figure  22  shows  reliability  i n c l u d e d i n the computer for  the  two cases  results  Q  s t u d i e d , the  factor  size  effects results  f o l l o w i n g information  : 3-parameter W e i b u l l .  Scale = 6738.0 Mpa L o c a t i o n = 3514.0 Mpa Shape = 3.97 Mean = 9660.0 Mpa  (2) f  with  f o r the  i n the program.  program. In o b t a i n i n g the  been used f o r the random v a r i a b l e s .  (1) E  different  : 3-parameter Weibull  Scale = 33.845 Mpa L o c a t i o n = 0.0 Shape = 7.8559 Mpa F i f t h p e r c e n t i l e = 15.87 Mpa  (3) f. : 3-parameter Weibull  has  Scale = 29.861 Mpa Location  = 4.03 Mpa  Shape = 2.911 Mpa  (4) d = dead load v a r i a b l e : Normal Mean = 1.0 Standard d e v i a t i o n = 0.15  (5) / = l i v e load v a r i a b l e : normal Mean = 0.75 Standard d e v i a t i o n = 0.15  Figure  22.  Reliability the  r e s u l t s with reliability  size  effect  program.  included  in  75 6.6.1. D i s c u s s i o n of r e s u l t s From the  r e s u l t s of F i g u r e  performance f a c t o r order  j3  =  4.0  = 0.75,  Figure  reliability  index  Considering  the two  22  index  0 of  f o r a l l slenderness  shows  for  noted that f o r a  a reliability  i s achieved  considered.  21 i t i s  the  a  small  same  increase  performance  cases, a performance  ratios of  the  factor  factor 0  the  =  <p . cr 0.75  Cr  can be taken  as a reasonable  value t o be  i n c l u d e d i n the  c u r r e n t design p r a c t i c e s f o r columns of any l e n g t h . The  procedure  reliability  outlined  column. A  in  the  servive  study  load e f f e c t Civil  f o r beams  taking  i s currently  Engineering  of  U n i v e r s i t y of B r i t i s h Columbia. I t  w i l l be of i n t e r e s t  f u r t h e r research, t o i n t e g r a t e the  model developed here  this  study.  the  over the l e n g t h of the  of  not take  for  account  reliability  Department  chapter into  i n t o account the d u r a t i o n of done  this  a n a l y s i s of columns does  the d u r a t i o n of load e f f e c t l i f e of the  in  been the for to  7. CONCLUSIONS AND  7.1.  RECOMMENDATIONS  CONCLUSIONS From the r e s u l t s of part  that  the  finite  deformations and  one  element  non-linear  of t h i s study,  analysis,  and  reliable  input  e l a s t i c i t y , compressive and including  size  parameters k for  the  effects  = 20.0  c  maximum  load  r e s u l t s . However, k short and some  and  c  that  = 5.0,  does have  intermediate  large  columns,  for  model  model does r e q u i r e  information  on  modulus  the  fairly  of  results  size  the computer  Pmax agree  effect  predictions  well  with  test  significance influence and  should  be known  for with  accuracy. For the r e l i a b i l i t y  it  fc  including  t e n s i l e s t r e n g t h s . The  show k  seen  m a t e r i a l p r o p e r t i e s , can  wood column behaviour s a t i s f a c t o r i l y . The accurate  i t is  i s observed that the  results  i n part two  of t h i s  study,  c u r r e n t performance f a c t o r s <l> ,  as  hr  given  i n CAN3-086.1-M84 (1984),  are more c o n s e r v a t i v e  what t h i s model p r e d i c t s . A value of 0p = 0.75 a reasonable that i f  one  t h i s new  f o r a l l slenderness value of  <f>^ i s  appears to be  r a t i o s . I t i s estimated adopted i n  the  design p r a c t i c e , i t w i l l  give r i s e to a r e l i a b i l i t y  of the order of 4.0.  a lower 0 i s r e q u i r e d , a  #p should be  If  introduced  f o r short and  T h i s p o i n t s to a d e f i c i e n c y  i n the 76  than  current index  /3  different  intermediate  "column formula"  columns. giving  77 the slenderness should  reflect  adjustment  factor K .  Idealy, t h i s  c  the changes due t o slenderness  factor  i n such a  that the same <j>^ - 0 r e l a t i o n s h i p be obtained  way  f o r a l l column  lengths.  7.2.  RECOMMENDATIONS It  i s recommended that a performance f a c t o r <f>^ = 0.75 be  used i n the c u r r e n t r a t i o s . However, there  i s need  p a r t i c u l a r , the  design  prior to  p r a c t i c e f o r a l l slenderness  to  do  adopting  more  research  this  research  should  in this  cover  e f f e c t s , and the case of c o r r e l a t e d  recommendation, area.  duration  of  In load  variables; neither  of  which has been i n c l u d e d i n the a n a l y s i s . The a p p l i c a t i o n of the R a c k w i t z - F i e s s l e r a l g o r i t h m involved to cases  some  be  r e q u i r e s a l l the v a r i a b l e s  u n c o r r e l a t e d . However,  or more  of the  in  of  the  will  depend on strength f  combined a x i a l  the modulus c  of  and the t e n s i l e  v a r i a b l e s are p a r t i a l l y with, u s i n g  f o r example  literature  [12],  and l a t e r a l  elasticity E , o  the  loads  f c  must be  the procedures a v a i l a b l e  before  using  the  or will  compression  s t r e n g t h f . For lumber,  c o r r e l a t e d and t h i s  be  problems  i n t h i s r e p o r t , the s t r e n g t h of beams, columns  beam-columns under  algorithm.  practical  intervening variables  c o r r e l a t e d . For example, i n the context discussed  some  these dealt i n the  Rackwitz-Fiesler  78 There i s not effects in  both  enough data a v a i l a b l e tension  compression,  size  hence  further  p r a c t i c a l as w e l l as t h e o r e t i c a l study i s necessary  i n order  to come up w i t h a r e a l i s t i c  and  at present on  design recommendation a p p l i c a b l e  to lumber of a l l grades and s p e c i e s .  REFERENCES 1. Chen,W.F., and Astuta,T., Vol.  II,  Space  Company, New  2.  Behaviour  and  Design,  Columns,  McGraw-Hill  Book  York, 732p.  Zienkiewicz,0.,  Engineering  1976, Theory of Beam  1971,  Science.  The  Finite  McGraw-Hill  Book  Element  method  Company,  in  London,  England.  3. Timoshenko S., and Woinowsky-Krieger  S.,  1959, Theory  p l a t e s and s h e l l s . McGraw-Hill Book Company, New  of  York.,  pp.  396-428.  4. Malhotra, S.K., of  and Mazur, S.J., 1970, B u c k l i n g  s o l i d timber Columns, Engrg.J.  Strength  13(A-4), I-VII.  5. Larsen,H.J., and T h i e l g a a r d , E. (1979)., L a t e r a l l y timber columns  , J.  6. Bleau, R.,  (1983).,Etude de R e s i s t a n c e  bois  de  Qualite  Struct.  Commerciale,  U n i v e r s i t y of Shebrooke the  Div.,  ASCE,  "Thesis  loaded  105(7) 1347-1363.  des Colonnes  presented  en  to  the  at Quebec, i n p a r t i a l f u l f i l m e n t  of  requirements f o r the degree of Doctor of P h i l o s o p h y " .  7. Buchanan, A.  H.,  S t r e n g t h Model 79  and design methods  for  80 bending and  axial  load  interaction  in  timber  members.,  "Thesis p r e s e n t e d to the U n i v e r s i t y of B r i t i s h Columbia, Vancouver,  B.C.,  in p a r t i a l  f u l f i l m e n t of the  at  requirements  f o r the degree of Doctor of P h i l o s o p h y " .  8. Zann, J . J . , bending and  (1985)., S t r e n g t h of  compression.,  "USDA  paper FPL 391., Washington,  9.  Hasofer,  A.M.,  and  lumber under  Forest  Service  Lind,  N.C.,  1974.,  Mechanics D i v i s i o n , ASCE, V o l . 100, No. EM1,  Rackwitz,  R.,  R e l i a b i l i t y Under Structures  "Exact  Fiessler,  Combined Load  B.  Engineering 111-121.  "Structural Computers  and  index  1979., "A d i s c u s s i o n on the a p p l i c a t i o n of concept  12. Der K i u r e g h i a n , A.,  Mechanics  1978.  Sequences",  to  wood  structures",  J o u r n a l of C i v i l E n g i n e e r i n g , V o l . 6, No.  Incomplete  pp.  and  , V o l . 9, pp. 489-494.  11. F o s c h i , R.O., the s a f e t y  and  Research  D.C".  I n v a r i a n t Second Moment Coe Format", J o u r n a l of  10.  combined  Canadian  1, pp. 51-58.  1986., " S t r u c t u r a l R e l i a b i l i t y Under  P r o b a b i l i t y Information", J o u r n a l of , ASCE, V o l . 112, No.  1, pp. 85-104.  Engineering  81  A P P E N D I X  A  1. Program COLUMN.FOR (source code)  82  2.  96  Sample Input/Output F i l e s  3. Program RBETA.FOR  (source code)  4. Sample Input/Output F i l e  99  119  PROGRAM COLUMN.FOR  83  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 V e r s i o n 2.0 * * * 2 A u g u s t , 1987 * * * * A PROGRAM FOR THE CALCULATION OF THE ULTIMATE LOAD ON 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. 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. * * * ************************************************************ 1  IMPLICIT REAL*8(A-H,0-Y) CHARACTER*20 NAMED 1 ,NAMEA1 ,NANS DIMENSION I X ( 2 1 , 4 ) , F ( 8 ) , N B C ( 2 1 ) , T K O ( 6 7 2 ) , X E ( 8 ) 1, R ( 8 4 ) , X O ( 8 4 ) , X ( 8 4 ) , B ( 8 4 ) , B 1 ( 8 , 8 ) , B 2 ( 8 , 8 ) , B 3 ( 8 , 8 ) , B 4 ( 8 , 8 ) 2, B 5 ( 8 , 8 ) , B 6 ( 8 , 8 ) , B 7 ( 8 , 8 ) , B 8 ( 8 , 8 ) , B 9 ( 8 , 8 ) , Y ( 5 ) , R E ( 8 ) , X P ( 8 4 ) 3, Q ( 2 0 ) , I Q ( 2 0 ) , E S T R ( 7 ) , F I ( 7 ) COMMON/C1/GAP(5),GAW{5),EN 1 ( 8 , 5 ) , E M 1 ( 8 , 5 ) , E M 2 ( 8 , 5 ) , NGAUSS COMMON/C2/DIFP, NINT COMMON/C3/DEFL,PDEFL ********************************************************** * DEFINE VARIABLES *  ********************************************************** NELEM NJBC NBC(I) IX 1 2 3 4 EN 1 EMI,EM2 GAP GAW NGAUSS NITER TOP EPSLON FC  = = = = =  = = = = = = = = =  NO OF ELEMENTS NO OF JOINTS WITH B.C. NO OF B.C. AT NODE I B.C. CODE U UX W WX INTERPOLATION FUNCTIONS FOR u INTERPOLATION FUNCTIONS FOR w CORDINATE AT GAUSS POINT CORRESPONDING WEIGHT NO OF GAUSS POINTS MAX. NO OF ITERATIONS TOLERANCE FOR LOAD STEP TOLERANCE FOR SOLUTION VECTOR MATERIAL STRENGTH IN COMPRESSION  * * * * * * * * * * * * * * * * *  c c c c c c c c c c c c c c c c c c c c c c 6 7 8  43 44 65  75  * FT MATERIAL STRENGTH IN TENSION * * EO * MOE OF THE MATRIAL (RANDOM) = * EN SLOPE OF THE STRESS-STRAIN CURVE IN COMPR. * = * SPAN * MEMBER LENGTH = * W * WIDTH OF SECTION = * H * DEPTH OF SECTION = * E * ECCENTRICITY OF AXIAL LOAD * NEQ * = NO OF EQUATIONS TO BE SOLVED * NJOINT = * NO OF NODES = * NDOF * NO OF VARIABLES PER NODE * NODEL * = NO OF NODES PER ELEMENT * NDIMB * = NO OF VARIABLES PER NODE * LBW,LHB = * HALF BANDWIDTH INCLUD. THE DIAG. * NA * NO OF UNKNOWNS FOR TOTAL PROBLEM = = * RE * ELEMENTAL LOAD VECTOR * R * STRUCTURE LOAD VECTOR = * B * GLOBAL LOAD VECTOR RETURNED = * TKO * GLOBAL TANGENT MATRIX = * XE * ELEMENT DISPLACEMENT VECTOR = * X * GLOBAL SOLUTION VECTOR * B1,..B9 = * ARRAYS FOR TEMPORARY STORAGE ********************************************** ************  84  WRITE( * , 6 ) FORMAT(' ENTER DATA F I L E NAME '/) READ(* ,8) NAMED1 WRITE( *,7) FORMAT(' ENTER OUTPUT F I L E NAME '/) READ(* ,8) NAMEA1 FORMAT(A) OPEN (1, F I L E = NAMED1,STATUS='OLD') OPEN ( 2 , F I L E = 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 READ (1,*) I Q ( I ) , Q ( l ) E = 0.0D0 I F (NP.NE.O) READ(1,*) E DO 65 I = 1, NJOINT 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) ) CONTINUE NEQ = NDOF*NJOINT NODEL = 2  C C * READS MATERIAL STRENGTH IN COMPRESSION (FC) AND TENSION ( F T ) , C BOTH CORRESPONDING TO THE SPECIFIED CROSS-SECTION AND THE C REFERENCE SPAN SREF. XKC AND XKT ARE THE WEIBULL S I Z E EFFECT C SHAPE PARAMETERS IN COMPRESSION AND TENSION RESPECTIVELY. C READ(1,*) FC,FT,SREF,XKC, XKT  C C  C  C  79 790 791 799  761  760 792 C C 3773  NDIMB = NODEL*NDOF LBW = NDIMB LHB = LBW NA = LBW*NEQ READ MOE AND SLOPE m OF CURVE READ(1,*) EO,EN READ PROBLEM S I Z E 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 * ADJUST STRENGTHS TO THE ACTUAL VOLUME FC = FC *(SREF/SPAN)**(1.0/XKC) FT = FT *(SREF/SPAN)**(1.0/XKT) OBTAIN SHAPE FUNCTIONS AND DERIVATIVES : N1,M1,M2 CALL SHAPE(DEL) WRITE(*,79) FORMAT(' TOLERANCE FOR PMAX ? '/) READ(*,*) TOP WRITE(*,790) FORMAT(' TOLERANCE FOR CONVERGENCE? '/) READ(*,*) EPSLON WRITE(*,791) FORMAT(' MAX. NUMBER OF ITERATIONS? '/) READ(*,*) NITER WRITE(*,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 FQ1 = 0.0D0 FQ2 = 1.ODO FQ3 = FQ2 NFLAG = 0 DO 792 J = 1, NEQ X P ( J ) = 0.0D0 START CALCULATIONS FOR TRIAL LOAD LEVELS CALL TIME(ZIM) ' ZIM0 = ZIM 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  4000 4001 C 80 C  81  82 83 87 180 185  C C 777 84 85 C  86 88 C  IF (NINT.EQ.1.AND.NP.EQ.O) WRITE(*,4001) FQ FORMAT(//' SOLUTION FOR P =',E15.6,* :'/) FORMAT(//' SOLUTION FOR LATERAL LOAD FACTOR=',E15.6,*:'/) I N I T I A L I S E ARRAYS DO 80 J = 1, NEQ XO(J) = X P ( J ) R ( J ) = 0.DO EXTERNAL LOAD VECTOR R IF (Q0.EQ.0.0D0) GO TO 87 DO 81 J = 1 , 8 RE(J)=0.D0 CONTINUE RE(3) = FQ*Q0*DEL RE(4) = FQ*Q0*DEL**2/3.D0 RE(7) = R E ( 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) + R E ( J J ) CONTINUE CONTINUE I F (NQ.EQ.0) GO TO 185 DO 180 J = 1, NQ JS = ( I Q ( J ) - I ) * N D O F + 3 R(JS) = R(JS) + Q(J)*FQ 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 BEGIN ITERATIONS AT THE TRIAL LOAD LEVEL CONTINUE DO 84 1 = 1 , NA TKO(I) = O.0D0 DO 85 K = 1, NEQ B(K) = -R(K) DO 645 I E = 1, NELEM I N I T I A L I Z E ARRAYS DO 88 I = 1 , 8 F ( I ) = 0.0D0 DO 86 J = 1 , I B 1 ( I , J ) = 0.0D0 B 2 ( I , J ) = 0.0D0 B 3 ( I , J ) = 0.0D0 B 4 ( I , J ) = 0.0D0 B 5 ( I , J ) = 0.0D0 B 6 ( I , J ) = 0.0D0 B 7 ( I , J ) = 0.0D0 B 8 ( I , J ) = 0.ODO B 9 ( I , J ) = 0.ODO CONTINUE CONTINUE PICK ELEMENT SOLUTION FROM GLOBAL VECTOR DO 90 J J = 1 , 8  87  .90  91 C  93 C  96  98 99 100 101 C C  K = ( I E - 1)*NDOF + J J X E ( J J ) = XO(K) CONTINUE DO 101 K = 1, NGAUSS Y(K) = 0.D0 DO 91 1=1, 8 Y(K) = Y ( K ) + X E ( I ) * E M 1 ( I , K ) CONTINUE 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) B 3 ( 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) CONTINUE DO 100 L = 1, NGAUSS 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) 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 B 5 ( 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) B 6 ( I , J ) = B6(I,J)+DEL*0.5D0*AR*(EN 1 ( I , K ) - G A P ( L ) * 1 H*0.5D0*EM2(I,K))*(E0+EN*E0)*FAC*Y(K)*EM1(J,K)* 2 GAW(K)*GAW(L) B 7 ( 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) B 8 ( 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 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) CONTINUE CONTINUE CONTINUE OBTAIN ELEMENT TANGENT MATRIX 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  88  102 105 645 C  1080  108  110 111 C C C  112  EKT = B I ( I , J ) + B 2 ( I , J ) + B 3 ( l , J ) + B 4 ( 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 CONTINUE CONTINUE CONTINUE INTRODUCE BOUNDARY CONDITIONS DO 1 1 1 U O = 1 , NJOINT IF (NBC(IJO).EQ.O) GO TO 111 DO 110 J = 1, NBC(IJO) I I = ( U O -1)*ND0F + I X ( I J O , J ) LBW1 = LBW - 1 DO 108 K = 1, LBW1 J J = I I - LBW + K IF (JJ.LE.O) GO TO 1080 I J = (JJ-1)*(LBW-1) + I I TKO(IJ) = 0.0D0 J J = II + K IF (JJ.GT.NEQ) GO TO 108 I J = (II-1)*(LBW-1) + J J TKO(IJ) = 0.0D0 CONTINUE I J = ( I I - 1)*(LBW-1) + I I TKO(IJ) = 1.0D0 B ( I I ) = 0.0D0 CONTINUE CONTINUE SOLUTION OF THE SYSTEM CALL DECOMP(NEQ,LBW,TKO,IERROR) IF(IERROR .EQ. 1) GO TO 377 4 CALL SOLVN(NEQ,LBW,TKO,B) DO 1 1 2 1 = 1, NEQ X(I) = XO(I)-B(I) 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@ XO(I) = X ( I ) GO TO 777 I ERROR = 1 GO TO 3774 WRITE(2,900) NITER, P FORMATC NO CONVERGENCE IN',13,'ITERATIONS AT P=',E13.6/) GO TO 901 (  115 430 431 900 C C C C 118  AFTER CONVERGENCE, OBTAIN STRESSES AND STRAINS AT THE CURRENT LOAD LEVEL CONTINUE EMAXP = O.0DO EMAXN = 0.0D0 SUME = 0.0D0  89  500  501  505  506  507 508  509 C C C 510  511 512  515 516  518 520  DO 550 I E = 1, NELEM DO 500 J = 1, 8 K = (IE-1)*NDOF +J X E ( J ) = X(K) CONTINUE DO 540 K = 1, NGAUSS FACTOR = 0.0 DO 501 I = 1, 8 FACTOR = FACTOR + X E ( I ) * E M 1 ( I , K ) EPLUS = 0.5D0 * FACTOR**2 EMINUS = EPLUS DO 505 I = 1 , 8 EPLUS = EPLUS + (EN 1 ( I , R ) - H * 0 . 5 D 0 * E M 2 ( I , K ) ) * X E ( I ) EMINUS = EMINUS + ( E N 1 ( I , K ) + H * 0 . 5 D 0 * E M 2 ( I , K ) ) * X E ( I ) 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 EPOS = EPLUS IF(EMINUS.GT.EPOS) EPOS=EMINUS ENEG = 0.0D0 GO TO 530 EPOS = EPLUS ENEG = EMINUS GO TO 510 EPOS = 0.0D0 ENEG = EPLUS IF (DABS(EMINUS).GT.DABS(ENEG)) ENEG = EMINUS GO TO 530 EPOS = EMINUS ENEG = EPLUS * FINDS THE POSITION OF THE NEUTRAL AXIS ESTR(1) = EMINUS F l ( 1 ) = -1.0D0 ESTR(NG2) = EPLUS F l ( N G 2 ) = 1.0D0 DO 512 L = 1, NGAUSS SUM = 0.5*FACTOR**2 DO 511 I = 1,8 SUM = SUM + ( E N 1 ( I , K ) - G A P ( L ) * H / 2 . 0 * E M 2 ( I , K ) ) * X E ( I ) ESTR(L+1) = SUM F l ( L + 1 ) = GAP(L) CONTINUE DO 515 I = 1, NG1 PROD = E S T R ( I ) * E S T R ( I + 1 ) IF (PROD.LE.0.0D0) GO TO 516 CONTINUE XN = F I ( I ) - E S T R ( I ) * ( F I ( I + 1 ) - F I ( I ) ) / ( E S T R ( I + 1 ) - E S T R ( I ) ) IF (ESTR(I).EQ.O.0DO) GO TO 518 IF ( E S T R ( I ) . L T . 0 . 0 D 0 ) HN = (1.0D0 - XN)*H/2.0D0 IF (ESTR(I).GT.0.0D0) HN = (1.0D0 + XN)*H/2.0D0 GO TO 520 IP (ESTR(I+1).LT.0.0D0) HN = (1.0D0 + XN)*H/2.0D0 IF ( E S T R d + 1 ) .GT.0.0D0) HN = ( 1 . 0D0 - XN)*H/2.0D0 SUME = SUME + (HN/H)*(E0*EPOS)**XKT*GAW(K)  90 530 538 540 550  560  563 564  565  377 4  8888 8889 8890 8810  5650 5655 833 7330 7331 8334  8338 8336  IF(EPOS.LT.EMAXP) GO TO 538 EMAXP = EPOS IF(DABS(ENEG).LT.DABS(EMAXN)) GO TO 540 EMAXN = ENEG CONTINUE CONTINUE SMAXP = E0*EMAXP SMAXN = E0*EMAXN I F (DABS(SMAXN).LE.FC) GO TO 560 SMAXN = SMAXN - ( ( E 0 + EN*E0)*EMAXN + F C * ( 1 . 0 + EN)) I F (SUME.EQ.0.ODO.OR.SMAXP.EQ.0.ODO) GO TO 563 SUME = SUME/(2.0*NELEM*(XKT+1.0)*SMAXP**XKT) F T T = FT * SUME**(-1.0D0/XKT) GO TO 564 F T T = FT I F (SMAXP.GE.FTT) GO TO 3774 DEFL = 0.0D0 DO 565 I E = 1, NELEM J = (IE-1)*NDOF + 3 I F ( D A B S ( X ( J ) ) . G T . D A B S ( D E F L ) ) DEFL = X ( J ) CONTINUE J = NEQ - 1 I F ( D A B S ( X ( J ) ) . G T . D A B S ( D E F L ) ) DEFL = X ( J ) I F (NP.EQ.O) PDEFL = FQ3 I F (NP.NE.0) PDEFL = P3 CONTINUE I F (NINT.EQ.0) GO TO 8810 I F (IERROR.EQ.1) WRITE(*,8888) I F (IERROR.EQ.0.AND.SMAXP.LT.FTT) WRITE(*,8889) SMAXP IF (IERROR.EQ.0.AND.SMAXP.GE.FTT) WRITE(*,8890) SMAXP FORMAT(' IERROR=1,FAILS (DIVERGENCE OR SINGULAR MATRIX)'/) FORMAT(' IERROR=0 SMAXP = ',E15.6,' SURVIVES'/) FORMAT(' IERROR=0 SMAXP = ',E15.6,' FAILS'/) CONTINUE I F (NP.EQ.O) GO TO 4500 I F (IERROR.EQ.1) GO TO 7330 I F (SMAXP.GT.FTT) GO TO 7331 I F (SMAXP.EQ.FTT) GO TO 7337 PI = P3 I F (SUME.EQ.0.ODO.OR.SMAXP.EQ.0.ODO) GO TO 5650 SMAX1 = SMAXP*SUME**(1.0D0/XKT) GO TO 5655 SMAX1 = SMAXP DO 833 J = 1, NEQ XP(J) = X(J) GO TO 8334 P2 = P3 GO TO 8334 P2 = P3 NFAIL = 1 SMAX2 = SMAXP*SUME**(1.0D0/XKT) I F (P1.EQ.0.ODO) GO TO 8338 TOLP = (P2-P1)/P1 IF (TOLP.LE.TOP) GO TO 7338 GO TO 8336 IF (P2.LE.0.1D0) GO TO 7338 IF (NFAIL.EQ.1) GO TO 8340 P3 = (PI + P2)/2.0  91 8340 7337  7338  7339 570  683 4500  4331 4580 4833 4330 4334 5338 4337  4338 4339 670  901  GO TO 3773 P3 = PI + (P2-P1)*(FT-SMAX1)/(SMAX2-SMAX1) GO TO 3773 P = P3 PP = P3 PAV = P3 GO TO 7339 I F (P1.EQ.0.ODO) P2 = 0.0D0 P = P2 PP = PI PAV = (P+PP)/2.0 CALL TIME(ZIM) ZIM = ZIM - ZIMO WRITE(2,570) PP,P,PAV,SMAXP,SMAXN,DEFL,PDEFL,SLAMDA 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 FORMAT(' TIME =',F7.1,' SECS.'/) WRITE(*,570) PP,P,PAV,SMAXP,DEFL,PDEFL,SLAMDA GO TO 901 I F (IERROR.EQ.1) GO TO 4330 I F (SMAXP.GT.FTT) GO TO 4330 I F (SMAXP.EQ.FTT) GO TO 4337 IF (NFLAG.EQ.1) GO TO 4331 FQ1 = FQ2 FQ2 = 2.0D0*FQ2 GO TO 4580 FQ1 = FQ3 DO 4833 J = 1,NEQ XP(J) = X(J) GO TO 4334 NFLAG = 1 FQ2 = FQ3 I F (FQ1.EQ.0.ODO) GO TO 5338 TOLP = (FQ2-FQ1)/FQ1 I F (TOLP.LE.TOP) GO TO 4338 I F (NFLAG.EQ.0) FQ3 = FQ2 IF (NFLAG.EQ.1) FQ3 = (FQ1+FQ2)/2.ODO GO TO 3773 P = FQ3 PP = FQ3 PAV = FQ3 GO TO 4339 P = FQ2 PP = FQ1 PAV = (P+PP)/2.0 CALL TIME(ZIM) ZIM = ZIM-ZIM0 WRI TE ( 2 , 67 0 ) PP , P , PAV, SMAXP, SMAXN ,'DEFL , PDEFL 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 CONTINUE  3131  C C  c c*  c 4  c 5  C 10  350 150  CONTINUE CLOSE (1,STATUS='KEEP') CLOSE (2,STATUS='KEEP') STOP END  92  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 GAP( 1 -0.774596669241483D0 0.0D0 GAP (2 -GAP(1) GAP (3 GAW( 1 0.555555555555556D0 GAW(2 0.888888888888889D0 GAW(1) GAW(3 GO TO 10 * * * 4 POINT GAUSSIAN INTEGRATION GAP( 1 = -0 861136311594053DO GAP (2 -0.339981043584856DO 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 *** 5 POINT GAUSSIAN INTEGRATION GAP( 1 = -0 906179845938664D0 GAP (2 -0.538469310105683DO GAP (3 0.0D0 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) I N I T I A L I S E S EN1,EM1,EM2 DO 150 I L = 1 8 . DO 350 IK = 1 NGAUSS EN 1 ( I L , IK) • • 0.0D0 EM1(IL,IK) • • 0.0D0 E M 2 ( I L , I K ) : 0.0D0 CONTINUE CONTINUE DO 250 I = 1 , NGAUSS EN 1(1 ,1 ) = (-0.75D0+0.75D0*GAP(I)**2)/DEL EN 1(2 ,1 ) = (-1.D0-2.D0*GAP(I)+3.D0*GAP(I)**2)*0.25D0 EN 1(5 ,1 ) = (0.75D0-0.75D0*GAP(I)**2)/DEL EN 1 ( 6 , 1 ) = (-1.D0+2.D0*GAP(I)+3.D0*GAP(l)**2)*0.2 5D0 E M I ( 3 , 1 ) = (-0.7 5DO+0.7 5DO*GAP(I)**2)/DEL E M I ( 4 , 1 ) = (-1.D0-2.D0*GAP(I)+3.D0*GAP(l)**2)*0.25D0  93  250  EM1(7,1 EM1(8,I EM2(3,1 EM2(4,1 EM2(7,1 EM2(8,I 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 I F ( A A ( 1 ) . L E . 0 . D 0 ) 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 = A A ( I J D ) KO = 1 I F ( J . G T . L H B ) 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 I F ( I ERROR.EQ.1) RETURN A A ( I J D ) = DSQRT(SUM) DO 568 I = 1, KB II = J+I KO = 1 IF (II.GT.LHB) KO=II-KB " SUM = AA(IJD+I) I F ( I . E Q . K B ) 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 A A ( 6 7 2 ) , S ( 8 4 )  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 =  675 680 C 685  690 699  1-1  KO = 1 I F ( I . G T . L H B ) KO=I-KB SUM = S ( I ) II = LHB*I-KB DO 67 5 K = KO, 11 IK = KB*K+I-KB SUM = SUM-AA(IK)*S(K) S ( I ) = SUM/AA(II) CONTINUE BACKWARD SUBSTITUTION 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 ( N I ) II = LHB*NI-KB DO 690 K = I 1, KO IK = KB*NI+K-KB SUM = SUM-AA(IK)*S(K) S(NI) = SUM/AA(II) CONTINUE RETURN END  C SUBROUTINE CONVRG(XO,X,IER,NEQ,EPSLON,ITER) THIS SUBROUTINE CHECKS THE CONVERGENCE OF SOLUTION VECTOR IMPLICIT REAL*8(A-H,0-Y) COMMON/C2/DIFP,NINT DIMENSION X O ( 8 4 ) , X ( 8 4 ) IER = 0 PARXO = O.0D0 PARDIF = 0.0D0 PARX = 0.0D0 DO 602 1 = 1 , NEQ PARXO = PARXO + X O ( I ) * * 2 PARX = PARX + X ( I ) * * 2 602 PARDIF = PARDIF + ( X ( I ) - X O ( 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  C* C  95  604 605  RETURN RETURN 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  FILES.  SAMPLE INPUT DATA FILE FOR AXIAL COMPRESSION 10 5 1 10 1 0 0.0 -0.001 2 1 2 1 3 11 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)  100  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 S I Z E EFFECTS CONSIDERED) 15 A u g u s t , 1987  A PROGRAM FOR THE EVALUATION OF THE R E I A B I L I T Y 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 L I M I T . END LOAD I S APPLIED CENTRALLY. LATERAL LOADS CAN BE DISTRIBUTED ALONG THE LENGTH OF THE MEMBER THE PROGRAM FINDS THE R E L I A B I L I T Y INDEX BETA FOR A * BEAM-COLUMN TAKING INTO ACCOUNT 5 RANDOM VARIABLES * WHICH CONSTITUTE THE LOAD AND MATERIAL RESISTANCES * PROBLEM DATA I S READ FROM UNIT OUTPUT IS STORED I N UNIT #2.  * * * * *********************************************************  #1  * MAXIMUM OF 10 VARIABLES * * MAXIMUM OF 20 ELEMENTS * I M P L I C I T REAL*8(A-H,0-Z) REAL*8 INVNPR,NORMPR DIMENSION X ( 1 0 ) , Y ( 1 0 ) , U ( 1 0 ) , D E L T A ( 1 0 ) , S I G ( 1 0 ) 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 L O C ( 1 0 ) , M U ( 1 0 ) , N N ( 1 0 ) , NNN(10) INTEGER*2 ICODEOO) INTEGER*2 M X C ( 1 0 ) , 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  *  c c c c c c c c c c  * * * * * * * * * *  Q  * * *  DEFINE  VARIABLES  *  **********************************************************  FCN W H GAMA 1 ALFD ALFL EMIN NELEM = NJBC NGAUSS =  COMPRESSIVE STRENGTH, FIFTH PERCENTILE WIDTH OF SECTION DEPTH OF SECTION RATIO OF NOMINAL DL TO LL DEAD LOAD FACTOR L I V E LOAD FACTOR MODULUS OF E L A S T I C I T Y , MEAN VALUE NO OF ELEMENTS NO OF JOINTS WITH B.C. NO OF INTEGRATION POINTS  * * * * * * * * * *  101  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * NO OF B.C. AT NODE I = * B.C. CODE * = U = * UX * = W 3 = * WX 4 N = NO OF RANDOM VARIABLES FOR TOTAL PROB. * * EN1 ,EM1,EM2 = INTERPOLATION FUNCTIONS = CORDINATE AT GAUSS POINT * GAP = CORRESPONDING WEIGHT GAW * = NO OF * NELEM ELEMENTS = NO OF GAUSS POINTS * NGAUSS = MAX. * NITER NO OF ITERATIONS = TOLERANCE FOR LOAD • * TOP = TOLERANCE FOR SOLUTION VECTOR * EPSLON = MATERIAL STRENGTH IN COMPRESSION * FC = MATERIAL STRENGTH IN TENSION * FT = MOE OF THE MATRIAL * EO = SLOPE OF THE STRESS-STRAIN CURVE * EN = MEMBER LENGTH * SPAN = WIDTH OF SECTION W * = DEPTH OF SECTION * H = ECCENTRICITY OF AXIAL LOAD * E = NO OF EQUATIONS TO BE SOLVED * NEQ = NO OF NODES NJOINT * = NO OF VARIABLES PER NODE * NDOF = NO OF NODES PER ELEMENT * NODEL * SREF = REFERENCE SPAN = SIZE EFFECT SHAPE PARAMETER (COMP .) XKC = SIZE EFFECT SHAPE PARAMETER (TENS .) XKT = NO OF VARIABLES PER NODE NDIMB * LBW,LHB = HALF BANDWIDTH INCLUD. THE DIAG. * = NO OF UNKNOWNS FOR TOTAL PROBLEM * NA ***************************** *•* ***************************  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  Q C  NBC(I) = IX 1 2  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,*) ( I X ( N J , J V ) , J V = 1 , N B C ( N J ) ) 992 2 CONTINUE  C READ(1,*) N, ( I C O D E ( I ) , I = READ(1,*) (MXC(I), I = 1,N) DO 7779 I = 1,N MEX(I) = 0 IF (MXC(I).EQ.0) GO TO 7779  1,N)  7779 7780 7784 7782  GO TO 7780 CONTINUE GO TO 7782 WRITE(*,7784) FORMAT(' ENTER EXPONENTS FOR DISTRIBUTION OF EXTREMES'/) READ(*,*) ( M E X ( I ) , 1=1,N) 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,*) A V E R ( I C ) , S T D ( I C ) GO TO 9 C C LOGNORMAL (CODE=2) C 12 READ(1,*) A V E R ( I C ) , S T D ( I C ) GO TO 9 C C WEIBULL ( CODE=3 ) C 13 READ(1,*) L O C ( I C ) , S C A L E ( I C ) , S H A P E ( I C ) GO TO 9 C C GUMBEL EXTREME TYPE I ( CODE=4 ) C 14 READ(1,*) B ( I C ) , A ( I C ) 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 ) I F ( I C O D E ( I ) . N E . 3 ) GO TO 805 IF (X(I ) .GT.LOCd ) ) GO TO 805 WRITE(*,1270) 1270 FORMAT (' CHANGE I N I T I A L 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)  103  C C C C  2080 3080 4080 C C  CALC VECTORS EN1, EM1, CALL SHAP(DELT)  EM2  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 IF (CC .GT. CA) GO TO 3080 CK = CK1 GO TO 4080 CK = CK2 CONTINUE OBTAIN NOMINAL DESIGN LOAD PLN = R*W*H*FCN*CK/(ALFD*GAMA1+ALFL)  C W R I T E ( 2 , 1 0 8 0 ) ( I C O D E ( 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 F 1 X ( X ) , AND F 2 X ( 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 ( F 2 X ( I ) . L E . 0 . 0 D 0 ) GO TO 68 DSIG = DLOG(CONST) - Y ( I ) * Y ( I ) * 0 . 5 D 0 - D L O G ( F 2 X ( I ) ) IF (DSIG.LT.-709.0D0) GO TO 865 S I G ( I ) = DEXP(DSIG) GO TO 87 865 S I G ( 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 + S I G ( I ) * S I G ( I ) * D E L T A ( I ) * D E L T A ( I ) SUM = DSQRT(SUM) DO 20 I = 1 ,N NN(I) = - S I G ( I ) * D E L T A ( I ) / S U M 20 NNN(I) = DABS(NN(I)) C C CALC BETA SDMU=0.0D0 SDX = 0.0D0 DO 25 1=1,N SDMU = SDMU + D E L T A ( I ) * M U ( I ) 25 SDX = SDX + D E L T A ( 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 = S I G ( 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 ) = S I G ( 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,' f p = ' , E l 3 . 6 / ) 51 FORMAT(5X,'BETA = ',F10.3) GO TO 900 65 I F (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 I T S ' , / , 2 ' LOCATION PARAMETER.',/, 3 ' TRY NEW I N I T I A L POINT') 1341 FORMAT(' I ERR 1 =',12,' NEGATIVE OR ZERO DENSITY F 2 X ( l ) ' , / , 1 ' TRY NEW I N I T I A L 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 S C ( N ) , S K ( N ) , A ( N ) , B ( N ) , X ( N ) , A V E R ( 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) F 1 X ( I ) = NORMPR(RATIO) F2X(I)=DEXP(-0.5D0*RATIO*RATIO)/(STD(l)*DSQRT(PI 2 ) ) I F (MXC(I).EQ.0) GO TO 20 CALL EXTR (F 1X (I ) , F2X (I ) , MXC (I ) , M E X ( 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 ) . L E . 0 . 0 ) GO TO 99 PARAM = (DLOG(X(I)) - MU)/SDP F 1 X ( I ) = NORMPR(PARAM) POW = DEXP(-0.5D0*PARAM*PARAM) F 2 X ( I ) = POW/(SDP*X(I)*DSQRT(PI2)) IF (MXC(I).EQ.O) GO TO 2 0 CALL E X T R ( F I X d ) . , F 2 X ( l ) , MXC (I ) , M E X ( I ) ) GO TO 20  106 C WEIBULL C 13 IF ( X ( I ) . L E • L O C ( I ) ) GO TO 99 POW = - ( ( X ( I ) - L O C ( l ) ) / S C ( I ) ) * * S K ( I ) POW = DEXP(POW) F 1 X ( 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 IF (MXC(I).EQ.0) GO TO 20 CALL E X T R ( F 1 X ( I ) . F 2 X ( I ) M X C ( I ) , M E X ( I ) ) GO TO 20 C C GUMBEL EXTREME TYPE I C 14 POW = - A ( I ) * ( X ( I ) - B(I ) ) POW = DEXP(POW) F 1 X ( I ) = DEXP(-POW) F 2 X ( I ) = A(I)*POW * F 1 X ( I ) IF (MXC(I).EQ.0) GO TO 20 CALL E X T R ( F 1 X ( I ) F 2 X ( I ) , M X C ( I ) , M E X ( I ) ) 20 CONTINUE RETURN 99 IERR=1 RETURN END  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 FUNCTION NORMPR(X) * NORMAL PROBABILITY INTEGRAL IMPLICIT REAL*8(A-H,0-Z) REAL*8 NORMPR DIMENSION E ( 1 6 ) , H ( 1 6 ) PI = 2.ODO * DSQRT(DATAN(1.0D0)) IF (DABS(X).GT.5.0D0) GO TO 20 0.989400934991650EO E( 1 0.944575023073233EO E(2 0.865631202387832E0 E(3 0.755404408355003E0 E(4 0.617876244402644E0 E(5 0.458016777657227E0 E(6 0.281603550779259E0 E(7 0.095012509837637E0 E(8 0.027152459411754E0 H(1 0.062253523938648E0 H(2 0.095158511682493E0 H(3 0.124628971255534E0 H(4 0. 1 49595988816577E0 H(5  107  1  10  20  25  H(6) = 0.169156519395003E0 H(7) = 0.182603415044924E0 H(8) = 0.189450610455068E0 DO I I = 1,8 E{ 1 7-1 ) = - E ( I ) 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) CONTINUE ERF = Y * S/PI NORMPR = (1.0D0 + ERF)/2.0D0 RETURN I F (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 I F (X.LT.0.0D0) NORMPR = S/2.ODO RETURN I F (X.GT.0.0D0) NORMPR = 1.ODO IF (X.LT.0.0D0) NORMPR = 0.ODO RETURN END  C C C C  5  80 20  FUNCTION INVNPR(Y) * INVERSE NORMAL PROBABILITY * 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 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 INVNPR = 0.0 RETURN INVNPR=X2 RETURN END  C SUBROUTINE COLUMN(XW,N,PAV) IMPLICIT REAL*8(A-H,0-Z) DIMENSION F ( 8 ) , T K 0 ( 6 7 2 ) , X E ( 8 ) , X W ( N ) 1,R(84),X0(84),X(84),B(84),B1(8,8),B2(8,8),B3(8,8),B4(8,8)  108  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  43 44  2, B 5 ( 8 , 8 ) , B 6 ( 8 , 8 ) , B 7 ( 8 , 8 ) , B 8 ( 8 , 8 ) , B 9 ( 8 , 8 ) , Y ( 5 ) , R E ( 8 ) , X P ( 8 4 ) 3, Q ( 2 0 ) , I Q ( 2 0 ) , E S T R ( 7 ) , F I ( 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)  ********************************************** * EN I,EM1 ,EM2 = INTERPOLATION FUNCTIONS * GAP = CORDINATE AT GAUSS POINT * GAW = CORRESPONDING WEIGHT * = NO OF ELEMENTS NELEM = NO OF GAUSS POINTS * NGAUSS * NITER = MAX. NO OF ITERATIONS = TOLERANCE FOR LOAD * TOP = TOLERANCE FOR SOLUTION VECTOR * EPSLON = MATERIAL STRENGTH IN COMPRESSION * FC * = MATERIAL STRENGTH IN TENSION FT * EO = MOE OF THE MATRIAL = SLOPE OF THE STRESS-STRAIN CURVE * EN = MEMBER LENGTH * SPAN * W = WIDTH OF SECTION = DEPTH OF SECTION * H * E = ECCENTRICITY OF AXIAL LOAD * = NO OF EQUATIONS TO BE SOLVED NEQ = NO OF NODES * NJOINT = NO OF VARIABLES PER NODE * NDOF = NO OF NODES PER ELEMENT * NODEL = REFERENCE SPAN SREF * = S I Z E EFFECT SHAPE PARAMETER (COMP. ) * XKC = S I Z E EFFECT SHAPE PARAMETER * XKT (TENS. ) * NDIMB = NO OF VARIABLES PER NODE * LBW,LHB = HALF BANDWIDTH INCLUD. THE DIAG. = NO OF UNKNOWNS FOR TOTAL PROBLEM * NA ******************************************************** 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 READ(1,*) I Q ( I ) , Q ( l ) 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  761  760 792 C C 3773  4000 4001 C 80 C C C  81  * 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 FQ1 = 0.0D0 FQ2 = 1.0D0 FQ3 = FQ2 NFLAG = 0 DO 792 J = 1, NEQ X P ( J ) = O.ODO START CALCULATIONS FOR TRIAL LOAD LEVELS 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 FORMAT(//' SOLUTION FOR P =',E15.6,' :'/) FORMAT(//* SOLUTION FOR LATERAL LOAD FACTOR= , E1 5. 6 , ':'/) I N I T I A L I S E ARRAYS DO 80 J = 1, NEQ XO(J) = X P ( J ) R ( J ) = 0.D0 1  EXTERNAL LOAD VECTOR R RE = ELEMENT LOAD VECTOR IF (QO.EQ.O.ODO) GO TO 87 DO 81 J = 1 , 8 RE(J)=0.DO 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) + R E ( J J )  110 82 83 87  180 185  C C  CONTINUE CONTINUE I F (NQ.EQ.O) GO TO 185 DO 180 J = 1 , NQ JS = ( I Q ( J ) - 1 ) * N D O F + 3 R(JS) = R(JS) + Q(J)*FQ 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  BEGIN ITERATIONS AT THE TRIAL LOAD LEVEL CONTINUE DO 84 I = 1 , NA 84 TKO(I) = O.ODO DO 85 K = 1, NEQ 85 B(K) = -R(K) DO 645 I E = 1, NELEM C I N I T I A L I Z E ARRAYS DO 88 I = 1 , 8 F ( I ) = O.ODO DO 86 J = 1 , I B l ( I , J ) == 0 .0D0 B 2 ( I , J ) == 0 .0D0 B 3 ( I , J ) == 0 .0D0 B 4 ( I , J ) == 0 .0D0 B 5 ( I , J ) == 0 .0D0 B 6 ( I , J ) == 0 .0D0 B 7 ( I , J ) == 0 .0D0 B 8 ( I , J ) == 0 .0D0 B 9 ( I , J ) == 0 .0D0 86 CONTINUE 88 CONTINUE C PICK ELEMENT SOLUTION FROM GLOBAL VECTOR DO 90 J J = 1 , 8 K = ( I E - 1)*NDOF + J J X E ( J J ) = XO(K) 90 CONTINUE DO 101 K = 1, NGAUSS Y(K) = 0.D0 DO 91 1=1 , 8 Y(K) = Y ( K ) + X E ( I ) * E M 1 ( 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 777  111 DO 100 L = 1 , NGAUSS 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.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 B 5 ( I , J ) = B5(I,J)+DEL*0.5D0*AR*(EN 1 ( I , K ) - G A P ( 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) B 6 ( I , J ) = B6(I,J)+DEL*0.5D0*AR*(EN 1 ( I , K ) - G A P ( L ) * 1 H*0.5D0*EM2(I,K))*(E0+EN*E0)*FAC*Y(K)*EM1(J,K)* 2 GAW(K)*GAW(L) B 7 ( 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) B 8 ( 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) B 9 ( 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 = B 1 ( I , J ) + B 2 ( I , J ) + B 3 ( I , J ) + B 4 ( 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 T K O ( I J ) = 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 + I X ( U O , J ) LBW1 = LBW - 1 DO 108 K = 1, LBW1 J J = II - LBW + K IF ( J J . L E . 0 ) GO TO 1080 I J = ( J J - 1 ) * ( L B W - 1 ) + II T K O ( I J ) = 0.0D0 1080 J J = 11 + K IF (JJ.GT.NEQ) GO TO 108  C  1  108  110 111 C C  I J = {11- 1)*(LBW-1) + J J T K O ( I J ) = O.ODO CONTINUE I J = ( I I - 1)*(LBW-1) + I I T K O ( I J ) = 1.0DO B ( I I ) = O.ODO i CONTINUE CONTINUE  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 I F (ITER.EQ.NITER) GO TO 431 I F (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 I E = 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 + X E ( I ) * E M 1 ( I , K ) EPLUS = 0.5D0 * FACTOR**2 EMINUS = EPLUS DO 505 I = 1, 8 EPLUS = EPLUS + (EN 1 ( I , K ) - H * 0 . 5 D 0 * E M 2 ( I , K ) ) * X E ( I ) EMINUS = EMINUS + (EN 1 ( I , K ) + H * 0 . 5 D 0 * E M 2 ( I , K ) ) * X E ( 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  11 3  507 508  509 C C 510  511 512  515 516  518 520 530 538 540 550  560  563 564  565  GO TO 530 EPOS = EPLUS ENEG = EMINUS GO TO 510 EPOS = O.ODO ENEG = EPLUS I F (DABS (EMINUS ) .GT.'DABS (ENEG) ) ENEG = EMINUS GO TO 530 EPOS = EMINUS ENEG = EPLUS * FINDS THE POSITION OF THE NEUTRAL AXIS ESTR(1) = EMINUS F l ( 1 ) = -1.0D0 ESTR(NG2) = EPLUS F l ( N G 2 ) = 1.0D0 DO 512 L = 1, NGAUSS SUM = 0.5*FACTOR**2 DO 511 I = 1,8 SUM = SUM + (EN1(I,K) - G A P ( L ) * H / 2 . 0 * E M 2 ( I , K ) ) * X E ( I ) ESTR(L+1) = SUM F l ( L + 1 ) = GAP(L) CONTINUE DO 515 I = 1, NG1 PROD = E S T R ( I ) * E S T R ( I + 1 ) I F (PROD.LE.O.ODO) GO TO 516 CONTINUE XN = F I ( I ) - E S T R ( I ) * ( F I ( 1 + 1 ) - F I ( I ) ) / ( E S T R ( I + 1 ) - E S T R ( l ) ) I F (ESTR(I).EQ.O.ODO) GO TO 518 I F (ESTR(I).LT.O.ODO) HN = (1.0D0 - XN)*H/2.0D0 I F (ESTR(I).GT.O.ODO) HN = (1.0D0 + XN)*H/2.0D0 GO TO 520 I F (ESTR(1+1).LT.O.ODO) HN = (1.0D0 + XN)*H/2.0D0 I F (ESTR(I+1).GT.O.ODO) HN = (1.0D0 - XN)*H/2.0D0 SUME = SUME + (HN/H)*(E0*EPOS)**XKT*GAW(K) IF(EPOS.LT.EMAXP) GO TO 538 EMAXP = EPOS IF(DABS(ENEG).LT.DABS(EMAXN)) GO TO 540 EMAXN = ENEG CONTINUE CONTINUE SMAXP = E0*EMAXP SMAXN = E0*EMAXN IF (DABS(SMAXN).LE.FC) GO TO 560 SMAXN = SMAXN - ( ( E 0 + EN*E0)*EMAXN + FC*(1.0 + E N ) ) I F (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 FTT = FT I F (SMAXP.GE.FTT) GO TO 3774 DEFL = O.ODO DO 565 I E = 1, NELEM J = (IE-1)*NDOF + 3 IF ( D A B S ( X ( J ) ) . G T . D A B S ( D E F L ) ) DEFL = X ( J ) CONTINUE J = NEQ - 1 IF ( D A B S ( X ( J ) ) . G T . D A B S ( D E F L ) ) DEFL = X ( J )  11 4  3774  8888 8889 8890 8810  5650 5655 833 7330 7331 8334  8338 8336 8340 7337  7338  7339 4500  4331 4580 4833  IF (NP.EQ.O) PDEFL = FQ3 I F (NP.NE.O) PDEFL = P3 CONTINUE IF (NINT.EQ.0) GO TO 8810 IF (IERROR.EQ.1) WRITE(*,8888) IF (I ERROR. EQ.0.A.ND. SMAXP. L T . FTT) WRITE (*, 8889 ) SMAXP IF (IERROR.EQ.0.AND.SMAXP.GE.FTT) WRITE(*,8890) SMAXP FORMAT(' IERROR = 1,FAILS (DIVERGENCE OR SINGULAR MATRIX)'/] FORMAT(' IERROR = 0 SMAXP = *,E15.6,' SURVIVES'/) FORMAT(' IERROR = 0 SMAXP = ',E15.6,' FAILS'/) CONTINUE IF (NP.EQ.O) GO TO 4500 IF (IERROR.EQ.1) GO TO 7330 I F (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 SMAX1 = SMAXP DO 833 J = 1, NEQ XP(J) = X ( J ) GO TO 8334 P2 = P3 GO TO 8334 P2 = P3 NFAIL = 1 SMAX2 = SMAXP*SUME**(1.0D0/XKT) I F (PI.EQ.0.ODO) GO TO 8338 TOLP = (P2-P1)/P1 IF (TOLP.LE.TOP) GO TO 7338 GO TO 8336 I F (P2.LE.0.1D0) GO TO 7338 I F (NFAIL.EQ.1) GO TO 8340 P3 = (PI + P2)/2.0 GO TO 3773 P3 = P1 + (P2-P1)*(FT-SMAX1)/(SMAX2-SMAX1) GO TO 3773 P = P3 PP = P3 PAV = P3 GO TO 7339 I F (PI.EQ.0.ODO) P2 = 0.ODO P = P2 PP = PI PAV = (P+PP)/2.0 CONTINUE GO TO 901 I F (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 FQ1 = FQ3 DO 4833 J = 1,NEQ XP(J) = X ( J )  4330 4334 5338 4337  4338 4339 901  GO TO 4334 NFLAG = 1 FQ2 = FQ3 IF (FQ1.EQ.O.ODO) GO TO 5338 TOLP = (FQ2-FQ1)/FQ1 IF (TOLP.LE.TOP) GO TO 4338 IF (NFLAG.EQ.O) FQ3 = FQ2 IF (NFLAG.EQ.1) FQ3 = (FQ1+FQ2)/2.0D0 GO TO 3773 , P = FQ3 PP = FQ3 PAV = FQ3 GO TO 4339 P = FQ2 PP = FQ1 PAV = (P+PP)/2.0 CONTINUE RETURN END  C C*  SUBROUTINE SHAP(DELT) 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 *** 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 *** 4 POINT GAUSSIAN INTEGRATION 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 *** 5 POINT GAUSSIAN INTEGRATION 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) I N I T I A L I S E S EN 1 ,EM1 ,EM2 DO 150 I L = 1, 8 DO 350 IK = 1, NGAUSS -  C  C 4  C 5  C 10  116  350 150  250  EN 1 ( I L , I K ) = O.ODO E M 1 ( I L , I K ) = O.ODO E M 2 ( I L , I K ) = O.ODO CONTINUE 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 E N 1 ( 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) CONTINUE RETURN END  C SUBROUTINE DECOMP(NN,LHB,AA,IERROR) THIS SUBROUTINE DECOMPOSES A MATRIX USING CHOLESKY 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 I F ( A A ( 1 ) . L E . 0 . D 0 ) 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 = A A ( I J D ) KO = 1 I F ( J . G T . L H B ) 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 I F ( I ERROR.EQ.1) RETURN A A ( I J D ) = DSQRT(SUM) DO 568 I = 1, KB II = J + I KO = 1 IF (II.GT.LHB) KO=II-KB SUM = AA(IJD+I) I F ( I . E Q . K B ) GO TO 565 DO 540 K = KO, J1 JK = KB*K+J-KB IK = KB*K+II-KB C* C  540 565 568 590  SUM = SUM-AA(IK)*AA(JK) AA(IJD+I) = SUM/AA(IJD) CONTINUE 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 A A ( 6 7 2 ) , S ( 8 4 ) C FORWARD SUBSTITUTION KB = LHB-1 S(1) = S ( 1 ) / A A ( 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 ) I I = 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 ( N I ) II = LHB*NI-KB DO 690 K = I 1, KO IK = KB*NI+K-KB 690 SUM = SUM-AA(IK)*S(K) S ( N I ) = 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 X O ( 8 4 ) , X ( 8 4 ) IER = 0 PARXO = O.ODO PARDIF = O.ODO PARX = O.ODO DO 602 I = 1, NEQ PARXO = PARXO + X O ( I ) * * 2 PARX = PARX + X ( I ) * * 2  1 18  602 1002 606  603 604 605  PARDIF = PARDIF + ( X ( I ) - X O ( I ) ) * * 2 IF (NINT.EQ.1) WRITE(*,1002) PARXO, PARX, PARDIF FORMAT(' NORMX0=',E13.6,'NORMX=',E13.6,'NORMDIF=',E13.6/) I F (ITER.EQ.0) GO TO 606 IF (PARDIF.GE.DIFP) GO TO 605 DIFP = PARDIF I F (PARX0.EQ.0.ODO) GO TO 603 DIF = DSQRT(PARDIF/PARX0) I F (DIF.LE.EPSLON) GO TO 604 IER = 1 RETURN RETURN IER = 2 RETURN END 1  C  6644  1202  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 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 D E L T A ( I ) = (PU1 - PU2)/(0.02D0*XW(I)) IF ( l . G E . ( N - 2 ) ) GO TO 1202 GO TO 6644 DELTA(N-1) = -PLN*GAMA1 DELTA(N) = -PLN RETURN END  SAMPLE INPUT/OUTPUT FILE  1 20  SAMPLE  INPUT DATA F I L E FOR R E L I A B I L I T Y 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 11 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 CODES : 3 3 BETA = 5.136 ITERATIONS = 4 TOLB = 0.0100 VECTOR XO VECTOR X S E N S I T I V I T Y COEFFS. L = 3.2 </> = 0. 6  OUTPUT F I L E  4693.6 3878.8 0.8282  7.053 32.302 0.0000  8.538 30.358 0.0000  1.025 1.3053 0.3963  0.881 1.0553 0.3963  p  EXPLANATIONS Vector  Xo : I n i t i a l  Vector  X  value  f o r the  : C o o r d i n a t e s of t h e most point (design point)  Sensitivity  Solution  (trial)  coefficients  corresponding to  likely  variables. failure  S e n s i t i v i t y of (3 t o e a c h of the v a r i a b l e s . In t h i s c a s e j3 i s most s e n s i t i v e t o 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 ) . 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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062671/manifest

Comment

Related Items