NON-LINEAR ANALYSIS FOR TRANSVERSELY POST- TENSIONED TIMBER BRIDGE By Cai,Shen B.Eng. & M.Eng. Beijing Institute of Civil Engineering and Architecture A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1991 © Cai, Shen, 1992 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of Civil Engineering The University of British Columbia Vancouver, Canada Date ^ DE-6 (2/88) ABSTRACT Timber bridges have been very important in North America due to an abundant resource of the material and the relatively unsophisticated requirement for equipment and skilled labour. The transversely post-tensioned laminating method can improve the loading distribution capacity of the timber bridge structure and elongate the bridge service life. To develop a method for its analysis it is necessary to consider the non-linear behaviour of between-beams friction in the bridge structure. A method and corresponding computer program (PTB) has been developed for the non-linear analysis for the transversely post-tensioned timber bridge. The finite strip element method is used in the analysis. Wheel loadings are idealized as patches with static uniformly distributed loadings. The beam to beam friction parameters were obtained from tests. The relative movements between beams, in three-directions, stresses and deformations of the structure are obtained by the PTB program. As an application, a bridge with post-tensioned T-beams made of Parallam is considered. ii Table of Contents ABSTRACT^ ii List of Tables^ vi List of Figures^ vii Acknowledgement 1 INTRODUCTION 1 1.1 Background ^ 1 1.2 Literature Review ^ 3 1.3 2 Objectives of the Thesis ^ 4 NON-LINEAR FINITE ELEMENT ANALYSIS 6 2.1 6 Structural Model ^ 2.1.1^Middle surface displacement functions for the flange of the T-beam 6 2.1.2^Displacement functions for the centroid of the T-beam web 9 2.1.3^Degrees of freedom for the joint element 2.2 2.3 ^ 11 Strain Energy of the Structural System and Virtual Work ^ 13 2.2.1^Strain energy in the flange of the T-beam element ^ 14 2.2.2^Strain energy in the web ^ 18 2.2.3^Strain energy in the web-flange connector ^ 20 2.2.4^Virtual work done by non-conservative forces in the joint element 23 Global Equation ^ iii 26 2.3.1^Stiffness contribution from one T-beam element ^27 2.3.2^Stiffness contribution from one joint element ^29 2.3.3^Consistent load vectors ^30 2.3.4 System of equations corresponding to one T-beam (Element Matrix) 32 2.4 Global Equations and Solution Method ^ 2.4.1^Characteristics of the global equation 33 ^33 2.4.2 The incremental procedure ^ 34 2.4.3 Newton-Raphson method ^ 35 2.4.4 Jacobi iteration method ^39 3 FRICTION PARAMETERS FOR THE SPRING MODEL ^40 3.1 The Friction Test ^40 43 3.2 Spring Model ^ 49 ANALYTICAL RESULTS^ 4.1 Beam Test ^ 49 52 4.2 T-beam Analysis^ 4.2.1^Analysis for post-tensioning force only ^54 4.2.2^Analysis for the one patch of vertical loading ^55 4.2.3^Wheel loading analysis ^61 5 CONCLUSIONS Bibliography Appendices ^ 73 ^ 74 ^ 76 A SHAPE FUNCTIONS AND DERIVATIVES iv ^ 76 B THE PTB PROGRAM^ 79 B.1 PTB Program Specifications ^ 79 B.2 Program Structure ^ 80 C PTB USER'S MANUAL^ 85 D PTB SOURCE CODE AND I/O FILE ^ 92 D.1 PTB.FOR ^ 93 D.2 PTB.DAT ^ 147 D.3 PTB.OUT ^ 148 D.4 PTBB.OUT ^ 151 List of Tables 3.1 Data of Friction Test Samples ^ 3.2 Statistical Results of the Friction Test 4.1 Properties of the Cantilever Beam ^ 50 4.2 Status of the Springs in the Cantilever Beams ^ 50 4.3 Properties of T-beam Structure ^ 53 4.4 Analytical Results: Post-tensioning Force Only ^ 55 vi ^ 42 43 List of Figures 1.1 Analytical Structural Model ^ 5 2.1 The Structural Model ^ 7 2.2 The T-beam Model and The Joint Model ^ 7 2.3 The Distribution of the Discrete Springs ^ 12 2.4 The Constitutive Relationship of the Springs ^ 12 2.5 The Vechicle and Post-tensioning Load ^ 30 2.6 The Global Equation System ^ 34 2.7 The Incremental Method ^ 36 2.8 Newton-Raphson Method ^ 36 3.1 The Friction Test Setup ^ 41 3.2 Two Kinds of Friction Test ^ 41 3.3 Friction Test Recording Curves(1) ^ 44 3.4 Friction Test Recording Curves(2) ^ 45 3.5 Friction Test Recording Curves(3) ^ 46 3.6 Idealized Friction Test Curves ^ 47 3.7 Same Slip Assumption ^ 48 4.1 Two Cantilever Beam Structure ^ 50 4.2 Cantilever Beam Structure Displacements 4.3 Built-in Solid Beams ^ 51 4.4 Three T-beam Structure ^ 52 vi i ^ 51 4.5 Three T-beam Structure with Post-Tensioning Force Only ^ 54 4.6 Displacements of T-beam Structure under Post-Tensioning Force Only^56 4.7 T-beam Structure One Patch Loading Case ^ 57 4.8 Bending Stress of T-beam Structure One Patch Loading Case ^ 58 4.9 Displacements of T-beam Structure One Patch Loading Case ^ 58 4.10 Relative Displ. of T-beam Structure One Patch Loading Case ^ 59 4.11 Upward Deflection of the T-beam ^ 60 4.12 Displacements of T-beam under One Patch Loading ^ 60 4.13 T-beam Structure CASE 1 ^ 62 4.14 Bending Stress in T-beam Structure CASE 1 ^ 63 4.15 Displacements in T-beam Structure CASE 1 ^ 63 4.16 Relative Displacements in T-beam Structure CASE 1 ^ 64 4.17 Bending Stress under Moving Loading ^ 64 4.18 Displacements under Moving Loading ^ 65 4.19 T-beam Structure CASE 2 ^ 65 4.20 Bending Stress in T-beam Structure CASE 2 ^ 66 4.21 Displacements in T-beam Structure CASE 2 ^ 66 4.22 Relative Displacements in T-beam Structure CASE 2 ^ 67 4.23 The Effect of the Stiffness of the Spring ^68 4.24 T-beam Structure CASE 3 ^ 69 4.25 Bending Stress in T-beam Structure CASE 3 ^ 69 4.26 Displacements in T-beam Structure CASE 3 ^ 70 4.27 Relative Displacements in T-beam Structure CASE 3 ^ 70 4.28 Effect of the Fourier Terms on Bending Stress 4.29 Effect of the Fourier Terms on Displacement ^ ^71 71 4.30 Effect of the Fourier Terms on Relative Displacement in Joint 1 ^ 72 viii 4.31 Effect of the Fourier Terms on Relative Displacement in Joint 2 ^ 72 B.1 PTB Program Flow Chart ^ ix 81 Acknowledgement The author takes this opportunity to gratefully acknowledge the guidance of his advisor Dr. Ricardo 0. Foschi throughout this research and thesis preparation. His suggestions and advice during our discussions proved invaluable. The author also expresses his sincere gratitude to Dr. D.L. Anderson and Dr. H.G.L. Prion for reading this thesis and suggesting improvements. The financial support in the form of a Research Assistantship from the National Engineering Research Council of Canada is gratefully acknowledged. Also acknowledged are contributions by MacMillan Bloedel Research Limited Centre such as the use of their laboratories in the friction testing of Para1] am, and for providing other data on this composite product. Chapter 1 INTRODUCTION 1.1 Background Timber has been very important in North American bridge construction due to its abundance and the relatively unsophisticated requirements for equipment and skilled labour in comparison to steel and concrete bridges. More than a thousand timber bridges with short and medium spans have been built across Canada. Two kinds of bridges, using either nailed-laminated or glue-laminated construction have been most commonly designed. In the first half of this century, with the high quality and high strength of steel and reinforced concrete materials being used more and more, timber as a material for bridge construction was progressively ignored. Structural engineers were reluctant to use timber in bridges, since the frequent failure of such bridges in the past convinced engineers that these structures could not be expected to last more than 40 years [1], leading to expensive replacements. In 1973, the Ontario Ministry of Transportation and Communications began a laminated timber bridge testing program to evaluate the load carrying capacity of this kind of bridges. The results pointed out that the load carrying capacity and service life were a direct function of the 'tightness' of the structural system. 'Tightness' was defined as the ability of the structure to prevent relative movement between components and prevent the entry of harmful agents between adjoining interfaces [2]. The tightness also reflects 1 Chapter 1. INTRODUCTION ^ 2 the loading distribution ability of the structure. In turn, this ability is an indication of the bridge load carrying capacity. The bridge is a structure subjected to vehicle wheel loading, which is relatively concentrated. Slips between components of the deck (cover) would occur and form gaps. Water and incompressible materials (such as stones) would enter these gaps. The incompressible materials would force the gaps to remain open and make the deterioration of the structure worse, resulting in local deck failures after only a relative short service life. Increasing the tightness of the bridge means increasing its load distribution capacity of the bridge and, consequently, increasing its loading carrying capacity. The transversely post-tensioned laminating process, developed for this purpose, was first used in Ontario, in 1976, for the rehabilitation of a nailed-laminated bridge cover. Comparison of the load test results obtained before and after the transverse post-tensioning indicated the effectiveness of the procedure in improving the load distribution and decreasing the deflection in the bridge [2]. Since then, the transverse post-tensioning method has been accepted in new bridge design and construction. Timber bridges can be constructed in the four seasons of the year. This is especially important in Canada because of its winter weather. Their economic advantage has also been proven by cost statistics for short span bridges [3]. A method is needed to analyze the transversely post-tensioned timber bridge structure and obtain the relative interface movement under load. Such analysis could be used for the study of the reliability of the bridge for different performance criteria. The development of such an analysis is the main purpose of this research. Chapter 1. INTRODUCTION^ 3 1.2 Literature Review Previous research has developed different models for the two or three-dimensional analysis of stiffened plates. Cheung [4] developed the finite strip approach to analyze thin rectangular plates with two opposite edges simply supported . A trigonometric series function was used in the approximation of the deflection of the plate between the two simply supported edges, while a finite element representation was used in the cross-direction. For this case, the finite strip element method required less computer memory and increases the computational efficiency when compared to standard finite elements techniques. Foschi [5] used a combination of Fourier series and finite elements to analyze a type of timber floor system. This consists of a cover (deck) and the beams. The cover is fastened (nailed) to the timber beams to form a composite T-beam finite strip. The width of the T-beam flange is equal to the spacing between the beams. The flange (cover) is considered as an orthotropic thin plate in the analysis. FAP (Floor Analysis Program) was developed based on this model. The program can be used to consider gaps in the cover perpendicular to the span. However, the consideration of gaps caused coupling between the different Fourier terms in the global equations, requiring more computer memory and longer solution times. Thompson, Goodman and Vanderbilt developed the computer program FEAFLO [6],[7], to analyze the same wood floor system, taking account of the composite behaviour between cover and beams. The method considered the floor as a system of crossing Tand rectangular beams, including the effects of slip between layers owing to fastener deformation. Taylor [8] used the computer program ORTHOP to analyze the transversely posttensioned Hebert Creek timber bridge, which assumed no transverse flexural stiffness in Chapter 1. INTRODUCTION^ 4 the structural system while the transverse load transfer was due to vertical shear only. The friction coefficient between the timber interior surfaces was assumed to be about 0.5. Onate and Suarez [9] used Mindlin's plate theory to establish an analytical model taking into account the effect of transverse shear deformation. In their model, the transverse section was discretized by one dimensional finite elements and the longitudinal behaviour of the structure was defined by Fourier series. The element was called the simple two noded strip element with one single integrating point. Harik and Salamoun [10] also adopted the strip element to analyze the stiffened orthotropic rectangular plates. They idealized the stiffened plate as a system of plated strips and beam segments rigidly connected to each other. All these analyses only performed a linear elastic study of the structural system. Actually, the system response is nonlinear and the non-linearities due to friction between components must be considered in a more general analysis. 1.3 Objectives of the Thesis In order to obtain a general theoretical analysis, the non-linear properties for the transversely post-tensioned timber bridge must be considered. On the basis of the program FAP by Foschi [5], a non- linear, transversely post-tensioned timber bridge analysis computer program, PTB(Post-tensioned Timber Bridge), has been developed. The structure is shown in Figure 1.1. Nonlinear springs, representing the nonlinear frictional behaviour between T-beams, are inserted between the T-beam elements. The analysis is limited to static loading. The thesis objective is the development of a method for the calculation of the relative displacement between two adjacent T-beams, under the vehicle wheel loading. The method and the PTB program can be used to obtain the relationship between maximum relative Chapter 1. INTRODUCTION ^ 5 Figure 1.1: Analytical Structural Model displacement and post-tensioning force, permitting the calculation of the force required to achieve a target maximum relative displacement. The Newton-Raphson method [14] and the incremental loading method [15] are used to solve the non-linear problem. Since the non-linear behaviour of the friction springs produces coupling between the Fourier terms, the Jacobi iteration method [5] has been used to solve the coupled global equations. As an example, this Thesis considers a post-tensioned bridge with T-beams made with Parallam, a composite wood product manufactured by MacMillan Bloedel Ltd. of Vancouver, B.C.. Chapter 2 NON-LINEAR FINITE ELEMENT ANALYSIS 2.1 Structural Model The bridge structural model consists of two types of structural elements, ie. the T-beam and the joint between beams shown in Figure 2.1. The joint element, as mentioned in Chapter 1, is assumed to represent the non-linear friction behaviour between the interfaces of the adjacent beams. It is assumed that all the T-beam elements remain linearly elastic, with small deformations. The joint elements have non-elastic material properties. For the T-beam element, its flange can be assumed fastened to the web to produce an assembly capable of composite action, behaving as a stiffened plate under loading [5]. The flange may also be assumed to be rigidly connected to the web. Orthotropic thin plate theory [11] has been used in the analysis, to consider the different elastic properties of the flange in the directions perpendicular and parallel to the beams. A semianalytical procedure [12] has been used here, assuming in the x-direction a Fourier series for the displacement function, and a finite element approximation along the y-direction. The T-beam and the joint elements are illustrated in Figure 2.2 . The displacements of the T-beam and the joint element are matched at their common nodes. 2.1.1 Middle surface displacement functions for the flange of the T-beam These displacements are represented as follows: 6 Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ X Joint Element T-beam Element flange S web H Figure 2.1: The Structural Model S Y 1 Joint Element T-beam Element Figure 2.2: The T-beam Model and The Joint Model 7 Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS 1) in the z-direction N nrx 8 ) (2.1) nrx L ) (2.2) E F3n(y)sin( n Lr x ) (2.3) W(X,y) =^Fin (y)sin( n=1 2) in the x-direction N u(x,y) = 3) in the y-direction EF 2n ( y )cos( n=1 N v(x,y) = n=1 where: N = the number of terms used in the Fourier series L = the span of the bridge These displacements satisfy simply supported conditions at x = 0 and x = L. The functions Fin (y), F2n(Y) and F3n (y) can be expressed in terms of polynomials in the non-dimensional variable e = 2y/s, where: s = the spacing of the T-beams. A 5t h order polynomial is used for Fin(Y), and 4 th order polynomials for F2n (y) and F3n (y), with degrees of freedom associated with nodes 1, 2 and 3 in Figure 2.2. Thus, if {8 n } is the vector of nodal displacements, F in (e), F2,() and F3n (e) and their corresponding derivatives can be written as: F1 n (S) = {M0(e)}T {s F2 n(0 = { 114-3 ()} T { bn} F3n(e) = { ill5 ()} T { Sn} with corresponding derivatives: Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ dFin(e ) 1 9 = {MAW {an} (2.7) d 2 F1n(e) = {M2(e)} T { 8n} ck 2 (2.8) ck dF2n(e) ck dF3n() c14. = { 1114(e)} T { 8n} (2.9) = { 1116(e)} T { 6n} (2.10) The vectors {M1 (01 to {M6 ()} are given in Appendix A. 2.1.2 Displacement functions for the centroid of the T-beam web These are expressed as follows: 1) in the z-direction "W( N E W4nsin ( n ) :"--- n=1 2) in the x-direction rx ) .t N U(x) = E u4n cos(n:x )^ (2.11) (2.12) n=1 3) in the y-direction N V(x) = E V4n sin(7)^ n=1 4) rotation about x-axis 8(x) = N Ee nrx 4n sin( -) (2.13) (2.14) n=1 The vector {8,j, associated with the n th term of the Fourier series, is shown in Equation 2.15. Each T-beam element has 19 degrees of freedom. Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 10 U in n ln • S yin Vln • S U2n V2n W 2n • S IQ = < (2.15) W4n U4n "Kin 9 4n • 3 W3n W 3n • 8 U3n 1 U 3n • S v 3n v an • s The prime (') denotes the derivative with respect to y. For example, wi n • s is the derivative of w(x,y)with respect to y at node 1 associated with the n th item of Fourier series. The rotation degree of freedom is multiplied by the beam spacing s to make all components in {6n} dimensionally consistent. Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 11 2.1.3 Degrees of freedom for the joint element Similar to the T-beam element, the degrees of freedom for the joint element are associated with the each term of the Fourier series. win W ln • 5 Uln V in = (2.16) W2n w2n •£ 2n V2n The joint element consists of discrete three-dimensional springs at both the top and the bottom of the flanges along the bridge span shown in Figure 2.3. The dots show the positions of the three-dimensional springs which connect the flanges of two adjacent T-beam elements. The three-dimensional spring can act in the x-, the y- and the zdirection. The constitutive relationships for the x-, y- and z-direction springs are illustrated in Figure 2.4, where: fi x = the friction coefficient of the T-beam flange in the x-direction; A z =the friction coefficient of the T-beam flange in the z-direction; = the normal stress in the T-beam flange caused by post-tensioning force; e = the spring separation along the x-direction; t = the width of spring's effective area; Ax = the elongation of the spring in the x-direction; Ay = the elongation of the spring in the y-direction; Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ Three-dimensional Spring Figure 2.3: The Distribution of the Discrete Springs Fx Fy ii,a,e t o^Ay .pay e t x-direction E; y-direction Figure 2.4: The Constitutive Relationship of the Springs 12 Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 13 Az = the elongation of the spring in the z-direction The elongation of the spring is represented by the relative displacement between the two adjacent T-beam flanges at the corresponding location of each spring. When the relative displacement of the spring in the x-direction or in the z-direction (see Equation 2.54 and Equation 2.56) exceeds the corresponding friction limit, the spring loses its stiffness in the corresponding direction, but maintains a constant force equal to the friction force. The spring in the y-direction is used as a contact spring, to permit flange separation but not overlapping. When the relative displacement is negative or zero, ie. the two flanges tend to overlap each other, the stiffness of the y-direction spring is set very high, to enforce no overlapping between them. Thus, when the relative displacement is positive, which means that there is a separation between two adjacent T-beam flanges at the spring position, the spring is inactive or has no stiffness and no force in it. When the spring in y-direction loses its stiffness, this spring will not make contribution to the stiffness matrix. The effects of the springs in the x- and the z-direction depend on the `contact condition', ie, if there is a separation, the springs in the x- and in the z-direction are not active. 2.2 Strain Energy of the Structural System and Virtual Work The energy variational approach and the principle of virtual work are to establish the global equations for the structure. In this section, the strain energy of each component and the virtual work done by the conservative and the non-conservative forces are presented in terms of nodal displacement degree of freedom. Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 14 2.2.1 Strain energy in the flange of the T-beam element The strain energy per unit area of flange can be obtained by applying small deflection orthotropic plate theory. ^ Ufu Kx 8 2 w 2 Ky 8 2 w 2^a2w 82 w 2 ( axe + 2 ( 0y 2^Kv axe )( ay e 02w D x au 2 D y av 2 ( )2 + ( ) + ( ) 2 0x 2 Oy 0x0y 07) 12 ( 0u0v 1 D G au +2K0 ax ay ) + 2 `. ay Ox (2.17) For the plate with the thickness d, the stiffnesses Kx ,Ky ,K„,KG, D.,D y ,Dv and DG are given as follows: 1) the flexural stiffness in the x-direction Kx = Ex d 3 12(1 — 2) the flexural stiffness in the y-direction Ky = En d3 12(1 — vxy vyz ) and Kv = vzy Kx 3) the torsional stiffness KG Gd3 - 12 4) the axial stiffness in the x-direction Dx = Ex d (1 — vxy vyx ) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 15 5) the axial stiffness in the y-direction Eyd Dy = ^ (1 — vxy vyx ) and D, = vxy D z 6) the shear stiffness in plane DG G•d in which, d =the thickness of the flange; Ex = the elastic modulus of the flange material in the x-direction; Ey = the elastic modulus of the flange material in the y-direction; vxy = the Poisson's ratio, strain in the x-direction when stress is applied in y-direction; vyx ---- the Poisson's ratio, strain in the y-direction when stress is applied in the x-direction; G = the shear modulus in the x-y plane. The total strain energy Uf in the element flange can be obtained by integrating Uft, over the whole area of the flange. L 2 Uf^ F Ufudxdy 0 (2.18) Substituting Equation 2.17 and derivatives of the w(x, y),u(x, y) and v(x,y) (shown in Appendix A) into Equation 2.18, the total strain energy in the flange of one T-beam element is obtained in terms of the vectors {5,J. The integration makes use of the orthogonality of the trigonometric functions in the interval (0,L): Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 16 oL j cos( nr m rx = { 2 m = n )cos( Lx L 0 Tit n ) L r sin( nr : )sin( ^ mia rx ) 2 m=n 0 n n y 2 = — , de = - dy, dy = -s-de s 2 Thus, Uf 8 = E Ufi^(2.19) i=1 where the Ufi components are given as follows: A^L Uf i^12. {1 Ufu ldX}ClY = -r 0 N ic n 4 7.4 .9 1. 1 {sn}T{m0(e)}{m0(0}T{sn}de E^ 8L 3^- 1 n=i (2.20) 1 Uf2 = L 2. 2 {1 :FC n=i 2 Ufu2dX}dy n} { (0}{ 1112 W} T { En}4 -1 { S T M2 S) 3 2/1 (2.21) Uf 3 I • L {fu 3dX}dy 2. =E n=1 - - --Kv( 77L 1 ) 2 ( .9 ) 121^{ 5.} T {M0(e)}{M2()} T {Sn}de -1 (2.22) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 17 Uf4^ Za {^U f u 4dX}cly J0 N^nr s L t1 2KG ( )2(- )— --1{8n}T{M1()1{M1(e)}T{an}ck 2 2 L (2.23) Uf5^Ufu5dX}ely ET n s L fi ( L'77r-)2() 2 J- 1 {871}T{M3V)}{M3(°}T{Sn}4 (2.24) Uf6^ 1{1 L 1 n=i 0 Ufu 6dX}cly /1 fri} T {M6(0}{M6(0} T { 8.}4 2 s 2 -1 (2.25) Uf7^a{ 0 Uf u7 dx}dy 2 L Dv(17r-)— n=1^L 2 { 8n} T {M3()}{M6(E)} T { 8.}4 (2.26) Uf^ 8 1=1 {0iL Ufusdx }dy 2 DG SL 2 4 . 1 1 {5n } T [XX]{8„}ck (2.27) ^ Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 18 Where: [xx ] = ( -2.9- ) 2 {m4 ()1{ 1114()} T + ( -77- ) 2 {m5(0}{ m5(4)} T +(7-)(){{1115( 0}{ 1114(e)} T + {m4(e)}{m5(0} T } ^ (2.28) 2.2.2 Strain energy in the web This is given in terms of the displacements U, V and W and the rotation O. Thus, Uj Ely^d2W ,^EIz^, 2 d 2 V 2^EA fL,dU 2^GIt 1,d8 , 2 = dx+^dx+^dx 2.29 ) dx+ ) (^) 2 Jo dx 2 ^2 Jo dx 2)^2 o dx ^2 o dx Where: Iy = (BH 3 /12)^the moment of inertia about the y-axis; ^/z = (H13 3 /12)^the moment of inertia about the z-axis; OHB 3 the torsional moment of inertia; A = BH, the T-beam web cross sectional area; E = the elastic modulus of the web; G = the shear modulus of the web. = the torsional coefficient dependent on the ratio (H/B) [13]. It can be assumed that the ratio E/G is fixed. For most wood products, this is a large number, eg. E/G = 17.0. Substituting the derivatives of the W(x), U(x) V(x) and 8(x) (see Appendix A ) into Equation 2.29 the strain energy of the web can be written as: ^ Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 19 4 ^ uj^ (2.30) i=1 in which EIy r( d2W ) 2 dx 2 Jo dx 2 Ejl ^EI 2( m7r^inrx 12sin ,s.n,mrx 2 yy^fL^ )d x L L ^" ( L ) ( L )a n=1 m=1 EIy v (winl or vi( L (2.31) 2 nn=1) L ) ( 2 EIz I, eV ( Uj2 2 o dx 2 vnv,4m ( nr )2 mr \ 2sin( nrx )sin( mrx )dx ^m=1Jo ^L^L I^L^L EIy ‘.1\r\ ■‘1‘r 2y EIz 2 )2dx I L ( L N E(V4n ) 2 ( 17^2 ) 4 (2 ) n=1 (2.32) EA iL ( dU )2dx 2 Jo dx EA \1v , 2.",r U j3 (nr y mr \sin( nrx , s . n , mrx ^ )d x ‘---d o^L ( L )^k L ) I ( L 2 n=1 L-'dm=1 IL u4nu4m EA N 2A n=1 ^) (2.33) = (2.34) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 20 2.2.3 Strain energy in the web-flange connector Three kinds of deformations have been considered in the connectors between the flange and the beam 1) The deformation of the connector parallel to the beam: ru2_ d ( aw 2) 1^H dW1 emu= [^2 Ox^[^2 dx (2.35) 2) The deformation of the connector perpendicular to the beam: Ay^[V2 d ( aw^2 )]^[V + I 19] 2 ay^2 (2.36) 3) The rotation deformation of the connector about the x-axis lv a2 = ( ay ) 9 (2.37) Where: d = the thickness of the flange H = the height of the web u 2 = the displacement in the flange at point 2 in the x-direction v 2 = the displacement in the flange at point 2 in the y-direction = the displacement in the flange at point 2 in the z-direction The flange is assumed connected by uniformly spaced discrete connectors along the longitudinal centerline of the beam. Thus, the total strain energy in the connectors in each T-beam is: UN = Au) ? + 4 1 (64 +^ ?. ( 2.38) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 21 Where: NA = the number of the total connectors in one T-beam; (Au) i = the deformation in the x-direction for the i th connectors; (Av) i = the deformation in the y-direction for the i th connectors; (00) i = the rotation deformation for the i th connector about the x-axis; k z = the stiffness of a single connector corresponding to the x-direction; ic y = the stiffness of a single connector corresponding to the y-direction; k e = the rotation stiffness of a single connector. Alternatively, an equivalent continuous connector model can be used to calculate the strain energy of the connectors per T-beam. UN = IL Au\2 ky 12ek o v , 2e^2e (A0)2]dx (2.39) where e = the connector spacing along the beam. Substituting Equation 2.1 - Equation 2.3 , Equation 2.11 - Equation 2.14 and the corresponding derivatives into Equation 2.35 - Equation 2.37, N rtrx d ) Fin (y = 0) cos( ^ ) L^2^L^L n=1 nrx )— Tir nir nrx cos( U4n cos ( 2 vv 4n L^ L )1 Au^E [F2n ( y = 0)cos( - nrx - (2.40) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 22 N d dFin (y = 0)^nrx Av = E [F3n (y = 0) I74n^e4n]sin( ^ ) 2^dy 2^L n=1 N^ d ,^ nrx E [v2. — 2 w 2n — — 94n]sin( L ) 2 n=1 —^ N AO =^[w2n n=1 nirx 04n]sin(— T (2.41) (2.42) Finally, introducing these three Equations into Equation 2.39, UN N L 1 =^E^[u2n —^— W4n(H d)] 2 n=1 2 2e d ,^H n +ky [V2n — 11)2 V4n 2 — n + k9 [w'2 n — 04,] 2 } n ^ (2.43) In order to express the strain energy of the connector in terms of the nodal degrees of freedom, we introduce a set of (19x1) vectors { e i } such that the i-th entry row is unit, but all other rows are zero. Substituting Equation 2.15 and {e 7 }T , {e 8 }T { e l3 }T into the Equation 2.43, the strain energy of the connector can be written as follows: UN N =E L 1 ke {k{X} + k y {Y} + —{0}} 2e n=1 2^ (2.44) Where: {X} = {8n } T [{e 7 } — {e n } — {e io }(H + d)] • [{e 7 } — {e n } — {e io }(H + d)] T {8n }^(2.45) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 23 {Y} = {an}T [{e8} [{e8} fegl§ {69}279 { e12 } — 13 _ —1 1— {613}' 2sd T 51 {an} {O} = { 5n} T Re9} — {e13}1. ^ Re91 — {en}i{sn} (2.46) (2.47) 2.2.4 Virtual work done by non-conservative forces in the joint element The joint element consists of 2N 8 'three dimensional' springs. The deformations of the i th `three-dimensional spring' in the x-, the y- and the z-direction are expressed, respectively, as follows. When a spring is at the top of the flange its z-coordinate is —0.5d, and its z-coordinate is 0.5d when it is at the bottom of the flange. Along the x-direction, u(x i , zi )^u(x i ) — z i 5w C aw (2.48) . with the spring elongation then being, Au(x i ,z i )^u2(xi,zi) — u 2 (x i ) zi ( 5w1 ) Ox xi zi (2.49) Similarly, along the y-direction, v(x i , z i )^v(x i ) — z i Caw Ow (2.50) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 24 with elongation Av(xi,zi)^v2(xi,zi)— vi(xi,zi) awl ^(aw2 z, Vi (Xi ) + V2(Xi) +^ ( — ay)^ay) (2.51) xi For the z-direction, w (x i , z i )^w(x,)^ (2.52) and the spring elongation is: Aw(x i ,^=-- w 2 (x i ,z i ) — =^w2(xi) ^ - (2.53) Substituting the displacement functions of the T-beam flange and the corresponding derivatives in Au(x i ,z i ), Av(x i ,z i ) and Aw(x i ,z i ) above, N^ nir^fir^ nrxi ^w1n ^Au(x i , z i ) E[—u1n+u2r, +z i^w in ]cos ( ) (2.54) n=1 Av(x i ,^= E[—vin n=1 nirxi V2n ZiWln ZiW2n]Sin( ^ ) N Aw(x i , zi) = E[—w1n + w2n]sin( nrxi ) (2.55) (2.56) n=1 which can be written in terms of the nodal degree of freedom vector. (2.57) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ N^N = n=1^n=1 EfAnl = EAdfs'l ^ 25 (2.58) Where: {{Au n } {An} = {Av n } 1 (2. 59) {Awn } z i (-yr )cos(Ti )^0^—sin(7-71) 0^7-9isin(n)^0 --cos ( mr: i )^0^0 [Br,]T = r 0^—sin( m i )^0 - Zi (n7 )C08( -7"--ri)^0^sin(7-i-) (2.60) 0^—isin(Ti)^0 cos ( "---1-7i )^0^0 1 0^sin(12T)^0 The forces in one 'three-dimensional' spring will be: N {F} = [D]{A} = [D]^[Bn ]{8,J,} n=1 (2.61) Where: [D] = EI 0 0 0 E; 0 0 0 E: = the stiffness modulus of the spring in the x-direction Ey = the stiffness modulus of the spring in the y-direction (2.62 ) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 26 E: = the stiffness modulus of the spring in the z-direction Assuming that the structural system is subjected to a virtual displacement, corresponding to a virtual change in the nodal displacement vector for the joint element {,51-*, and the arbitrary deformations of the spring can be obtained: N {A} - = E113.1{8,7,}* ^ (2.63) n=1 The internal virtual work in one 'three-dimensional' spring is, therefore, T47 = [{A}1 T {F} N = E({87.T)T[Bn]T{F} n=1 N N = E E (feLY) T [Bn] T [D][B.]{ 8;17,} (2.64) n=1772=1 The internal virtual work in one joint element is then the sum = 2N. ^ i=1 (2.65) where 2N8 = the number of the total three dimensional spring in one joint element. 2.3 Global Equation The global system of equations for the structural system is [G]{U} {P}^ in which: [G] = the global stiffness matrix (2.66) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 27 {U} = the global nodal displacement vector {P} = the global consistent load vector The global stiffness matrix, the global nodal displacement vector and the global consistent load vector are assembled by the total contributions from each individual element using the conventional finite element method. 2.3.1 Stiffness contribution from one T-beam element The stiffness contribution of one T-beam element can be derived by minimizing the total strain energy in the T-beam element. The total strain energy of the T-beam is the summation of the strain energy of the flange, the strain energy of the web and the strain energy of the connectors. U = Uf +Uj + UN 8^4^3 =E Ufi E u i + E uNi^ j i=1^i=1^i=1 (2.67) Take the first variation of the strain energy of the flange with respect to the {S n }, we get the stiffness contribution from the T-beam. 8 (Un) )71 + 8 (Uj )7/ + 8 (UN)n 8^4^3 a(Uf = = E(sum + E(suji )„ + E(SyNon n a(ifi)n^ i=i (2.68) K. n 4 r 4 s 1 4L3 [/if-M-0(0}{M0(0}Tda5.1 (2.69) 8 (Uf2). = - P^ ) L[f i{M2(e)}{M2(01T41{6.} . 3 l Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 28 (2.70) a(Uf3 )„ = —Kv ( L ) 2( 2 (2.71) 4U/4) n = 2KG(7)2(2L)[f1{M1()}{M1(e)}TAR8n} (2.72) ; L li { M3 ()}{ M3 (e)}7 41{ 8n} a ( Uf 5 )n = 221- (7) 2 [f (2.73) 8(Uf6) n = 2!-L{111{114-6()}{M6()}Tde]{an} (2.74) 8(Uf7) 71 = D, ( 74 [ 111{ m3 (4)}{ m6 (4)} Td 4ll sri l (2.75) 8(Ufs). = DG 2 ( )L[ { li ( s2 )2 m-4 ( ) 1{ 14 (0}T ig + 1 ( 7 )2{m ()}{111 ( )}Tde 1 5 5 nr 2^ ri 7r 2 + .1_ 1 (i);{M5(0}{M4(e)/ T 4 L i ( 1217);{M4(0}{M5(0} T ckif 8n} (2.76) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 4 ^nr E ^ W 4n( + ^I^Trir ^z V4 „(— ) 4 L )4 L + E L^ 2 L EA nr nr GIt u4„( )2L o4n(— ) 2 L 2^L^2^L 3 L E^— flz[X][.X] T 2e 1=1 29 ke ky[Y][Y]T — [0][0]T}{8n} (2.77) (2.78) [X] = [fe 7 1— {e n } — {e io }(H d)]^ (2.79) r^r^i H i len/ le131 . 0 [Y] [{e8} [ 0 ] = [ { e 9 } — { e 1 3 }]^ (2.80) (2.81) 2.3.2 Stiffness contribution from one joint element Because there are non-conservative forces in the joint element, the corresponding stiffness contribution is obtained by using the principle of virtual work. As shown in the previous section, the internal virtual work done by the spring forces in the joint element is: 2N. N N E E E ({8T)T[Bn]T[D][Bm] le„} (Wirt =i=1 n=1 m=1 8 (2.82) Since {8;1 }* is an arbitrary nodal displacement the stiffness contribution from the joint element corresponding to the n th and in th terms in Fourier series is: Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 30 X P(x,Y) PT PT y2 y1 Y Figure 2.5: The Vechicle and Post-tensioning Load 2N, BK[n,m]^[B„]T [D][B m ]^ (2.83) BK[n,m] is a 8 x 8 symmetric matrix, producing coupling between the Fourier terms. This coupling is reduced ( i.e. BK[n,m] approaches [0] if n m) when N, is large and all springs are elastic, since such limit approaches the orthogonality situation among the trigonometric shape functions. 2.3.3 Consistent load vectors Two kinds of loading are shown in Figure2.5. Vehicle wheel load Under a virtual displacement w (x,y), the potential of the vehicle wheel load p(x,y) is - Ui = 1 °2 I Y2 fulx,y)} T {p(x,y)}dxcly 01^yi ^ ^ Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ j(x2r2 j ^ r ^{ {Mo ()} 7' {S„}*} T {p(x,y)}sin(n )dxdy = Lx xi Y1 n=1 31 (2.84) From the definition of the consistent load we have 1V x2 Y2 >{{Sn}*}T{Rn} = Yi n=1^ {iv*(x,Y)}T{P(x,Y)}dxdY N - 1M°1 I Y2 IE{{84*}T{Mo(E)}{P(x,Y)}sin(nr 2 X )dxdy (2.85) Then: {R.}= f r{m0(0}{pcx,yilsin(n ^ )dxdy x 1^yi where: {R n } = consistent load corresponding to the n th (2.86) term in the Fourier series. If p(x, y) = P, a constant uniformly distributed load over the patch, {R.} = P^ 2 iziz2 1: {Mo(e)}sin( : 2 n x )dxd P s L^nirx cos( n 7 2 )} 166 {Mo ()lck 2 nr {cos( Li) (2.87) Where: 2y i^2y2 el —;^ = s^s Post-tensioning load The individual concentrated post-tensioning load is assumed to be applied, respectively, at node 1 and node 3 of the first and last T-beams. Under an arbitrary displacement e(x,y), the potential of the post-tensioning load PT(x i ) is NT UPT = Ely ( ,O} T {PT(xi)} i=1 (2.88) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 32 From the definition of the consistent load we have : E{{8 ,}1 T {R,} = Et{8 7 NT n }l T {M5()1 nrx (2.89) E PT (x )sin( ^)^ i n=1^n=1^i=1 NT^nrx . sin(^')PT(x L i )^ {R„,} 5 = E i=1 NT {R-.}18 — sin( L )PT(x ) E ^nrxi (2.90) (2.91) i i=1 where: PT(x i ) = the value of the it h pair of post - tensioning forces; NT = the total number of the post-tensioning forces; x i = the location of the i-th post-tensioning force. 2.3.4 System of equations corresponding to one T-beam (Element Matrix) The system of equations corresponding to one T-beam element has the following form: (2.92) [AKe]{51 = {Re} where: [AK (1,1)]^[0]^. [AK'] = [0]^[AK(2, 2)]^... [0]^[0] [0] [0] [AK(N, N (2.93) )] Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ = 33 (2.94) { 6- eN} {R7} {R e } (2.95) {RN} Each submatrix AK [n, n] is a 19 x 19 symmetric square matrix. N is the number of items in the Fourier series. {Sn} and {R„e } ( n = 1, 2 ... N ) are the 19 x 1 vectors. 2.4 Global Equations and Solution Method The characteristics of the global equations and the solution strategy will be discussed with an example, using only two terms in the Fourier series and three T-beams ( NJT=3 ). 2.4.1^Characteristics of the global equation The form of the global equations is shown in Figure 2.6 and Equation 2.96: [GK(1,1)] PK(2,1)] [GK(2, 1)] PK(2,2)]_ I {U1} {U2} = { P1 } (2.96) j^{P2} AK 1 (1, 1) is the stiffness contribution from the T-beam element 1 associated with the first term of the Fourier series. As mentioned in the previous section, it is a 19 x 19 symmetric matrix, as are each of the other AK i (n,n) ( n = 1, 2 ) matrices. BK i (n,m) is the stiffness contribution from the joint element i associated with the n-th and m-th terms in Fourier series. It is a 8 x 8 symmetric matrix. Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ n=1 ^ 34 n=2 n=1 n=2 P2 Figure 2.6: The Global Equation System GK(n, n)(n = 1, 2) is the (NJT * 19) 2 symmetric banded matrix with a band width of 19 made up of contributions from the AK(n, n) and BK(n, n); GK(1, 2) and GK(2, 1) only have the contributions from the joint elements. {Un }and {Pn } are the (NJT * 19 ) x 1 displacement vector and consistent load vector associated with the n th term of the Fourier series 2.4.2 The incremental procedure The incremental method [15) is used to solve the non-linear global equations. The initial values play a very important role in the non-linear analysis iteration method. The Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 35 incremental method solves the global equations by replacing the solution to [G(U)]{U} = {P} by successive solutions of [G(U) i ]{AU2} = {APi } The final solution is {U} = E 3i=1 {AU2 } Where: {AUi } = the incremental displacement of each step. {APi } = the incremental load The initial value used in order to get {Ui } is the solution { U j _ 1 } in the previous step. With in each step, the global stiffness matrix must be updated according to the current displacement. The final solution is the summation of all displacement increments. The process is illustrated in Figure 2.7. 2.4.3 Newton-Raphson method From the Figure 2.7 we could note that the more the number of steps, the more accurate the solution is. In order to decrease the error within each incremental load step, the Newton-Raphson method has been used in the PTB program. This method is illustrated in the Figure 2.8. The final solution in each incremental step {AUi } is obtained when the Newton-Raphson method converges. Thus, {AUi} = E{Au k k=1 } Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ A P2 API A U2 U1 Figure 2.7: The Incremental Method AP1 A A A U 1^ 4, pu2 Figure 2.8: Newton-Raphson Method 36 Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 37 (2.97) [G(U)]t{Auk} = {AQh}^ I I { Auk} II 2 < cid 1 1{ uk _1 } 112 where: [G(U)] t = the global tangent stiffness matrix which must be updated to the current displacements; {AC} = imbalance load vector; Eid {u k = a tolerance; - 1 } = Etil{Dui}; Il{uk}112 = {Ei= 1 uZ i } 1 / 2 the Euclidean norm of the vector {u k }. The derivations of the tangent stiffness matrix [G(U)] t and imbalance load vector {A4 k } are expressed below. At the element level [K(u)]{u}^{p}^ (2.98) then {R} = {p} — [K(u)]{u} ^ (2.99) To guess {ui_ i }, we get a residual force vector {R 1-1 } = {p} — [K(ui- 1 )]{ui- 1 } 4 0 ^ (2.100) We look for {u k } and make {Rk } equal zero {u k } = { uk_i} + {Auk}^ (2.101) Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ {Rk } 38 = {R(u k _ i + Au k)} = {R(u k )}^0 ^ (2.102) Expand {Rk } in Taylor series in the neighbourhood of {u k } and only take the first two terms. „^19{R} , {R(u k _ i + Auk)} = {R(uk-1)1 }{uk-i} • {Auk} = {0} (9{ u (NM- 9{ 14 1ft,^{Auk} k-1 ( ( 2. 1 0 3) {R(uk-i)} {p} — [K(uk_i)]{uk-il (2.104) where: (9{R} a{u} Ifuk^= ig(U)lt [g(u)] t is the element tangent stiffness matrix corresponding t o the { uk _ 1 } and {u k _ 1 } is the current displacement. {p} - [K(uk_i)]{uk_il = {q(u)} {q(u)} is the imbalance load vector corresponding to displacement {u k _ 1 } at the element level. At the global level Assemble all element tangent stiffness matrices [g(u)] t and imbalance load vectors {q(u)}by conventional finite element method, we get the global tangent stiffness matrix [G(U)]t and the global imbalance load vector {.64k}. Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS ^ 39 2.4.4 Jacobi iteration method Within each iteration of the Newton Raphson method, the Jacobi iteration method is used to solve the global system of Equation 2.97. This avoids storing the entire global matrix, requiring less computer memory . If the diagonal matrices are dominant, the rate of the convergence of this method is very high. The procedure can be written as follows: N { Al k};n = [CK(m,m)] t {{AQic}ni — E[GK(m,n)t{Auk}in 1} (n m) (2.105) n.1 with starting vector {Au k };„ = [G K (m, m)]t {AQ k}n (2.106) The iterative procedure is stopped when the following convergence condition is satisfied: Il{Auk}z„, — {Auie}n 1 112 < 611{Auk}in,71112^(2.107) m = 1,2,...,N where: E is the error tolerance (which could be .001 for example ). Theory presented in this Chapter has been incorporated into computer program PTB given in Appendix B. Chapter 3 FRICTION PARAMETERS FOR THE SPRING MODEL 3.1 The Friction Test In order to determine the friction characteristics of Parallam, an experiment was conducted at MacMillan Bloedel Limited Research Centre. Four kinds of surface textures were studied: • wet and preservative treated, • dry and preservative treated, • wet and untreated, • dry and untreated. The test setup is shown in Figure 3.1. The dimensions of the top specimen were 63.5mm x 63.5mm x 38mm and the size of the bottom specimen was 88.9mm x 88.9 x 38mm Both were cut from rough-sawn planks. It should be pointed out that the surface of the preservative treated specimens was much rougher due to the chemical treatment process. The specimens classified as `wet' were soaked in water for 24 hours before test. All test specimens were measured before being soaked in water and weighted just after being taken out of water, and wiped clean of surface water. For each kind of surface texture, both the friction coefficient for the sliding direction perpendicular and parallel to Parallam fibres were determined, as shown in Figure 3.2. 40 Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL ^41 weight sliding clamp ^ "`- top specimen bottom specimen Figure 3.1: The Friction Test Setup direction sliding Parallel direction of sliding Perpendicular Figure 3.2: Two Kinds of Friction Test ^ Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL ^42 Specimen Type Table 3.1: Data of Friction Test Samples Average Weight (gr) Number of Specimens Dry, Treated, Parallel Dry, Treated, Perpendicular Wet, Treated, Parallel Wet, Treated, Perpendicular 110.99 111.16 141.15 140.08 20 19 17 14 Dry, Untreated, Parallel Dry, Untreated, Perpendicular Wet, Untreated, Parallel Wet, Untreated, Perpendicular 92.36 92.75 135.18 137.47 15 12 19 12 Weight of the Clamp Applied Weight 109.91 5000.0 Table 3.1 gives the data of the samples used in each type of tests. In the perpendicular direction test, the sliding direction was perpendicular to the Parallam fibres. In the parallel direction test, the sliding direction was parallel to the fibres. The force FH required to produce sliding between the surfaces of the specimens was proportional to the force applied normal to the plane of motion. The ratio is defined as the coefficient of friction. FH =^ Fv (3.1) In all tests, an Instron 4210 testing machine was used to provide a 0.4in/min. constant sliding velocity. The horizontal force FH vs. sliding movement relationship through the test was monitored and recorded, and typical curves are shown in Figure 3.3 to 3.5. From these Figures, it can be seen that the curve of the horizontal force vs. the movement, after sliding, is not smooth, reflecting the stop-and-go induced by surface roughness. The results and statistical data are presented in Table 3.2. Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL Table 3.2: Statistical Results of the Friction Test Classified Specimens a of it Ct, of it E ( N/mm ) tt Dry, Treated, Parallel .655 .115 .076 16.902 Dry, Treated, Perpendicular .831 .080 .096 17.378 Wet, Treated, Parallel .046 .833 .038 12.618 Wet, Treated, Perpendicular .059 12.353 .878 .068 Dry, Untreated, Parallel .031 .083 11.183 .373 Dry, Untreated, Perpendicular .419 .041 .098 10.891 Wet, Untreated, Parallel .804 .052 .065 14.900 Wet, Untreated, Perpendicular .811 .048 .059 12.339 ^ a- of E 2.895 1.922 3.197 3.863 1.751 1.992 2.523 1.943 43 CC, of E 0.171 .111 .254 .313 .157 .183 .170 .158 = friction coefficient (mean); E = stiffness ( initial slope of the FH vs. displacement curve)(mean); or = standard deviation; C, = coefficient of variation. 3.2 Spring Model As previously mentioned, a 'three-dimensional spring' model is assumed in the non-linear analysis and its constitutive relations are shown in Figure 2.4 of Chapter 2. The friction coefficients A x and A z are directly adopted from the friction testing. The stiffnesses E; and E: ,on the other hand, are derived from the friction testing results using the 'same slip' assumption. We can idealize the friction test results in Figure 3.6. In the spring model let the elastic displacement limits in the x and the z-direction be, respectively, AL, and Afin, . Thus, fixcryet Ali en = ^ E; (3.2) Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL ^44 Figure 3.3: Friction Test Recording Curves(1) Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL ^45 F(1, ) Dry Treated Penendicular 2 DEtn). ) 0 5 1.0 F(lb) 4 Dry'Untreate:: Parallel 0.5 Figure 3.4: Friction Test Recording Curves(2) Chapter 3. FRICTION PARA1VIETERS FOR THE SPRING MODEL ^46 ^Ff Its F 10^ fi 1- met Treated _Perpendicular' Wet Treated Parallel, D(in ^ 0.5^1.0 0.5 Figure 3.5: Friction Test Recording Curves(3) 1.0 Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL ^47 Figure 3.6: Idealized Friction Test Curves = A zEo- su et (3.3) In the parallel direction of the friction test: Aix = FH El (3.4) In the perpendicular direction of the friction test: z = FH (3.5) z The 'same slip' assumption implies that = Al. Ati m = Aix^ from which, (3.6) (3.7) Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL ^48 0= Figure 3.7: Same Slip Assumption El = E: = Ex t crYet FT' EtzcrYet Fir (3.8)^- (3.9) E; and E.; are the stiffnesses used in the 'three dimensional' spring; Ext and E are the stiffnesses obtained from the friction test. The relations between the spring stiffness and the stress ay is clearly shown in Figure 3.7. We can note that E; and E; are proportional to the normal stress cry , however the maximum friction force in the spring cannot exceed the shear strength of the material rs et, where Ta is the shear strength of the flange and et is effective area of the spring. Chapter 4 ANALYTICAL RESULTS 4.1 Beam Test In order to test the spring model described in Chapter 2, a simple beam structural problem, shown in Figure 4.1, is now considered. The structure consists of two builtin beams connected by four springs. At the top, there are one y-direction and one z-direction spring. The other two are assumed in the y and z-directions at the bottom of the beams. The left beam is subjected to one vertical concentrated force at midspan. The properties and the dimensions of this structure are shown in the Table 4.1 The structure was analyzed by a specific program using the same spring model theory, which is used in the PTB program. When P > 0.0, this axial load produces tension in both beams. The results given in the Figure 4.2 show the displacements of node 2 and 3 vs. the tensioning force P. The status of the four springs as a function of the post-tensioning force P is given in the Table 4.2 In Table 4.2, F means the force in the spring is beyond the friction limit; E means that the deformation of the spring is within the elastic range, less than the friction limit Atim• From Figure 4.2 and Table4.2, we can note that increasing the post-tensioning force P influenced the status of the springs, with an associated change in load distribution within the structure. Finally, when the P reached 100, 000(N), the two y-direction springs were 49 50 Chapter 4. ANALYTICAL RESULTS^ 4 L L B Figure 4.1: Two Cantilever Beam Structure Table 4.1: Properties of the Cantilever Beam BEAM 200 (mm) width of the beam ( B ) 400 (mm) depth of the beam ( H )^. 4000 (mm) span of the beam ( L ) 14000 (MPa) elastic modulus of the beam (Eb) elastic modulus of the y-dir. spring (E y ) 1.0E10 (MPa) elastic modulus of the z-dir. spring (Er ) 17.4 (MPa) .65 friction coefficient (/./ z ) LOADING vertical concentrated loading ( V ) -20,000 (N) Table 4.2: Status of the Springs in the Cantilever Beams STATUS OF THE SPRINGS Post Force (N) y-dir. at top z-dir. at top y-dir. at bottom z-dir. at bottom open not active contact F 0.0 open not active F contact 10 F open not active contact 100 open not active contact F 1,000 not active E open contact 10,000 E contact contact E 100,000 E contact E contact 1000,000 Chapter 4. ANALYTICAL RESULTS ^ VERTICAL DISPL. vs. POST TENSIONING FORCE " 4.00 51 . NODE 2 --B-NODE 3 2.00 CD 0.00 a) 8-2.00 co -4.00 - aCD -6.00 - 'C^o^ CD -8 10.8E+00 1.0E+02^1.0E+04^1.0E+06^1.0E+08^1.0E+10 Post-Tensioning Force ( N ) Figure 4.2: Cantilever Beam Structure Displacements L/2 3 L Figure 4.3: Built-in Solid Beams Chapter 4. ANALYTICAL RESULTS^ 52 Joint 1^Joint 2 1 Z y 2 3 Figure 4.4: Three T-beam Structure active ( full contact between the beams) and the deformations of the two z-direction springs ( in terms of relative displacement at the interface of the beams ) were all within the elastic deformation limit. When the P reached 100, 000(N) the displacement of node 2 was decreased by several times in comparison with that for P = 0.0(N). The displacement of node 2 must then be nearly equal to the displacement of node 2 in Figure 4.3. The structure in Figure 4.3 is a solid beam and all its properties are the same as those for the case of Figure 4.2. This was independently verified with a structural analysis of the beam in Figure 4.3. 4.2 T-beam Analysis The configuration and the other properties of the structure are presented in Figure 4.4 and Table 4.3. The bridge is assumed built with Parallam T-beam (dry, treated conditions). Since a sine series in the x-direction is used to approximate the displacement w(x,y) Chapter 4. ANALYTICAL RESULTS^ Table 4.3: Properties of T-beam Structure Parameter Notation Value T-beam span L 12000 ( mm ) T-beam spacing 1000 ( mm ) s web height II 500 ( mm ) web width B 100 ( mm ) the ratio of elastic modulus REG 17.0 to shear modulus of the web flange height 100( mm ) d friction coefficient in x-direction 0.50 F. friction coefficient in z-direction 0.65 fiz three dim. spring number in one joint element 42 2N, spring elastic modulus in the x-direction Ex 16.9 (N/mm) spring elastic modulus in the y-direction EY 1.0E10 (N/mm) spring elastic modulus in the z-direction Etz 17.4 (N/mm) web elastic modulus in the x-direction 14000.0 (MPa) EL eb web elastic modulus in the y-direction EYweb 1400.0 (MPa) flange elastic modulus 14000.0 (MPa) E flange post-tensioning force spacing 1200 (mm) 53 Chapter 4. ANALYTICAL RESULTS^ 54 PT L PT Figure 4.5: Three T-beam Structure with Post-Tensioning Force Only (deflection) of the structure, simply supported boundary conditions are satisfied at x = 0 and x = L. Other edges are free. The analysis considered five loading cases and was done by running program PTB. Since all cases are symmetric about midspan in the x-direction. four Fourier terms with order 1,3,5 and 7 were chosen in the analysis. The objective was to study the effect of the post-tensioning force in the loading distribution within the structure. The program PTB can give either nonlinear or linear elastic results, the latter implying that all springs keep their initial stiffness throughout the loading. Studying the relative displacement in z-direction between the flanges of adjacent the T-beams constitutes the main objective of this analysis. 4.2.1^Analysis for post-tensioning force only The structure and loading condition are shown Figure 4.5 Three T-beams were subjected to ten pairs of post-tensioning forces along the span. Chapter 4. ANALYTICAL RESULTS ^ 55 Table 4.4: Analytical Results: Post-tensioning Force Only PTB PROGRAM SAME SLIP ( LINEAR) POST FORCE (N) 1.00 10,000 Maximum Displacement ( T-beam 1) ( mm ) -0.57952E-05 -0.57726E-01 Maximum BEND. STRESS ( T-beam 1) ( MPa ) -0.14474E-05 -0.14681E-01 Maximum Displacement (T-beam 2 ) ( mm ) -0.58064E-05 -0.58461E-01 Maximum BEND. STRESS ( T-beam 2 ) ( MPa ) -0.11604E-05 -0.11126E-01 PTB PROGRAM SAME SLIP (NON-LINEAR ) POST FORCE ( N) 1.00 10,000 Maximum Displacement ( T-beam 1) ( mm ) -0.57952E-05 -0.57726E-01 Maximum BEND. STRESS ( T-beam 1 ) ( MPa ) -0.14474E-05 -0.14681E-01 Maximum Displacement (T-beam 2 ) ( mm ) -0.58064E-05 -0.58461E-01 Maximum BEND. STRESS ( T-beam 2 ) ( MPa ) -0.11604E-05 -0.11126E-01 From Table 4.4, we can note that the results from the non-linear analysis and those from the linear analysis are identical. This means all springs were in the elastic range at the different levels of the post-tensioning forces. Figure 4.6 gives the relative deformations of the joint element 1 in z-direction along the span of the T-beam. The negative deformation of the joint element indicates that the flange of T-beam 2 moved down less than the flange of T-beam 1 did. 4.2.2 Analysis for the one patch of vertical loading The load case is shown in Figure 4.7. This is a symmetric problem. Along the 12m long span 10 pairs of post-tensioning force were acting. Each pair of the post-tensioning force was equal to PT. The structure carried one patch of vertical loading equal to 500 mm x 300 mm x 0.2 (N/mm 2 ) = 30,000 (N) at midspan of T-beam 2. For the symmetric Chapter 4. ANALYTICAL RESULTS ^ 56 5E-08 0 -5E-08 -1E-07 E • -1.5E-07 -2E-07 0 as -2.5E-07 -3E-07 .S^-3.5E-07 0 E Uco 0. cn . 2,000^4,000^6,000^8,000^10,000^12,000 The x-coordinate along the span (mm) PT = 1.0(N) 0 77^-0.0002 Ti • 4.0004 -0.0006 -0.0008 -0.001 -0.0012 0^2,000^4,000^8,000^8,000 The x-coordinate along the span (mm) 10,000 12,000 PT = 10,000(N) Figure 4.6: Displacements of T-beam Structure under Post-Tensioning Force Only Chapter 4. ANALYTICAL RESULTS^ 57 Figure 4.7: T-beam Structure One Patch Loading Case problem , we only plot the stresses and the displacements of the T-beam 1 and T-beam 2. Figure 4.8 to Figure 4.10 show the effect of the post-tensioning force on, respectively, the maximum bending stresses in the webs and the maximum vertical displacements in beams. When PT = 0.0, the total vertical loading is carried by beam 2 only. The maximum relative displacement between beam 1 and beam 2 ( ie. the maximum deformation of the Joint element 1 in the z-direction ) was more than 20mm. With the increase in the post-tensioning forces, the link between beam 1 and 2 became tighter and, consequently, the load sharing behaviour of the structure improved dramatically. When PT approached 10,000 (N), the non-linear analysis results were identical to those from the linear elastic analysis. When PT reached 100,000 (N), the maximum relative displacement, in z-direction, between the flanges of beam 1 and beam 2, was negligible when compared to the 20mm for PT = 0.0. • Chapter 4. ANALYTICAL RESULTS ^ 58 co a. 110 T-beam 1 T-beam 2 E °3 6 (D co .5 4 o u) - ED 2 . (.7) ...O. • 0 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07 The Post-tensioning Force ( N ) Figure 4.8: Bending Stress of T-beam Structure One Patch Loading Case • 25 T-beam 1 20 T-beam 2 E X15 co co 10 5 5 A l 0 cp 76' 9- iE+00 LJ 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06 ^ The Post-tensioning Force ( N ) Figure 4.9: Displacements of T-beam Structure One Patch Loading Case 1E+07 Chapter 4. ANALYTICAL RESULTS^ 59 Joint Element 1 20 a) U SI 0. 15 2 .6 10 5 > fa' cl" a+oo cc * 1E+01^1E+02^1E+03^1E+04^1E+05 The Post-tensioning Force ( N ) 1 E+06 1E+07 Figure 4.10: Relative Displ. of T beam Structure One Patch Loading Case - In the example shown in Figure 4.9, upward (negative) deflections of the T-beams can be seen to be developing as the post-tensioning force becomes large. The reason for this behaviour is explained as follows: when the flange of the T-beam subjected to post-tensioning forces in the y-direction, the flange will expand in the x-direction due to Poisson's effect. Since the flange is connected to the T-beam web by the connectors with non-zero stiffnesses, shear forces develop to maintain compatibility. The direction of these forces as shown in Figure 4.11, is such that they cause an upward deflection which, shown in Figure 4.9, becomes more important for post-tensioning forces greater than 1E + 05N. Figure 4.12 is the maximum displacements of the T-beams caused only by one patch loading. Those displacements were measured from the upward deflection caused by posttensioning forces. See Figure 4.11. It is very obvious that the post-tensionging force can improve the loading sharing capacity of the whole structure. Chapter 4. ANALYTICAL RESULTS^ 60 * PT Flange Web 4^4^4 4^4^+ Y^PT Flange Internal Shear Force Web 4 Figure 4.11: Upward Deflection. of the T beam - --u) 25 ^ szt T-beam 1 20 - T-beam 2 a) E s.; 2 m15 a) -c w oo cn 5 _ ^ U C) 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06 The Post-tensioning Force ( N ) Figure 4.12: Displacements of T-beam under One Patch Loading 1 E+07 Chapter 4. ANALYTICAL RESULTS ^ 61 4.2.3 Wheel loading analysis Four patches of vertical loading simulating a vehicle with four wheel were applied to the structure. According to the different locations of the four wheels we have three cases to study. CASE 1 The wheel loadings (each patch carried 500 mm x 300 mm x 0.1 (N/mm 2 ) = 15,000(N)) were at the centre lines of the T-beam 1 and T-beam 2 respectively. The separation in x-direction between the two patches was 2000mm. The structure is shown in Figure 4.13. PTB obtained the displacement and stress plots shown in Figure 4.14 to Figure 4.16. At the beginning ( PT = 0.0 ), the maximum bending stress in beam 1 and beam 2 were the same. All beams acted almost independent of each other, with little load sharing. Beam 3 did not share the wheel loading applied to beam 1 and beam 2. When the posttensioning forces were increased, the contribution from beam 3 to the sharing of loading increased, producing a corresponding difference between the maximum bending stresses of beam 1 and beam 2. For the maximum relative displacement between the adjacent flanges of the T-beams, the effect of the post-tensioning force was very obvious. At the beginning, the maximum relative displacement between the beam 1 and beam 2 was very small due to their similar loading and boundary conditions. But the maximum absolute deformation of the joint 2 in z-direction was very large. After the post-tensioning force got to 100,000 (N), the value of the maximum deformation of the joint element 2 in z-direction was negligible in comparison with that obtained at the PT = 0.0 as shown in Figure 4.16. Keeping the post-tensioning force at PT = 100, 000(N), the four patches of wheel loading were then moved from x = 0.0 to x = 6000mm. The Figure 4.17 and Figure 4.18 62 Chapter 4. ANALYTICAL RESULTS ^ PT L PT / Joint 1 1 ^ Joint 2 2^3 Figure 4.13: T-beam Structure CASE 1 give us the non-linear analytical results from PTB program. When the vehicle moved to midspan, the displacement and the bending stress of the T-beam achieved the peak values. CASE 2 The four patches of vehicle loading were acting at the edges of T-beam 1 and T-beam 2 as shown Figure 4.19. The deformations in joint elements were correspondingly more than those in CASE 1 at same level of the post-tensioning force. But the trend of the posttensioning force effect was the same. When PT approached 10,000 (N), the tightness of and load distribution in the structure was improved dramatically. The analytical results are shown from Figure 4.20 to Figure 4.22. The initial stiffnesses Ex and Ez of the springs have an influence on the magnitude of the relative displacements. Figure 4.23 shows the effect on the maximum relative ^ Chapter 4. ANALYTICAL RESULTS^ ni 0 ^.0 2 100 ^ E 6° co — c nsion 63 to 8^ 6 T-beam 1 4- T-begm 2 w co E) 2 T-beam 3 ^ri 0 ^ 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07 The Post-tensioning Force ( N ) 2 Figure 4.14: Bending Stress in T-beam Structure CASE 1 L. s ens' ni^or -g- 25 ^ cn^- T-beam 1 ^20a^ E 15 m 10 T-beam 2 T-beam 3 - - 5 to 0 g (+00 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07 The Post-tensioning Force ( N ) Figure 4.15: Displacements in T-beam Structure CASE 1 Chapter 4. ANALYTICAL RESULTS^ 20 ^ 0 A Joint Element 1 a) Joint Elgment 2 15 O 64 10 cn o a) > 5 'o; cr X 2 1E+01^1E+02^1E+03^1E+04 4+00 ^ 1E+05 The Post-tensioning Force ( N ) 1E+06^1E+07 Figure 4.16: Relative Displacements in T-beam Structure CASE 1 n dinC E io a) T-bepm 1 T-Own 2 T-b_eFn 3 I- 6 " 65 0) 4 C 2 '15 C CO cd 2 .... .. ... o o^0 0 ^ 0^1000^2000^3000^4000^5000^6000 Position of the Moving Loading(mm) Figure 4.17: Bending Stress under Moving Loading Chapter 4. ANALYTICAL RESULTS^ E 25 MS T-beam 1 T-beArri 2 T-beam 3 N .ro 20 46 65 15 E a) MI 5 a o^ 0 1000^2000^3000 4000^5000 Position of the Moving Loading(mm) Figure 4.18: Displacements under. Moving Loading PT L Joint 1^Joint 2 1 ^ 3 Figure 4.19: T-beam Structure CASE 2 6000 ^ Chapter 4. ANALYTICAL RESULTS^ a. a10 ^ a) T-beam 1 8 ^ T-ber 2 T-beam 3 6 •c 4 U) 2 2 ei)^ "ti 00^ a) co - 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07 The Post-tensioning Force ( N ) 2 Figure 4.20: Bending Stress in T-beam Structure CASE 2 St4t6titibiliff' ^25 ^b 66 ^ T-beam 1 T-beam 2 15 T-beam 3 a) 10 co 2 5 - U) 0 ^ g g ^ 2 11+00 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07 The Post-tensioned Force ( N ) Figure 4.21: Displacements in T-beam Structure CASE 2 Chapter 4. ANALYTICAL RESULTS^ 67 •••••■■• Joint Element 1 Joint Element 2 1E+01^1E+02^1E+03^1E+04 1E+05 The Post-tensioning Force ( N ) t3^ 1E+06^1E+07 Figure 4.22: Relative Displacements in T-beam Structure CASE 2 vertical displacement if the horizontal spring stiffness Ex and the vertical spring stiffness Ez are changed from Ex = 16.9N/mm and Ez = 17.4N/mm to Ex = 169N/mm and Ez = 174N/mm respectively. It is seen that these parameters should be known with some accuracy if the maximum relative displacement is to be estimated at intermediate posttensioning forces. However, the level of the post-tensioning force which makes the whole structure work as a unit is less affected by the uncertainty in Ex and E. Determination of Ex and Ex should include a larger experimental sample than used in this thesis. CASE 3 In this case, the vehicle loadings were next the center lines of the T-beams as shown Figure 4.24. In addition to the maximum bending stress, maximum displacement in the beam webs and the maximum deformation of the joint elements in the z-direction shown respectively in Figure4.25 through Figure4.27, the program PTB was used to study the •• Chapter 4. ANALYTICAL RESULTS^ 68 E E 60 Joint Element 1 c.) 50 a) _c Softer Spring Joint Element 1 40 Ex= 16.9 N/mm . 30 75. y 20 • Ex = 169 N/mm 76 10 E z = 174 N/mm _ a) CC Stiffer Spring E z = 17.4 N/mm - 1E+01^1E+02^1E+03^1E+04^1E+05 2^The Post-tensioning Force ( N ) 1E+06^1E+07 Figure 4.23: The Effect of the Stiffness of the Spring convergence of the solution as the number of terms in the Fourier series was increased. As discussed in Chapter 2, the larger the number of Fourier series terms included, the greater the coupling in the system of equations. Since the structure and the loading cases were symmetric in x-direction, only odd number terms in the Fourier series were included. From Figures 4.28 to 4.31, it can be seen that the convergence is quite good with only two Fourier terms (n = 1,3), being sufficient to obtain satisfactory answers. Chapter 4. ANALYTICAL RESULTS ^ 69 Figure 4.24: T-beam Structure CASE 3 VS. 0-10 T-beam 1 —aT-beam 2 QC — a) co co 6 •••••... T-beam 3 E 4 C) C C 2 co gQ 2 it+oo 1E+01^1E+02^1E+03^1E+04 1E+05^1E+06 The Post-tensioning Force ( N ) Figure 4.25: Bending Stress in T-beam Structure CASE 3 70 Chapter 4. ANALYTICAL RESULTS^ 25 ^ E T-beam 1 — T•beem 2 20,a^ T-beam 3 15 0 co 010 a. U, 0^ 1E+00 1E+01^1E+02^1E+03^1E+04 1E+05 ^ The Post-tensioning Force ( N ) 1E+06 Figure 4.26: Displacements in T-beam Structure CASE 3 E 35 ^ Joint Demerit 1 —tt— Joint Element 2 O 300^ 0 • _c 0 a u , 25. 20 15 O 10 Co T1 5 cc ii+oo • 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06 The Post-tensioning Force ( N ) Figure 4.27: Relative Displacements in T-beam Structure CASE 3 71 Chapter 4. ANALYTICAL RESULTS^ 0 8.5 E 8 a) co 7.5 (/) 7 ci) 6.5 co c 6 „ 5.5 co a+oo 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06 1E+07 The Post-tensioning Force ( N ) Figure 4.28: Effect of the Fourier Terms on Bending Stress E 20 as co 18 6" .06 a) 1 .0 14 - Ci 12 . 0 1Q i+cio 1E+01 1E+02^1E+03^1E+04 1 E+05 1E+06 The Post-tensioning Force ( N ) Figure 4.29: Effect of the Fourier Terms on Displacement 1E+07 Chapter 4. ANALYTICAL RESULTS ^ 72 C C) E 25 5 Number of Terms 7 20 1 15 "e5 a0 10 Te 3 0 5 . 5 0 1E+00 2 1E+01 1E+02^1E+03 1E+04 1E+05^1E+06 The Post-tensioning Force ( N ) Figure 4.30: Effect of the Fourier Terms on Relative Displacement in Joint 1 Number of Terms 7 —a1 3 5 1E+01^1E+02^1E+03^1E+04 The Post-tensioning Force ( N ) 1 E+05 1 E+06 The Effect of the Fourier Terms Figure 4.31: Effect of the Fourier Terms on Relative Displacement in Joint 2 Chapter 5 CONCLUSIONS A non-linear analysis for the post-tensioned timber bridge has been developed and implemented in the computer program PTB. The analysis takes into account the frictional contact between the bridge beams. The relative movements between beams, the stresses and deformations of the structure can be obtained as a function of the post-tensioning force by PTB. It is clear that post-tensioning forces increases the stiffness of the structure, improving load sharing and increasing the load carrying capacity of the structure. As part of a further study, it is suggested that more friction tests be done to confirm the 'same slip' assumption, and that reliability analyses be carried out to develop design criteria. The structural analysis presented here can form the basis for such probabilistic investigation. 73 Bibliography [1] Csagoly, Paul F., and Taylor, Raymond J. 1980. " A Structural Wood System for Highway Bridges ", IABSE Periodica 4/1980 IABSE Proceedings pp. 35-80. [2] Research and Development Division Ontario Ministry of Transportation and Communications, 1979. " Transverse Post-Tensioning of Longitudinally Laminated Timber Bridge Decks ", Research Report RR 220 [3] Verna, J.R., Graham, J.F., Shannon, P.H. 1984. " Timber Bridges: Benefits and Costs ", J. Structural Engr., ASCE, 110(7): 1563-1571 [4] Cheung, Y.K. 1968." Finite Strip Method in the Analysis of Elastic Plates with Two Opposite Simply Supported Ends ", Proc., Inst. Civ. Engrs., 40(5):1-7 [5] Foschi, Ricardo 0. 1982. " Structural Analysis of Wood Floor Systems ", Structural Division, Proceedings of ASCE, 108(ST7): 1557-1574 [6] Thompson, E.G., Goodman, J.R., and Vanderbilt, M.D. 1975. " Finite Element Analysis of Layered Wood Systems ", Journal of the Structural Division, ASCE, 101(12): 2659-2672 [7] Thompson, E.G., Vanderbilt, M.D., and Goodman, J.R., 1977. " FEAFLO: A Program for the Analysis of Layered Wood Systems ", 7: 237-248 [8] Taylor, Raymond J. and Csagoly, Paul F. 1980. " Rehabilitation of Wood Highway Bridges in Ontario Research Report Ministry of Transportation and Communications 74 Bibliography^ 75 [9] Onate, E. and Suarez, B. 1983. " A Unified Approach for the Analysis of Bridges, Plates and Axisymmetric Shells Using the Linear Mindlin Strip Element " Computers and Structures 407-426 [10] Harik, Issam E. and Salamoun, Ghassan L. 1988. " The Analytical Strip Method of Solution for Stiffened Rectangular Plates ", Computers and Structures 29(2): 283-291 [11] Timoshenko, S.P. and Woinowsky-Krieger,S., 1959."Theory of Plates and Shells", McGraw - Hill Book Co., New York, N.Y. [12] Zienkiewicz, 0.C., 1971." The Finite Element Method in Engineering Science ", McGraw-Hill Book Co., London,U.K. pp. 254-273 [13] " Wood Joist Floors ", 1978 Structural Research Report No.19, Civil Engineering Department, Colorado State Uni., Fort Collins, Colo., USA. [14] Cook, Robert D., Malkus,David S., Loesha,Michaele, 1989. " Concepts and Applications of Finite Element Analysis " Third Edition University of WisconsinMadison [15] Dhatt, Gouri., Touzot, Gilbert., 1984. " The finite Element Method Displayed ", Chichester [ West Sussex ] New York: Wiley Appendix A SHAPE FUNCTIONS AND DERIVATIVES Polynomial Shape Functions Vector {M0 } All components are zero except: mom= 42^1 M0( 2 ) = —8 (e — e e4 e) M0( 9 ) = 2 — 2 e + e) m0(10,1_42 +e M0(14) = E2^4E3^E41 1 M0( 15 ) —8 ( — e e3^e5) Vector {M3 } All components are zero except: M3 4e 1 ( 3 ) =^+ 76 e3 Appendix A. SHAPE FUNCTIONS AND DERIVATIVES^ M3( 4 ) = 8^ e + 63 - 6 4 ) M3(7) = 1 — 26.2 + 4.4 1 M3 (16) = — (36 + 46 2 — 6 3 — 26 4 ) 4 M3(17) 8 4^42 + + 4 4) Vector {M5 } All components are zero except: M5 (5) = M3 (3)^M5(6) = M3 (4)^M5(8) = M3 (7) M5 (18) = M3 (16)^M5(19) = M3 (17) Vector {M1} {M2 } {M4 } dM0 (k) dMg(k) ; M2(k) — M1(k) —^ de^de M4(k) d M3( k) ; dM5(k). M6(k) = de Where:^k = 1,2, ... ,19. Derivatives of the Displacement Functions aw(x,y)^N^ X = Etillo(4)1 T {8,2 }( n7 r )COS( nr ) az^n=1^L^L 77 78 Appendix A. SHAPE FUNCTIONS AND DERIVATIVES (92 w 1 ^E{iV10()}T{871}( 121-) 2 sin( 717) x2 9) n=1 aW n(X' Y) =^{Mi(0}T{bn}( ) Sin ( n7rX (-/Y^n=1^s^L N 2 nrx 52w (x ' Y) = E{M2()}T{Sn}( ) 2sin ( s^L ) 5 2^n=1 192W(X y) , 2 =E N{Mi()}T{Sn}(- )( OxOy r Ou(x,y)^_^ Ox n=1 1 1V1 au(x,y) \lT nrx )c°s( ) ^L^L fis l(L n )sn. 2 n=i 014(0/T{871}( )ws( (717 7 nrx av(x,y) N^nrx = E{m5(E)} {8n}(L^L )cos( ^ ) Ox^n=1 av(x,y) Oy N — 2^ nrx E{M6(0}T{8n}(-s )sin L n=1 d 2 U (x) N dx 2 n=1 d 2 W(x) dx 2 d 2 V(x) dx 2 d 2 0(x) dx 2 Un ( 117 ) 2 cos( n^ ) =^ 7 n=1 wn ( )2 sin (7 L-1- ) 2 sin(-1-7) —E 1%7 Vn ( -7 n=1 ) ) Appendix B THE PTB PROGRAM B.1 PTB Program Specifications The non-linear analysis for the transverse post-tensioned bridge structure can be done with the PTB computer program. The wheel loadings are idealized as patches of static loading. PTB evaluates the response of the structure subjected to the transverse posttensioning forces and the wheel loadings in terms of : • the maximum deflection of the T-beam web; • the maximum bending stress in the T-beam web; • the maximum deflection of the T-beam flange; • the maximum stress of the T-beam flange; • the maximum deformations (relative displacements ) of the joint elements in the x, the y and the z-direction. The stiffness matrices and the consistent load vectors were derived in Chapter 2. PTB can be run in an IBM 386 PC computer or a UNIX SUN work station. The program limits the size of the problem which can be solved, as follows: • MAX. NUMBER OF T-BEAMS = 10 ( MJT) • MAX. NUMBER OF JOINT ELEMENTS = 9 (MSP) 79 Appendix B. THE PTB PROGRAM^ 80 • MAX. NUMBER OF 3-D SPRINGS IN ONE JOINT ELEMENT = 126 (INCLUDING SPRINGS AT BOTH TOP AND BOTTOM) (MNP = 126 * 3 = 378 ) • MAX. NUMBER OF FOURIER TERMS = 4 ( MFT) (*) • MAX. NUMBER OF LOADED AREAS = 12 ( MLD) • MAX. NUMBER OF BOUNDARY CONDITIONS = 20 ( MBC) • MAX. NUMBER OF POST-TENSIONING FORCES PAIRS = 20 ( MPF) ^ (* ) If the problem is symmetric about midspan in the x-direction, 'four terms' in the Fourier series means that the terms with order 1,3,5 and 7 are chosen. If the problem is not symmetric in the x-direction, 'four terms' in the Fourier series means that the terms with order 1,2,3 and 4 are chosen. According to the capacity of the computer, the size of the problem can be expanded by changing the corresponding parameters in the PTB source code. B.2 Program Structure The program PTB consists of one main program and a number of subroutines. Each subroutine performs its specific function. The flow chart of the Main Program of PTB is given below. PTB uses common blocks to convey data between the main program and subroutines and between the subroutines as well. SUBROUTINE DATA This subroutine inputs data for the structure in free format through a data file PTB.DAT, which includes the size of the structure, the number of the terms in the 81 Appendix B. THE PTB PROGRAM^ PTB PROGRAM FLOW CHART / INPUT DATA / ASSEMBLE INCREMENTAL WHEEL LOADING VECTOR AND GLOBAL STIFFNESS MATRIX BOUNDARY CONDITIONS SOLVE EQUATIONS CORRESPONDING n th TERM - OF FOURIER SERIES TO GET DISPLACEMENTS OF THE STRUCTURE USING JACOBI ITERATION METHOD CHECK ALL THREE-DIMENSIONAL SPRINGS CALCULATE IMBALANCE LOAD VECTOR ASSEMBLE TANGENTIAL GLOBAL STIFFNESS MATRIX CALCULATE INCREMENTAL DISPLACEMENTS OF THE STRUCTURE ( STEP *) CHECK ALL THREE-DIMENSIONAL SPRINGS IF NO IF YES CHECK DISPLACEMENT CONVERGENCY IF NO IF YES / OUTPUT RESULTS / Figure B.1: PTB Program Flow Chart Appendix B. THE PTB PROGRAM^ 82 Fourier series, material properties, vehicle wheel loadings, transverse post-tensioning forces, number of three-dimensional springs in the joint elements, boundary conditions etc. This subroutine can print out the input data as a double check. SUBROUTINE GENMTX This subroutine calculates the integrals in Equation 2.20-Equation 2.27 using a sixpoint Gaussian quadrature to evaluate the contributions from the T-beam flanges to the stiffness matrix. SUBROUTINE C OVC ON This subroutine calculates the stiffness matrix of one T-beam element ESTIF(I,J) excluding the stiffness contributions of the T-beam web. ESTIF(I,J) is 19 x 19 symmetric matrix corresponding to each term of the Fourier series. SUBROUTINE STIF12 This subroutine creates the 12 x 12 stiffness matrix from one joint element with respect to each term of Fourier series. One joint element consists of 2 x N, three dimensional springs. N, has been input from the data file PTB.DAT. In this subroutine STIF12, subroutine BX, BY and BZ are called to evaluate the stiffness contributions from the x-, the y- and the z-direction springs according to the status of each spring respectively. SUBROUTINE CHECK This subroutine checks every spring in each joint element, according to the current displacement and the constitutive relationship, to calculate the stiffness matrix and the internal forces of the joint elements. SUBROUTINE LOAD This subroutine calculates the consistent load vector XF(n,IJ) according Equation 2.86 for wheel loading and Equation 2.90 and Equation 2.91 for post-tensioning forces. In XF(n,IJ) where: Appendix B. THE PTB PROGRAM^ 83 n = the n-th term in Fourier series IJ = 19 x NJT NJT = the number of the T-beam elements in the structure SUBROUTINE SLOAD According to the checking results obtained from SUBROUTINE CHECK, this subroutine creates the internal force vector components from the 'yielded' springs. SUBROUTINE ASMBLY This subroutine assembles the T-beam stiffness components, connector stiffness components and joint element stiffness components to form a global stiffness matrix by calling SUBROUTINE COVCON, and SUBROUTINE STIF12. From the Equation 2.96 we can note two kinds of matrices GK(n,n) and GK(n,m). The lower triangle of the diagonal global stiffness matrix GK(n,n) is stored column by column in the two dimensional matrix GSTIF(n,IJ). The components of the off-diagonal global stiffness matrix GK(n,m) are stored column by column in the two dimensional matrix ASTIF(nm,JK). SUBROUTINE GBC and SUBROUTINE ABC These two subroutines apply the boundary conditions to the stiffness matrices GSTIF(n,JK), load vector XF(n,IJ) and stiffness ASTIF(nm,JK). To impose a specific boundary restrain to the structure, the corresponding columns and rows in the stiffness matrices are set to zero, and the corresponding components in the load vector is also zeroed, as while the diagonal entry in GK(n,n) is set to 1. SUBROUTINE DCOMP GK(n,n) is a positive-definite matrix. It has a unique decomposition form GK(n,n) = [1,][L] T , where [L] is a lower triangular matrix with positive diagonal elements. This sub- routine decomposes GK(n,n) via the Cholesky method. The lower half band of GK(n,n), Appendix B. THE PTB PROGRAM^ 84 including the diagonal, is stored column by column in the GSTIF(n,JK). Cholesky decomposition is a useful method in solving the system linear equation like GK(n,n){u} = {Q} when GK (n, n) is a positive-definite banded matrix. SUBROUTINE SOLVE and SUBROUTINE ITERATION On the basis of Cholesky decomposition the SUBROUTINE SOLVE obtains the final results, according to Equation 2.106. Then the main program calls SUBROUTINE ITERATION to solve Equation 2.105. The final answer will be obtained when the iteration converges. li{Auk}!,i^< Eil{Auk}i7n12 k =1,2, ... ,N and E is the error tolerance SUBROUTINE IMBALANCE This subroutine calculates the imbalance load for each iteration in Newton-Raphson method. See Figure 2.8 SUBROUTINE RESULT This subroutine outputs all responses of the structure in terms: • the maximum displacement of the T-beam web • the maximum bending stress of the T-beam web • the maximum displacement of the T-beam flange • the maximum stress of the T-beam flange • the maximum deformations of the joint element in the x-direction, the y-direction and the z-direction In all cases, the maximum response is obtained by comparing 19 values at equally spaced locations along the longitudinal span. Appendix C PTB USER'S MANUAL PTB: Post-Tensioned Timber Bridge Non-linear Analysis Program User's Manual Version 1.0 by Cai Shen Department of Civil Engineering, UBC Vancouver, B.C. Canada V6T 1W5 November, 1991 85 Appendix C. PTB USER'S MANUAL^ 86 The computer program PTB Version 1.0, which is written in FORTRAN 77, is developed for implementation in a DOS-based microcomputer or UNIX SUN workstation. To run the program the user creates a data file named PTB.DAT in accordance with the input instructions given below. All data are to be entered in free format ( i.e. by providing a space or comma between each data entry ). The capacity limitations of the PTB program are: • MAX. NUMBER OF T-BEAMS = 10 ( MJT) • MAX. NUMBER OF JOINT ELEMENTS = 9 ( MSP) • MAX. NUMBER OF 3-D SPRINGS IN ONE JOINT ELEMENT = 126 (MNP = 126* 3 = 378 ) • MAX. NUMBER OF FOURIER TERMS = 4 ( MFT) • MAX. NUMBER OF LOADED AREAS = 12 ( MLD) • MAX. NUMBER OF BOUNDARY CONDITIONS = 20 ( MBC) • MAX. NUMBER OF POST-TENSIONED FORCES = 20 ( MPF) PTB produces two output files. PTB.OUT is a comprehensive output of the deformations and stresses in the structure and PTBB.OUT is a summary of the analysis results including the relative interface movements. Appendix C. PTB USER'S MANUAL^ 1. Enter: NTITLE NTITLE = problem title ( limited to 60 characters ) 2. Enter: NM, NJT, ISYM NM = maximum order of sine/cosine terms in the Fourier series. There can be maximum of four terms. For x-direction symmetrical problem, NM can be up to 7 ( the four terms will be of 1, 3, 5 and 7 ). If NM = 5, there will be three terms, these orders are 1, 3 and 5. For x-direction non-symmetrical problem, NM agrees with the number of terms. NM = 4 means terms with order 1, 2, 3 and 4 will be involved in the program. NJT = number of T-beams ( NJT < 10) ISYM 1^if the problem is x-direction symmetrical. 0^if the problem is not x-direction symmetrical 3. Enter: XL, SPJT XL = the span of the T-beam SPJT = T-beam spacing ( the width of the T-beam flange ) 4. Enter: BJT, HJT BJT = the width of the T-beam web HJT = the height of the T-beam web 87 Appendix C. PTB USER'S MANUAL ^ 88 5. Enter: EJT(I) EJT(I) = the elastic modulus of the i-th T-beam web 6. Enter: REG REG = the ratio of T-beam web elastic modulus to its shear modulus( E/G ) 7. Enter: PTK, EMUX, EMUZ PTK = the thickness of the T-beam flange EMUX = the x-direction friction coefficient of the T-beam flange EMUZ = the z-direction friction coefficient of the T-beam flange 8. Enter: PEX, PEY, PG, PVXY, PVYX PEX = the elastic modulus of the T-beam flange in x-direction PEY = the elastic modulus of the T-beam flange in y-direction PG = the shear modulus of the T-beam flange PVXY = Poisson's ratio of T-beam flange, strain in x-direction while stress in y-direction PVYX = Poisson's ratio of T-beam flange, strain in y-direction while stress in x-direction 9. Enter: SPC, XIN, CKPAL, CKPER, CKROT SPC = flange connector spacing along the T-beam span XIN = distance between the end support and the first connector along the T-beam span CKPAL, CKPER = Appendix C. PTB USER'S MANUAL ^ 89 connector load-slip moduli, respectively in the directions parallel and perpendicular to the T-beam ( unit = force/length ) CKROT = connector rotation modulus 10. Enter: TAOSX, TAOSZ TAOSX = the shear strength of the T-beam flange in x-direction TAOSZ = the shear strength of the T-beam flange in z-direction 11. Enter: NAN, N , MNL 8 NAN = 1^PTB will do non-linear analysis 0^PTB will do linear analysis only N, = the number of the top three-dimensional springs in one joint element ( the total number of the springs in one joint element is equal to N, * 2 * 3 < MNP = 378 ) MNL = the number of the incremental steps (each incremental loading is equal to the wheel loading divided by MNL i.e. PPWLD/MNL) 12. Enter: EX(I) EX(I) = the initial elastic modulus of the three dimensional spring in x-direction for the I-th joint element. Appendix C. PTB USER'S MANUAL^ 90 13. Enter: EY(I) EY(I) = the initial elastic modulus of the three dimensional spring in y-direction for the I-th joint element. 14. Enter: EZ(I) EZ(I) = the initial elastic modulus of the three dimensional spring in z-direction for the I-th joint element. 15. Enter: NWLD NWLD = the number of the wheel loadings; if NWLD = 0 goto step 17 16. Enter: JTLD(I), X1LD(I), X2LD(I), Y1LD(I), Y2LD(I), PPWLD(I) JTLD(I) = the number of the T-beam on which the I-th wheel loading is applied X1LD(I) = X2LD(I) = Y1LD(I) = Y2LD(I) = the coordinates of the I-th wheel loading patch, all coordinates are local Tbeam coordinates ( I = 1, NWLD ) 17. Enter: NPTF NPTF = the number of the transversely post-tensioning forces if NPTF = 0 goto step 20 Appendix C. PTB USER'S MANUAL^ 91 18. Enter: ES, AS, DIT, RFAC ES, AS = the elastic modulus and area of the post-tensioning cable DIT = the tooth distance of the threaded anchor head of post-tensioning tendon RFAC = factor to compensate the loss of post-tensioning force 19. Enter: XPTF(I), PTF(I) XPTF(I) = the x-coordinate of the I-th pair of post-tensioning force PTF(I) =the magnitude of the I-th pair of post-tensioning force ( I = 1, NPTF ) 20 Enter: NBC NBC = the number of applied boundary conditions if NBC = 0 skip step 21 21. Enter: IBC(I,1), IBC(I,2) IBC(I,1) =the boundary conditions for I-th T-beam IBC(I,2) = the number of the constraint ( 1 to 19 ) ( I = 1, NBC ) Appendix D PTB SOURCE CODE AND I/O FILE PTB.FOR PTB.DAT PTB.OUT PTBB.OUT 92 Appendix D. PTB SOURCE CODE AND I/O FILE^ D.1 PTB.FOR C C^Post Tensioned Timber Bridge Non-Linear Analysis Program C^ Version 1.0 (Nov. 1991)^ ,k C^ Civil Engineering Department C^ University of British Columbia C^ 2324 Main Mall C^ Vancouver, Canada V6T-1W5 C C C C C^ C^ PROGRAM LIMITATIONS C-C C 1) MAXIMUM NUMBER OF T-BEAMS^ = 10 ( MJT ) C 2) MAXIMUM NUMBER OF JOINT ELEMENTS ^= 9 ( MSP ) C 3) MAX. No. OF 3-DSPRINGS IN ONE JOINT ELEMENT = 126 (MNP 126*3 ) C 4) MAXIMUM NUMBER OF FOURIER TERMS ^= 4 ( MFT ) C 5) MAXIMUM NUMBER OF LOADED AREAS^= 12 ( MLD ) C 6) MAXIMUM NUMBER OF BOUNDARY CONDITIONS ^= 20 ( MBC ) C 7) MAXIMUM NUMBER OF POST-TENSIONED FORCES^= 20 ( MPF ) C C^ PROGRAM VARIABLES AND OPTIONS C^ C NTITLE = THE TITLE OF THE PROBLEM C NM^= MAXIMUM ORDER FOR THE SINE & COSINE SERIES. NM CAN BE UP C^TO 7 FOR SYMMETRIC PROBLEMS (ISYM = 1). NM CAN ONLY BE UP C^TO 4 FOR. NON-SYMMETRIC PROBLEMS. C NJT^= NUMBER OF T-BEAMS C ISYM^= 1 FOR SYMMETRIC ABOUT X-DIC CASE; = 0 OTHERWISE. C XL^THE SPAN OF THE BRIDGE C SPJT^= THE FLANGE WIDTH OF THE T-BEAM C BJT^= THE WEB WIDTH OF THE T-BEAM C EMT^= THE WEB HEIGHT OF THE T-BEAM C EJT(I) = THE WEB ELASTIC MODULUS OF THE i-th T-BEAM C R.EG^= THE RATIO OF THE WEB ELASTIC MODULUS TO ITS SHEAR MODULUS C PTK^= THE FLANGE THICK OF THE T-BEAM C EMUX^= THE FRICTION COEFFICIENT OF THE FLANGE IN X-DIRECTION C EMUZ^= THE FRICTION COEFFICIENT OF THE FLANGE IN Z-DIRECTION C PEX^= THE ELASTIC MODULUS OF THE T-BEAM FLANGE IN X-DIRECTION C PEZ^= THE ELASTIC MODULUS OF THE T-BEAM FLANGE IN Z-DIRECTION C PG^= THE SHEAR MODULUS OF THE T-BEAM FLANGE C PVXY^= POISSON'S RATIO, STRAIN IN X-DIC WHILE STRESS IN Y-DIC C PVYX^= POISSON'S RATIO, STRAIN IN Y-DIC WHILE STRESS IN X-DIC C SPC^= NAIL SPACING ALONG THE LONGITUDINAL SPAN OF THE BRIDGE C XIN^= DISTANCE BETWEEN THE END SUPPORT AND THE FIRST NAIL 93 Appendix D. PTB SOURCE CODE AND I/O FILE ^ C^ALONG THE LONGITUDINAL SPAN OF THE BRIDGE C CKPAL,CKPER= NAIL LOAD-SLIP MODULI, RESPECTIVELY, IN THE DIC. C^PARALLEL AND PERPENDICULAR TO THE WEB OF THE T-BEAM C CKROT = NAIL ROTATION MODULUS C TAOSX^= THE SHEAR STRENGTH OF THE FLANGE IN X-DIRECTION C TAOSZ^= THE SHEAR STRENGTH OF THE FLANGE IN Z-DIRECTION C NAN^= 1 - TO DO NON-LINEAR ANALYSIS C^0 - TO DO LINEAR ANALYSIS ONLY C NP^= THE HALF NUMBER OF THREE DIMENSIONAL SPRINGS IN ONE JOINT ELE. C EX(I)^= THE INITIAL ELASTIC MODULUS OF THE X-DIRECTION SPRING C EY(I)^= THE INITIAL ELASTIC MODULUS OF THE Y-DIRECTION SPRING C EZ(I)^= THE INITIAL ELASTIC MODULUS OF THE Z-DIRECTION SPRING C NWLD^= THE NUMBER OF THE WHEEL LOADINGS C JTLD(I) = THE NUMBER OF THE T-BEAM SUBJECTED BY THE WHEEL LOADING C C X1LD(I) C X2LD(I) = THE COORDINATES OF THE i-th WHEEL LOADING C Y1LD(I) C Y2LD(I) C C PPWLD(I) = MAGNITUDE OF THE WHEEL LOADING C NPTF^= THE NUMBER OF THE TRANSVERSELY POST-TENSIONING FORCES C ES,AS^= THE MODULUS AND AREA OF THE POST-TENSIONING CABLE C DIT^= THE TOOTH DISTANCE OF THE CABLE C RFAC^= THE C XPTF(I) = THE X COORDINATE OF THE i-th POST-TENSIONING FORCE C PTF(I) = THE MAGNITUDE OF THE i-th POST-TENSIONING FORCE C NBC^= NUMBER OF APPLIED BOUNDARY CONDITIONS. C IBC(I,1) = BOUNDARY CONDITIONS FOR I-TH T-BEAM WEB. C IBC(I,2) = NUMBER OF THE CONSTRAINT (1 TO 19). C IMPLICIT DOUBLE PRECISION (A - H, C) -Z) PARAMETER (MJT = 10, MEQ = 190, MSZ =-- 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B01/ RM00(19,19), RM22(19,19), RMO2(19,19), RM20(19,19), 1 RM33(19,19), RM66(19,19), RM36(19,19). RM63(19,19), 2 RM44(19,19), RM45(19,19), RM54(19,19), RM55(19,19), 3 RM11(19,19) COMMON /B02/ AA(19,19). RM(6,19). RM1(6,19) COMMON /B03/ ETA(6), GWT(6) COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RIT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PEW, PDG COMMON /B08/ CKPAL, CKPER, CKROT, SPC,XIN COMMON /B09/ X1LD(MLD), X2LD(MLD), YiLD(MLD), Y2LD(MLD), 1^PWLD(MLD).PPWLD(MLD), JTLD(MLD), NWLD COMMON /B10/ XPTF(MPF), PTF(MPF),NPTF 94 Appendix D. PTB SOURCE CODE AND I/O FILE ^ COMMON /B11/ IBC(MBC,2), NBC COMMON /B12/ ESTIF(19,19) COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ),FF(MFT,MEQ),RFF(MFT,MEQ),SVEC(MFT,MEQ), 3^SVEC1(MFT,MEQ),SVEC2(MFT,MEq),XFX(MFT,MEQ),EK(8,8) COMMON /B14/ Xl{(8.8).YK(8,8),2K(8,8), 1^BX1(1,8),BX2(1,8),BY1(1,8).BY2(1,8), 2^B21(1.8),B22(1,8).BX1T(8,1),BX2T(8,1),BY1T(8,1), 3^BY2T(8,1),B21T(8,1),B22T(8,1),SXK(8,8), 4^SYK(8,8),SZK(8,8),EK12(12,12) COMMON /B15/ EX(MSP),EY(MSP).EZ(MSP),NAN,NP,MNL COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,F2,TOLL,NNA COMMON /B18/ SIGMAY,WDX1(MSP),WDX2(MSP), 1^WDY1(MSP),WDY2(MSP),WDZ1(MSP),WDZ2(MSP) WRITE (*,111) 111 FORMAT (///////////////////////////////////) WRITE (",20) 20 FORMAT (T10, 2'*^ 3'"^PTB^ 4'w^ *',/T10, *',/T10. *',/T10, 5'* Post Tensioned Timber Bridge Non•Linear Analysis *',/T10, 6'"^PROGRAM^ *',/T10, Version 1.0^",/T10, 7'w^ *',/T10, 8'w^Civil Engineering Department^'',/T10, 9'*^University of British Columbia^*',/T10, Vancouver, Canada^*',/T10, November 1991 ",/T10, OPEN(UNIT=1,FILE='ptb.dat',STATUS='old') OPEN(UNIT=2,FILE='ptb.out',STATUS.'unknown') OPEN(UNIT=3,FILE='ptbb.ont',STATUS='unknown') PI = 4.0D0 DATAN(1.0D0) PI2 = PI**2 PI4 PI"4 TOL = .001 CALL DATA WRITE(*,1) 1 FORMAT(1X,'THE DATA HAVE BEEN INPUT') CALL GENMTX 95 Appendix D. PTB SOURCE CODE AND I/O FILE^ NPT 6 * NP DC) 25 I = 1.NST DO 25 J = 1,NPT ICH(I,J) = 1 25 ICH1(I,J) = 0 NNA 0 DO 26 IK =1,NM.NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM =(1+IK)/ 2 DO 26 I = 1,NEQ SVEC2(IKM,I) = 0.D0 SVEC(IKM,I) = 0.D0 RFF(IKM,I) = 0.D0 26 FF(IKM,I) = 0.D0 MMM = MNL + 1 DO 350 LMN =1,MMM C IF ( LMN .EQ. 1 ) THEN WRITE(*,*) WRITE(*,*)' APPLY POST-TENSIONING FORCE ONLY' END IF IF ( LMN .NE. 1 ) THEN WRITE(*,*) WRITE(*,*)' APPLY WHEEL LOADING TO THE BRIDGE' END IF C DO 3 I = 1, NWLD 3 PWLD(I) = (LMN-1)*PPWLD(I)/ (MMM 1) CALL LOAD(1.2) DO 30IK = 1,NM,NSTEP IKO =IK IF (NSTEP .EQ. 2) IK0 = (IK+1)/2 DC) 30 I = 1,NEQ 30 XFX(IKO,I) = XF(IKO,I) CALL ASMBLY WRITE(*,') 8 FORMAT(BX,'THE GLOBAL STIFFNESS MATRIX HAS BEEN ASSMBLED') DO 9 I = 1, NWLD 9 PWLD(I) = PPWLD(I)/ (MMM 1) WRITE(*,*) IF (LMN .EQ. 1) THEN IF (NPTF .EQ.0) GOTO 350 CALL LOAD(1,0) 96 Appendix D. PTB SOURCE CODE AND I/O FILE ^ END IF IF (LMN .GE. 2) THEN IF (NWLD .EQ.0) GOTO 350 CALL LOAD(0,2) END IF WRITE(*,io) 10 FoRMAT(SX,'THE GLOBAL LOAD MATRIX HAS BEEN ASSMBLED') WRITE(*,*) DO 29 IK =1.NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM =(1-f-IK)/2 DO 29 I = 1,NEQ 29 XF(IKM,I) = XF(IKM,I) RFF(IKM,I) DO 60 IK = 1, NM, NSTEP Do 55 IN = 1,IK,NSTEP IKM = IK INM = IN IF (NSTEP .EQ. 2) THEN IKM = (IK + 1)/2 INM = (IN + 1)/2 END IF C IF (IN .NE. IK) THEN IKO = (INM - 1) * NFT IKM^(INM 1)/2 END IF IF ( IN .EQ. IK ) THEN IKO = IKM CALL GBC(IKO) CALL DCOMP(IKO) Do 33 I = 1, NEQ FF(IKO,I) = XF(IKM,I) 33 CONTINUE CALL SOLVE(IKO) DO 44 I = 1,NEQ 44 SVEC(IKO,I) = XF(IKM,I) END IF 55 CONTINUE 60 CONTINUE C USING ITERATING METHOD TO SOLVE THE COUPLING EQUATION CALL ITERATE C IF ( NAN .EQ. o)G0To 315 C CALL CHECK(0) ^ WRITE(*,") WRITE( - ,*) CHECK HAS BEEN FINISHED' 97 Appendix D. PTB SOURCE CODE AND I/O FILE^ IF ( NNA .NE. 0 ) THEN WRITE(',")'^NON-LINEAR ANALYSIS IS DOING ' WRITE(*,*) DO 240 IK = 1, NM, NSTEP IKM = IK IF (NSTEP .EQ. 2) THEN IKM = (IK -I- 1)/2 END IF DO 240 I = 1, NEQ 240 SVEC1(IKM,I) = SVEC(IKM,I) C ^ 245 CALL IMBALANCE(1) CALL ASMBLY D() 249 IK = 1,NM,NSTEP DO 249 IN = 1,IK,NSTEF IKM = IK INM = IN IF(NSTEP .EQ. 2) THEN IKM = (IK-I-1)/2 INM = (IN-f-1)/2 END IF IF (IN .NE. IK) THEN IKO = (INM - 1) 'NFT IKM INM" (INM 1)/2 END IF IF ( IN .EQ. IK ) THEN IKO = IKM CALL GBC(IKO) CALL DCOMP(IKO) DO 246 I =1,NEQ FF(IKM,I) = XF(IKM,I) 246 CONTINUE CALL SOLVE(IKM) DO 297 I = 1,NEQ 247 SVEC(IKM,I) = XF(IKM,I) END IF 249 CONTINUE CALL ITERATE C NLF = 1 DO 260 IK = 1, NM, NSTEP IKM = IK IF (NSTEP .EQ. 2) THEN IKM = (IK + 1)/2 END IF TP = 0.0 TP1 = 0.0 AMAX = 0.0 DO 250 I = 1.NEQ TP1 = TP1 + (SVEC1(IKM,I)"*2)+SVEC2(IKM,I)*x2 TOLL = 0.01 98 Appendix D. PTB SOURCE CODE AND I/O FILE^ 250 TP = TP -1-(SVEC(IKM,I)"'2) WRITE(*,*)'TOLL = ',TOLL T1' = DSQRT(TP) TP1 = DSQRT(TP1) EP = TP/TP1 WRITE(*,")'EP = NON LINEAR',EP,'IKM = ',IKM WRITE(*,*)'LMN = ',LMN WRITE(",*)'PWLD(I) = ',(LMN-1)*PPWLD(1)/ (MMM 1) NPA = 1 IF ( EP .GT. TOLL) NPA = 0 NLF = NLF NPA DO 260 I = 1,NEQ SVEC1(IKM,I) = SVEC1(IKM,I)^SVEC(IKM,I) 260 CONTINUE DO 280 IK = 1, NM, NSTEP IKM = IK IF (NSTEP .EQ. 2) THEN IKM = (IK -I- 1)/2 END IF DO 280 I = 1, NEQ 280 SVEC(IKM,I) = SVEC1(IKM,I) IF ( NLF .EQ. 1) GOTO 300 WRITE(*,*) WRITE(*,*)'^NON-LINEAR ANALYSIS IS CONTINUING' WRITE(*,*) GOTO 245 END IF C IF ( NNA .EQ. 0)THEN WRITE(*,')'^NO NON-LINEAR SPRING OCCURS' WRITE(*,*) END IF C^ 300 CONTINUE CALL IMHALANCE(0) DO 310 I =1,NEQ RFF(IKM,I) = XF(IKM,I) 310 CONTINUE 315 CONTINUE DO 320 IK =1,NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM =(1-1-IK)/2 DO 320 I = 1,NEQ 320 SVEC2(IKM,I) = SVEC(IKM,I)^SVEC2(IKM,I) 350 CONTINUE D() 380 IK =1,NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM =(1+IK)/2 D() 380 I = 1,NEQ 380 SVEC(IKM,I) = SVEC2(IKM,I) c^CALL CHECK(1) 99 Appendix D. PTB SOURCE CODE AND I/O FILE ^ CALL RESULT WRITE(",398) 398 FORMAT(1X,'THE RESULTS HAVE BEEN OUTPUT') WRITE(*,*) CLOSE(1) CLOSE(2) CLOSE(3) STOP END C C SUBROUTINE DATA READS INPUT DATA FROM DATA FILE. C SUBROUTINE DATA IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RIT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PDV, PDG COMMON /B08/ CKPAL, CKPEB., CKROT, SPC,XIN COMMON /B09/ X1LD(MLD), X2LD(MLD), Y1LD(MLD), Y2LD(MLD), 1^PWLD(MLD),PPWLD(MLD), JTLD(MLD), NWLD COMMON /B10/ XPTF(MPF), PTF(MPF),NPTF COMMON /B11/ IBC(MBC,2), NBC COMMON /B15/ EX(MSP ).EY(MSP ),EZ(MSP ),NAN,NP,MNL COMMON /B17/ ICH(MSP.MNP),ICH1(MSP,MNP ),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA COMMON /B18/ SIGMAY,WDX1(MSP ),WDX2(MSP), 1^WDY1(MSP),WDY2(MSP),WDZ1(MSP),WDZ2(MSP) CHARACTER*72 NTITLE C INPUT DATA. READ (1,10) NTITLE 10 FORMAT (A) READ (1,*) NM, NJT, ISYM NST=NJT-1 READ (1,*) XL, SPJT READ (1,*) BJT, HJT READ (1,') (EJT(I),I=1,NJT), REG BETA = (1.0 - (0.63 BJT / HJT)) / 3.0 AJT = BJT HJT RIY = BJT * (HJT*'3) / 12.0 100 Appendix D. PTB SOURCE CODE AND I/O FILE^ RIZ = HJT (BJT"3) / 12.0 RIT = BETA * HJT (BJT"3) DO 20 I = 1, NJT GJT(I) = EJT(I) / REG 20 CONTINUE READ (1,*) PTK,EMUX,EMUZ READ (1,") PEX, PEY, PG, PVXY, PVYX PKX = (PEX * PTK**3) / (12.0D0 (1.0D0 (PVXY*PVYX))) PKY = PKX PEY / PEX PKV = PVXY PKX PKG = PG * (PTK"3) / 12.0D0 PDX = (PEX * PTK) / (1.0D0 - (PVXY*PVYX)) PDY = PDX * PEY / PEX PDV = PVXY * PDX PDG = PG * PTK READ (1,*) SPC, XIN,CKPAL, CKPER, CKROT C READ(1,*) TAOSX,TAOSZ READ(1,*) NAN,NP,MNL IF ( NAN .EQ. 0)THEN WRITE(*,*)' ** ONLY LINEAR ELASTIC ANALYSIS IS REQUIRED "' WRITE(*,*) END IF READ(1,*) (EX(I), I=1,NST) READ(1,') (BY(I), I=1,NST) READ(1,*) (EZ(I), I=1,NST) c Jan. 20,1992 c^ceexxx = ex(1) c^eeezaa^ez(1) c^alcak READ (1,") NWLD DO 150 I = 1, NWLD READ (1,*) JTLD(I), X1LD(I), X2LD(I), Y1LD(I), Y2LD(I), PPWLD(I) 150 CONTINUE READ (1,*) NPTF IF (NPTF .NE. 0) THEN READ (1,')ES,AS,DIT,RFAC DO 160 I = 1, NPTF READ (1,1 XPTF(I), PTF(I) 160 CONTINUE END IF READ (1,*) NBC IF (NBC .NE. 0) THEN D() 170 I = 1, NBC 170^READ (1,*) (IBC(I,J),J=1.2) 180 CONTINUE 101 Appendix D. PTB SOURCE CODE AND I/O FILE ^ END IF C DDXL = (5.0'9.8)*EMUX/EX(1) DDZL (5.0'9.8)*EMUZ/EZ(1) EXL = TAOSX *0.5"PTIC*XL/(NP-1) /DDXL EZL = TAOSZ '0.5"PTK'XL/(NP-1) /DDZL FXL = EXL DDXL FZL = EZL DDZL C IF (NPTF .NE. 0)THEN AW = PTK*XL/NPTF FS = PTF(1) SIGMAY = FS/AW FX = SIGMAY EMUX *0.5"PTK*XL/(NP-1) FZ = SIGMAY EMUZ *0.5*PTK'XL/(NP-1) IF ( FX .GT. FXL ) FX = FXL IF ( FZ .GT. FZL ) FZ = FZL DO 190 I = 1,NST EX(I)^EX(I)*SIGMAY*0.5*PTK'XL/(NP-1)/(5.0*9.8) EZ(I) = EZ(I)*SIGMAY*0.5"PTK*XL/(NP-1)/(5.0'9.8) IF ( EX(I) .GT. EXL ) EX(I) = EXL IF ( EZ(I) .GT. EZL ) EZ(I) = EEL c Jan. 20.1992 IF ( EX(I) .LT. EEEXXX ) EX(I) = EEEXXX IF ( EZ(I) .LT. EEEZZZ ) EZ(I) = EEEZZZ c 190 CONTINUE IF (EX(1) .NE. 0.0) DDX = FX/EX(1) IF (EZ(1) .NE. 0.0) DDZ = FZ/EZ(1) IF ( DDX .GT. DDXL ) DDX = DDXL IF ( DDZ .GT. DDZL ) DDZ = DDZL IF (EX(1) .EQ. 0.0) DDX = 0.0 IF (EZ(1) .EQ. 0.0) DDZ = 0.0 END IF C IF (NPTF .EQ. 0)THEN 102 Appendix D. PTB SOURCE CODE AND I/O FILE^ SIGMAY = 0.0 FX = 0.0 FZ = 0.0 DO 191 I = 1,NST EX(I) = 0.0 EZ(I) = 0.0 c Malcmagag.4.1(.****.W.K..IK*****WWW.71c.*.4.M.*..*MeMe.X.WA. c Jan. 20,1992 IF ( EX(I) .LT. EEEXXX ) EX(I) EEEXXX IF ( EZ(I) .LT. EEEZZZ ) EZ(I) = EEEZZZ 191 CONTINUE DDX = 0.0 DDZ 0.0 END IF C ^ C TO CALCULATE THE NO. OF THE TURN OF THE POST TENSIONING CABLE C ^ IF (NPTF .NE. 0)THEN TN = ( RFAC * FS * NJT SPJT ( 1/(ES'AS) 1/(PEY"AW) ))/DIT END IF C ^ C SET PARAMETERS FOR PROBLEM SIZE AND TYPE. C^ — NEg = NJT 19 LHB = 19 NSZ = 19 NEQ C NA]. = 144 * NST IF (ISYM .EQ. 0) THEN NSTEP = 1 NFT = NM ELSE IF (ISYM .EQ. 1) THEN NSTEP = 2 NFT = (NM + 1) / 2 END IF C ECHOING INPUT DATA. WRITE(2,200) NTITLE WRITE(3,200) NTITLE 200 FORMAT (/, lx, 26('''), ' FLOOR ANALYSIS PROGRAM *, 26 ('"'), //, ' PROBLEM TITLE: ', A) WRITE(2,210) NJT, NFT, NM WRITE(3,210) NJT, NFT, NM 210 FORMAT (/,' NUMBER OF FLOOR JOISTS 1^' NUMBER OF FOURIER TERMS USED 103 Appendix D. PTB SOURCE CODE AND I/O FILE ^ 2^• MAX. ORDER OF FOURIER TERM^= 14 ) WRITE(2,220) 220 FORMAT (//,' PROPERTIES AND DIMENSIONS OF JOISTS') WRITE(2.230) XL, SPJT, BJT, HJT =^E12.5, /, 230 FORMAT (' JOIST SPAN 1 ' JOIST SPACING = E12.5, /, 2 ' JOIST WIDTH = E12.5, /, 3 JOIST DEPTH = E12.5 ) Do 245 I = 1, NJT WRITE (2,240) I, EJT(I), GJT(I) 240 FORMAT JOIST NO. =', 13, 2X, • EJT =', E12.5, 2X, 1^GJT =', E12.5) 245 CONTINUE WRITE(2,250) 250 FORMAT (//,' PROPERTIES AND DIMENSIONS OF PLATE COVER') WRITE(2,260) PTK 260 FORMAT (' COVER THICKNESS = E12.5) WRITE(2,270) PKX, PKY, PKV, PKG 270 FORMAT (' KX = E12.5, 2X, ' KY = E12.5, /, 1^' KV = E12.5, 2X, KG = E12.5 ) WRITE(2,280) PDX, PDY, PDV, PDG 280 FoRMAT (' DX = E12.5, 2X, ' DY^E12.5, /, 1^DV^E12.5, 2X, DG = E12.5 ) WRITE(2,340) 340 FORMAT (//, PROPERTIES FOR CONNECTORS') WRITE(2,350) CKPAL, CKPER, CKROT 350 FORMAT (' STIFFNESS PARALLEL TO JOIST ^ E12.5, I, 1^' STIFFNESS PERPENDICULAR TO .10IST = E12.5, /, 2^ROTATIONAL STIFFNESS FLANGE/JOIST = E12.5 ) WRITE(2,360) SPC 360 FORMAT SPACING BETWEEN CONNECTORS ^E12.5, //) IF (NWLD .GT. 0) THEN WRITE(2.370) 370 FORMAT (' APPLIED TRANSVERSE LOADING') WRITE(2,380) 380 FoRMAT (' JOIST', 6X, 'Xi', 12X, 'X2', 12X, 'Y1', 12X, 'Y2', 1^11X, 'LOAD') DO 390 I = 1, NWLD 390 WRITE (2,400) JTLD(I), X1LD(I), X2LD(I), Y1LD(I), Y2LD(I), 1^PPWLD(I) 400 FORMAT (1X, 12, 5(2X,E12.5)) END IF IF (NPTF .GT. 0) THEN WRITE (2,410) 410 FoRMAT (' POST TENSIONING FORCES') WRITE (2,420) 104 Appendix D. PTB SOURCE CODE AND I/O FILE ^ 420 FORMAT (' X-LOC', 6X, 'FORCE') DO 430 I = 1, NPTF 430 WRITE (2,435) XPTF(I), PTF(I) 935 FORMAT (2(2X,E12.5)) END IF IF (NPTF .EQ. 0)THEN WRITE(2,') WRITE(2,*)'NO POST TENSIONING FORCE' END IF WRITE (2,440) NBC 440 FORMAT (//' NO. OF BOUNDARY CONDITIONS ', 14) IF (NBC .EQ. 0) GO TO 470 DO 950 I = 1, NBC 450 WRITE (2,960) (IBC(I,J)..1=1,2) 460 FORMAT (' JOIST NO. = ', 13, 3X, ' B.C. CODE NO. = 13) 970 CONTINUE C NPP = NP*2 WRITE(2,') WRITE(2,*)'THREE DIMENSIONAL SPRING BETWEEN TWO T. BEAMS',NPP WRITE(2,*) WRITE(2,1'STIFF. OF THE SPRING IN X-DIC.' DO 480 I = 1,NST 980 WRITE (2,490) I,EX(I) 490 FORMAT ('^JOINT NO. = ', 13, 3X, ' EX = E12.5) WRITE(2,1'STIFF. OF THE SPRING IN Y-DIC.' DO 500 I = 1,NST 500 WRITE (2,510) I,EY(I) 510 FORMAT ('^JOINT NO. = ', 13, 3X, ' EY = E12.5) WRITE(2,*)'STIFF. OF THE SPRING IN Z-DIC.' DO 520 I = 1,NST 520 WRITE (2,530) I,EZ(I) 530 FORMAT ('^JOINT NO. = ', 13, 3X, ' EZ = E12.5) WRITE(2,540)EMUX 590 FORMAT(1X,' FRICTION COEFFICIENT OF THE COVER (X-DIC.) =',E12.5) WRITE(2,543)EMUZ 543 FORMAT(1X,' FRICTION COEFFICIENT OF THE COVER (Z-DIC.) =',E12.5) IF (NPTF .NE. 0)THEN WRITE(2,700)ES WRITE(2,705)AS WRITE(2,710)DIT WRITE(2,715)RFAC WRITE(2,720)TN 700 FORMAT(1X,'STIFFNESS OF THE POSTENSIONING CABLE = ',E12.5) 705 FORMAT(1X,'THE AREA OF THE POSTENSIONING CABLE = ',E12.5) 710 FORMAT(1X,'THE DISTANCE BETWEEN THE CABLE TEETH = ',E12.5) 715 FORMAT(1X,'THE RELAX FACTOR OF THE CABLE = ',E12.5) 720 FORMAT(1X,'THE TURN NO. REQUIRED = ',E12.5) 105 ^ Appendix D. PTB SOURCE CODE AND I/O FILE ^ END IF RETURN END C SUBROUTINE GENTMX COMPUTES THE INTEGRALS REQUIRED IN THE DEVELOPMENT C OF THE ELEMENT STIFFNESS MATRIX NUMERICALLY USING SIX POINT GAUSS C QUADRATURE. C MCMCW*..,K*W**.747.(*.***..4..A(.1C...c*W7(..(*WWXDF11(.7..Y4.Wak...2K7(..**..,..f .*7K.Ma.” SUBROUTINE GENMTX IMPLICIT DOUBLE PRECISION (A - H, 0 Z) COMMON /B01/ RM00(19,19), RM22(19,19), RMO2(19.19), RM20(19,19), 1^RM33(19,19), RM66(19,19), RM36(19,19), RM63(19,19), 2^RM44(19,19). RM45(19,19), RM54(19,19), RM55(19,19). 3^RM11(19,19) COMMON /B02/ AA(19,19), RM(6,19), RM1(6,19) COMMON /B03/ ETA(6), GWT(6) DIMENSION ETA2(6), ETA3(6), ETA4(6), ETA5(6) CALL ZERO(RM, RM1) ETA(1) = 0.93246951420315200 ETA(2)^0.66120938646626500 ETA(3)^0-23861918608331900 ETA(4) -ETA(3) ETA(5) = -ETA(2) ETA(6) = -ETA(1) GWT(1) = 0.17132449237917000 GWT(2) = 0.36076157304813900 GWT(3) 0.46791393457269100 GWT(4) = GWT(3) GWT(5) = GWT(2) GWT(6) = GWT(1) DO 10 I = 1, 6 ETA2(I) = ETA(I) ETA(I) ETA3(I) = ETA2(I) ETA(I) ETA4(I) = ETA3(I) * ETA(I) ETAS(I) = ETA4(I) ETA(I) 10 CONTINUE C MO AND M2 MATRICES. DO 20 I = 1, 6 RM(I,1) = ETA2(I) - (5.0D0 ETA3(I) / 4.000) 1^(ETA4(I) / 2.000) + (3.000 ETAS(I) / 4.000) RM(I,14) = ETA2(I)^(5.000 " ETA3(I) / 4.000) 1^- (ETA4(I) / 2.000) - (3.000 ETA5(I) / 4.000) RM(I,10) = 1.000 - (2.0D0 " ETA2(I)) -I- ETA4(I) RM(I,2) = (ETA2(I) - ETA3(I) ETA4(I) ^ETAS(I)) / 8.0D0 RM(I,15)^(-ETA2(I) ETA3(I)^ETA4(I)^ETA5(I)) / 8.000 RM(I,9) = (ETA(I) - (2.000 * ETA3(I)) ^ETA5(I)) / 2.000 106 ^ Appendix D. PTB SOURCE CODE AND I/O FILE ^ RM1(I,1)^2.0D0 - (15.0130 ETA(I) / 2.000) 1^- (6.000 ETA2(I)) + (15.000 ETA3(I)) RM1(I.14)^2.000 + (15.000 ETA(I) / 2.000) 1^- (6.0D0 ETA2(I)) - (15.0130 * ETA3(I)) RM1(I,10) = -4.0D0 + (12.0D0 * ETA2(I)) RM1(I,2) = (2.0130 - (6.000 * ETA(I)) - (12.000 ETA2(I)) 1^-f (20.000 ETA3(I))) / 8.0D0 RM1(I,15) = (-2.000 - (6.0130 * ETA(I)) + (12.000 ETA2(I)) 1^+ (20.000^ETA3(I))) / 8.0130 RM1(I,9) = ((-12.0D0 ETA(I)) + (20.0130 ETA3(I))) / 2.000 20 CONTINUE CALL DMAT(1) DO 40 I = 1, 19 DO 30 = 1, 19 RMOO(I,J) = AA(I,J) 30 CONTINUE 40 CONTINUE CALL DMAT(2) DO 60 I = I, 19 DC) 50 J = 1, 19 RM22(I,J) = AA(I,J) 50 CONTINUE 60 CONTINUE CALL DMAT(3) D080 I = 1, 19 DO 70 = 1, 19 RMO2(I,J) = AA(I,J) 70 CONTINUE 80 CONTINUE CALL DMAT(4) DO 100 I = 1, 19 DO 90 .1 = 1, 19 RM20(I,J) = AA(I,J) 90 CONTINUE 100 CONTINUE CALL ZERO(RM, RM1) C M3 AND M6 MATRICES. DO 110 I = 1, 6 RM(I,3) = (-3.000 * ETA(I) / 4.000) + ETA2(I) 1^+ (ETA3(I) / 4.000) - (ETA4(I) / 2.0D0) RM(I.16) = (3.000 ETA(I) / 4.000) + ETA2(I) 1 (ETA3(I) / 4.0130)- (ETA4(I) / 2.000) RM(I,7) = 1.000 - (2.000 * ETA2(I)) + ETA4(I) RM(I,4) = (-ETA(I) ETA2(I) + ETA3(I) ETA4(I)) / 8.000 RM(I,17) = (-ETA(I) - ETA2(I) ETA3(I) + ETA4(I)) / 8.000 RM1(I,5) = (-3.000 / 4.000)^(2.0D0 ETA(I)) 1^+ (3.000 * ETA2(I) / 4.0130) - (2.0D0 * ETA3(I)) RM1(I,18) = ( 3.000 / 4.000) 4 (2.000 ETA(I)) 1^- (3.000 * ETA2(I) / 4.0130) - (2.000 ETA3(I)) RM1(I,8) = (-4.000 * ETA(I)) + 4.000 ETA3(I) 107 Appendix D. PTB SOURCE CODE AND I/O FILE ^ RM1(I,6) = (-1.0D0^(2.01)0 ETA(I))^(3.0D0 ETA2(I)) 1^(4.0D0 ETA3(I))) / 8.01)0 RM1(L19)^(-1.0D0 (2.0D0 ETA(I))^(3.0D0 ETA2(I)) 1^(4.0D0 * ETA3(I))) / 8.01,0 110 CONTINUE CALL DMAT(1) DO 130 I = 1, 19 DO 120 J = 1,19 RM33(I,J) = AA(I,J) 120 CONTINUE 130 CONTINUE CALL DMAT(2) DO 150 I = 1, 19 DO 140 J^1, 19 RM66(I,J) = AA(I,J) 140 CONTINUE 150 CONTINUE CALL DMAT(3) DO 170 I = 1, 19 DO 160 J = 1, 19 RM36(I,J) = AA(I,J) 160 CONTINUE 170 CONTINUE CALL DMAT(4) DO 190 I = 1, 19 DO 180 J = 1. 19 RM63(I,J)^AA(I,J) 180 CONTINUE 190 CONTINUE C M5 AND M4 MATRICES. DO 200 I = 1.6 RM(I,5) = RM(I,3) RM(I,3) = 0.0D0 RM(I,18) = RM(I,16) RM(I,16)^0.0D0 RM(I,8) = nm(I,T) Rm(i, r) = 0.0D0 , RM(I,6) = RM(I,4) RM(I,4) = 0.0D0 RM(I,19) = RM(I,17) RM(I,17) = 0.0D0 RM1(1,3) = RM1(I,5) RM1(I,5) = 0.0D0 I1M1(1.16) = RMI(I,18) RMI(I,18) = 0.0D0 = RM1(I,8) RM1(I,8) = 0.0D0 RM1(I,4) = RM1(I,6) RM1(I,6) = 0.0D0 RM1(I,17) = RM1(I,19) 108 ^ Appendix D. PTB SOURCE CODE AND I/O FILE^ itm1(I,19) = 0.000 200 CONTINUE CALL DMAT(2) DC) 220 I = 1, 19 DC) 210 J = 1, 19 RM44(I,J) = AA(I,J) 210 CONTINUE 220 CONTINUE CALL DMAT(4) DC) 240 I = 1, 19 DC) 230 .1 = 1, 19 RM45(I,J) = AA(I,J) 230 CONTINUE 240 CONTINUE CALL DMAT(3) DC) 260 I = 1, 19 DC) 250 .1 = 1. 19 RM54(I,J) = AA(I,J) 250 CONTINUE 260 CONTINUE CALL DMAT(1) DO 280 I = 1, 19 DO 270 .1 = 1, 19 RM55(I,J) = AA(I,J) 270 CONTINUE 280 CONTINUE CALL ZER.O(RM, RM1) C M1 MATRIX. DC) 290 I = 1, 6 RM(I,1) = (2.000 * ETA(I)) - (15.000 ETA2(I) / 4.000) 1^- (2.000 * ETA3(I)) + (15.000 ETA4(I) / 4.000) RM(I,14) = (2.000 * ETA(I)) + (15.000 ETA2(I) / 4.000) 1^(2.0D0 * ETA3(I)) . (15.000 * ETA4(I) / 4.000) RM(010) = (•4.000 * ETA(I)) + (4.000 ETA3(I)) RM(I,2) = (( 2.0D0 ETA(I)) . (3.0D0 ETA2(I)) 1^(4.0D0 ETA3(I)) + (5.0D0 ETA4(I))) / 8.000 RM(I,15) = ((•2.000 * ETA(0) • (3.0D0 * ETA2(I)) 1^+ (4.000 * ETA3(I)) + (5.0D0 ETA4(I))) / 8.000 RM(I,9) = (1.000 - (6.000 ETA2(I)) 1^+ (5.0D0 ETA4(I))) / 2.000 290 CONTINUE CALL DMAT(1) DO 310 I = 1, 19 DO 300 .1 = 1, 19 RM11(I,J) = AA(I,J) 300 CONTINUE 310 CONTINUE RETURN END 109 Appendix D. PTB SOURCE CODE AND I/O FILE^ C SUBROUTINE ZERO CALLED BY GENMTX. C SUBROUTINE ZERO IMPLICIT DOUBLE PRECISION (A - H, 0 - Z) COMMON /B02/ AA(19,19), RM(6,19), RM1(6,19) DO 20 I = 1,6 D 0 10 .1 = 1, 19 RM(I,J)^0.0D0 RM1(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE RETURN END C SUBROUTINE DMAT CALLED BY GENMTX. c SUBROUTINE DMAT (II) IMPLICIT DOUBLE PRECISION (A - H, 0 - Z) COMMON /B02/ AA(19,19), RM(6,19), RM1(6,19) COMMON /B03/ ETA(6), GWT(6) DO 20 I = 1, 19 DO 10 .1 = 1, 19 AA(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE DO 50 I = 1, 6 DO 40 J = 1, 19 Fl = RM(I,J) IF (II .EQ. 2 .OR. II .EQ. 4) F1 = R.M1(I,J) IF (F1 .NE. 0.0) THEN DO 30 K = 1, 19 F2 = RM1(I,K) IF (II .EQ. 1 .OR. II .EQ. 4) F2 = RM(I,K) IF (F2 .NE. 0.0) THEN AA(J,K) = AA(J,K) + Fl * F2 GWT(I) END IF 30^CONTINUE END IF 40 CONTINUE 50 CONTINUE RETURN END C SUBROUTINE COVCON COMPUTES THE STIFFNESS MATRIX FOR THE PLATE COVER C AND CONNECTORS ASSOCIATED WITH A GIVEN T-BEAM STRIP. C^ mak M 1 SUBROUTINE COVCON (IN,IK) 11C.A.*.aRaleatalellt...”..ara, **.311■31C.*”....11(*.W.R 110 Appendix D. PTB SOURCE CODE AND I/O FILE^ IMPLICIT DOUBLE PRECISION (A H, 0 - Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, Moz = 1296) PARAMETER (MFT = 4. MLD = 12, MDC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B01/ RM00(19,19), RM22(19,19), RMO2(19,19), RM20(19,19), 1 RM33(19,19), RM66(19,19), RM36(19,19), RM63(19,19), 2 RM44(19.19), RM45(19,19), RM54(19,19), RM55(19,19), 3 RMii(19,19) COMMON /B03/ ETA(6), GWT(6) COMMON /B04/ NM, NSTEP, NFT, NJT.NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RIT COMMON /B07/ PTK, PEX, PEy, PG, PVXy, PVYX, PKX, PKy, PKV, PKG, PDX, PDY, PDV, PDG COMMON /B08/ CKPAL, CKPER, CKROT, SPC,XIN COMMON /B12/ ESTIF(19,19) COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA DIMENSION AD(19,19) DO 20 I = 1, 19 DO 10 J = 1, 19 ESTIF(I,.I) = 0.0 10 CONTINUE 20 CONTINUE C CONTRIBUTION FROM THE PLATE COVER. FAC = 0.0 IF (IN .EQ. IK) FAC = PKX (IK**4) PI4 * Sp.TT / (4.0*XL* - 3) DO 60 I = 1, 19 DO 50 J = 1, 19 AD(I,J) = ItMoo(I,J) 50 CONTINUE 60 CONTINUE CALL ADD(AD, FAC) C FAC = 0.0 IF (IN EQ. IK) FAC PKY * 4.0 XL / (SPJT**3) DO 100 I = 1, 19 DO so = 1, 19 AD(I,J) = RM22(I,J) 90 CONTINUE 100 CONTINUE CALL ADD(AD, FAC) C FAC 0.0 IF (IN Eq. IK) FAC = -PKV PI2 / (SPJT*XL) FAC1 FAC (IK"2) DO 140 I = 1, 19 DO 130 .1 = 1, 19 111 Appendix D. PTB SOURCE CODE AND I/O FILE^ AD(I,J) = KMO2(I,J) 130 CONTINUE 140 CONTINUE CALL ADD(AD, FAC1) C FAC1 = FAC " (IN"2) DO 160 I = 1. 19 DO 160 J = 1, 19 AD(I,J) = 150 CONTINUE 160 CONTINUE CALL ADD(AD, FAC1) C FAC = 0.0 IF (IN .EQ. IK) FAC = PKG 4.0 " PI2 * (IK"2) / (SPJT"XL) DO 200 I = 1, 19 DO 190 J = 1, 19 AD(I,J)^RM11(I,J) 190 CONTINUE 200 CONTINUE CALL ADD(AD, FAC) C PAC = 0.0 IF (Ili .EQ. IK) FAC = PDX SPJT PI2 * (IK**2) / (4.0*XL) DO 240 I = 1. 19 DO 230 J^1, 19 AD(I,J) = RM33(I,J) 230 CONTINUE 240 CONTINUE CALL ADD(AD, FAC) C FAC 0.0 IF (IN .EQ. IK) FAC = PDY * XL / SPJT DO 280 I = 1, 19 D0270 .1 = 1, 19 AD(I,J) = KM66(I,J) 270 CONTINUE 280 CONTINUE CALL ADD(AD, FAC) C FAC = 0.0 IF (IN .EQ. IK) FAC = -PDV *PI 2.0 FAC1 FAC " IK DO 320 I = 1, 19 DO 310 J= 1, 19 AD(I,J) = RM36(I,J) 310 CONTINUE 320 CONTINUE CALL ADD(AD, FAC1) C FAC1 =PAC IN D() 340 I = 1, 19 112 Appendix D. PTB SOURCE CODE AND I/O FILE^ DO 330 J = 1, 19 AD(I,J) = RM63(I,J) 330 CONTINUE 340 CONTINUE CALL ADD(AD, FAC1) C FAC = 0.0 IF (IN .EQ. IK) FAC = PDG. XL SP.1T / 4.0 FAC1 = FAC 4.0 / (SPJT*'2) DO 380 I = 1, 19 DO 370 J = 1,19 AD(I,J) = RM44(I,J) 370 CONTINUE 380 CONTINUE CALL ADD(AD, FAC1) C FAC1 = FAC 2.0 " IN PI / (SPJT*XL) DO 400 I = 1, 19 DO 390 .1 = 1, 19 AD(I,J) = RM45(I,J) 390 CONTINUE 400 CONTINUE CALL ADD(AD, FAC1) C FAC1 = FAC * 2.0 " IK * PI / (SPJT"XL) DO 420 I = 1, 19 DO 410 J = 1, 19 AD(I,J) = 11M54(I,J) 410 CONTINUE 420 CONTINUE CALL ADD(AD, FAC1) C FAC1 = FAC IK * IN PI2 / (XL*"2) DO 440 I = 1, 19 DO 430 J = 1, 19 AD(I,J) = RM55(I,J) 430 CONTINUE 440 CONTINUE CALL ADD(AD, FACT) C CONTRIBUTION FROM THE COVER-TO-STIFFENER CONNECTORS. DC) 460 I = 1, 6 DO 450 J = 1, 19 AD(I,J) = 0.0 450 CONTINUE 460 CONTINUE AD(1,7) = 1.0 AD(1,11) = -1.0 AD(1,10) = -IK " PI (HJT t PTK) / (2.0"XL) AD(2,7) = 1.0 AD(2,11) = -1.0 113 Appendix D. PTB SOURCE CODE AND I/O FILE^ AD(2,10) = -IN * PI * (HIT + PTK) (2.0 *X L) AD(3,8) = 1.0 AD(3,12) = -1.0 AD(3,9) = -PTK / (2.0*SPJT) AD(3,13) = -HIT / (2.0*SPJT) AD(4,8) = 1.0 AD(4,12) = -1.0 AD(4,9) = -PTK / (2.0"SPIT) AD(4,13) = -HIT / (2.0*SPJT) AD(5,9) = 1.0 / SPIT AD(5,13) = -1.0 / SPIT AD(6,9) = 1.0 / SPIT AD(6,13) = -1.0 / SPJT C IF (IN NE. IK) GO TO 510 SPAR = XL / (2.0*SPC) SPER = XL / (2.0'SPC) SROT = SPER GO TO 520 510 SPAR = 0.0 SPER = 0.0 SROT = 0.0 520 SPAR = SPAR * CKPAL SPER = SPER CKPER SROT = SROT CKROT C DO 540 I = 1, 19 DC) 630 J = 1, 19 ESTIF(I,J) = ESTIF(I.J) + (SPAR AD(1,I) AD(2,J)) 1^ (SPER * AD(34) * AD(4.J)) 2^ (SROT AD(5,I) AD(6,J)) 530 CONTINUE 540 CONTINUE RETURN END C C FUNCTION Z IS USED TO CALCULATE Ix( n,m,i) alea....acakag C FUNCTION Z(N, M, Z1, ZO, XL, NX, PI) IMPLICIT DOUBLE PRECISION (A - H, 0 -Z) IF (N .EQ. M) GO TO 20 S1 DSIN((N M)*PI*Z1/XL) S2 = DSIN((N MrPrZO/XL) S3 = DSIN((N M)*PI*Z1/XL) S4 = DSIN((N M)*PrZO/XL) DI = (Si - S2)/ (2.0*(N M)) D2 = (S3 - S4) / (2.0"(N+M)) IF (NX .EQ. 1) GO TO 10 Z = XL "(D1- D2) /PI RETURN 10 Z = XL * (Dl + D2 )/PI 3r. 114 Appendix D. PTB SOURCE CODE AND I/O FILE^ RETURN 20 S1 DSIN(2.0"N*PrZ1/XL) S2 = DSIN(2.0"N*PI'ZO/XL) DI = PI ( Z1 - Z0)/(2.0"XL) D2 = (S1 - S2) / (4.0*N) IF (NX .EQ. 1 ) GO TO 30 Z = XL * (DI - D2)/PI RETURN 30 Z = XL * (DI + D2)/PI RETURN END C C SUBROUTINE ADD, WHICH IS CALLED BY SUBROUTINE COVCON, BUILDS THE C STIFFNESS AND MASS MATRICES ASSOCIATED WITH THE PLATE COVER C AND THE CONNECTORS FOR A GIVEN T-BEAM STRIP. C SUBROUTINE ADD (AD, FAC) IMPLICIT DOUBLE PRECISION (A - H, 0 Z) COMMON /B12/ ESTIF(19,19) DIMENSION AD(19,19) IF (FAC .EQ. 0.0D0) RETURN DO 20 I^1, 19 DC) 10 J^1, 19 ESTIF(I,J) = ESTIF(I,J) + (FAC AD(I,J)) 10 CONTINUE 20 CONTINUE RETURN END C C SUBROUTINE ASMBLY ASSEMBLES THE DIAGONAL STIFFNESS SUB-MATRICES C FOR EACH FOURIER. TERM. GSTIF(IKO,JK) AND ASTIF(IKO,JK) C WW*Xc* ,..4.2..******.MC***7C21k.,F**21 , W.W.N[X, **..W.WWW*...*Mak*...,./.1C.**.” SUBROUTINE ASMBLY IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B06/ EJT(MJT), GJT(MJT), BJT, OUT, AJT, RIY, RIZ, RIT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PDV, PDG COMMON /B11/ IBC(MBC,2), NBC COMMON /B12/ ESTIF(19,19) COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ).FF(MFT,MEQ )..RFF(MFT,MEQ),SVEC(MFT,MEQ), 3^SVEC1(MFT,MEQ ),SVEC2(MFT,MEQ ),XFX(MFT,MEQ),EK(8,8) COMMON /B14/ XK(8,8),YK(8,8),ZK(8,8), 115 Appendix D. PTB SOURCE CODE AND I/O FILE ^ 1^BX1(1.8),DX2(1,8).BY1(1,8),BY2(1,8), 2^BZ1(1,8),BZ2(1,8),13X1T(8.1),BX2T(8,1),BYIT(8,1), 3^BY2T(8,1),BZ1T(8,1),HZ2T(8.1),SXK(8,8), 4^SYK(8,8),SZK(8,8),EK12(12.12) COMMON /B15/ EX(MSP),EY(MSP),EZ(MSP).NAN,NP,MNL COMMON /B17/ ICH(MSP.MNP),ICH1(MSP.MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL.NNA Do 60 IK = 1, NM, NSTEP DO 55 IN = 1,IK,NSTEP IKM = IK INM = IN IF (NSTEP .EQ. 2) THEN IKM = (IK + 1)/2 INM = (IN + 1)/2 END IF C IF (IN .NE. IK) THEN IKO = (INM 1) 'NFT IKM 'NM' (INM + 1)/2 DO 5 I = 1, NA1 5 ASTIF(IKO,I)= 0.0 DO 20 IE = 1,NJT IJ^(IE 1)*19 IF (IE .LE. NST ) THEN CALL STIF12(IE,IN.IK) C^FBC39.F NOV. 8,1991 C IF ( NBC .NE. 0 ) THEN CALL ABC(IKO,IE) END IF IJ = ( IE - 1 ) * 12 Do 17 J= 1,12 DO 17 K = 1,12 JK = ( IJ J 1 ) * 12 K ASTIF(IKO,JK) = ASTIF(IKO,JK) EK12(J,K) 17 CONTINUE C END IF 20 CONTINUE END IF IF ( IN .EQ. IK ) THEN IKO = IKM FAC2 = (IK**2) PI2 / (2.0D0 * XL) FAC4 = (IK**4) PI4 / (2.oD0 * (XL"3)) 116 Appendix D. PTB SOURCE CODE AND I/O FILE^ DO 25 I = 1, NSZ GsTiF(IK0,i) = o.oDo 25 CONTINUE C CONTRIBUTION FROM THE PLATE COVER AND CONNECTORS. CALL COVCON (IN,IK) Do 50 IE = 1, NJT IJ = (IE -^19 DO 90 = 1, 19 DO 30 K = 1, JK = LHB (IJ + K - 1) + J K + 1 GSTIF(IKO,JK) = GSTIF(IKO,JK) + ESTIF(J,K) 30^CONTINUE 40 CONTINUE C CONTRIBUTION FROM THE JOISTS. JK = LHB (IJ + 9) + 1 GSTIF(IKO,JK) = GSTIF(IKO,JK) + (EJT(IE) * RIY * FAC4) JK = LHB * (IJ + 10) + 1 GSTIF(IKO,JK) = GSTIF(IKO,JK) + (EJT(IE) * AJT * FAC2) JK = LHB (IJ + 11) + 1 GSTIF(IKO,JK) = GSTIF(IKO,JK) + (EJT(IE) " RIZ * FAC4) JK LHB^+ 12) + 1 GSTIF(IKO,JK) = GSTIF(IKO,JK) + (GJT(IE) RIT FAC2 / (sPJT*'2)) IF (IE .LE. NST ) THEN CALL STIF12(IE,IN,IK) IJ=(IE-1)*19+13 DO 45 .1=1,12 DO 45 K=i,J JK=LHB - (IJ+K-1)+J-K+1 GSTIF(IKO,JK) = GSTIF(IKO,JK) + EK12(J,K) 45 CONTINUE END IF 50 CONTINUE END IF 55 CONTINUE 60 CONTINUE RETURN END C C THIS SUBROUTINE GENERATES THE STIFNESS MATRICES OF THE JOINT C ELEMENT BETWEEN THE T-BEAM STRIP. C SUBROUTINE STIF12 (IE,IN,IK) IMPLICIT DOUBLE PRECISION (A - H, 0 - Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF 20) 117 Appendix D. PTB SOURCE CODE AND I/O FILE ^ PARAMETER (MSP = 9, MFO = 6, MNP 378) COMMON /B00/ PI, PI2. PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ.NA1 COMMON /B05/ XL, SPJT COMMON 0307/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY. PDV, PDG COMMON /1313/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ)2F(MFT,MEQ),RFF(MFT,MEQ),SVEC(MFT,MEQ), 3^SVEC1(MFT,MEQ),SVEC2(MFT,MEQ),XFX(MFT,MEQ),EK(8,8) COMMON /B14/ XK(8,8),YK(8.8),ZK(8,8), 1^BX1(1.8),BX2(1,8),13Y1(1.8),BY2(1,8), 2^13Z1(1,8),BZ2(1,8),13X1T(8,1),BX2T(8,1),BY1T(8,1), 3^BY2T(84),BZ1T(8,1),BZ2T(8,1),SXK(8,8), 4^SYK(8,8),SZK(8,8),EK12(12.12) COMMON /1315/ EX(MSP).EY(MSP ),EZ(MSP),NAN,NP,MNL COMMON /1317/ ICH(MSP,MNP).ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA M=IK N=IN DO 20 J=1,8 BX1(1,J)=0.0 BX2(1,J)=0.0 BY1(1,J)=0.0 BY2(1,J)=0.0 BZ1(1,J)=0.0 13Z2(1,3)=0.0 BX1T(J,1)=0.0 BX2T(J,1)=0.0 13Y1T(3,1)=0.0 BY2T(J,1)=0.0 BZ1T(J,1)=0.0 BZ2T(.7,1)=0.0 20 CONTINUE CALL BX(N,M,IE) CALL BY(N,M,IE) CALL BZ(N,M,IE) DO 60 1=1,8 DO 60 J=1,8 EK(I,J)=SXK(I,J)+SYK(I,J)+SZK(I,J) 60 CONTINUE CALL EXP RETURN END C C THE CONTRIBUTION OF THE SPRING IN X DIC. CALLED BY STIF12 C akalt”.“cTeacww..z...71ta.k..x......e.acsetkak.a.kale.**.alc.*NFWaRalc.Mle.**7.30**Te.e. SUBROUTINE BX (N,M,IE) IMPLICIT DOUBLE PRECISION (A - H, 0 - Z) PARAMETER (MIT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) 118 Appendix D. PTB SOURCE CODE AND I/O FILE ^ PARAMETER (MFT = 4, MLD = 12, MBC 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /1100/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PDV, PDG COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MF0.MOZ), 1^XF(MFT,MEQ).FF(MFT,MEQ),RFF(MFT,MEQ).SVEC(MFT,MEQ), 3^SVEC1(MFT,MEQ),SVEC2(MFT,MEQ),XFX(MFT,MEQ),EK(8,8) COMMON /B14/ XK(8,8),YK(8.8),ZK(8,8), 1^BX1(1,8),BX2(1,8),BY1(1,8).BY2(1,8), 2^BZ1(1,8),BZ2(1,8),BXIT(8,1),BX2T(8,1),BY1T(8,1), 3^BY2T(8,1),BZ1T(8,1),BZ2T(8,1),SXK(8,8), 4^SYK(8,8),SZK(8,8),EK12(12.12) COMMON /B15/ EX(MSP),EY(MSP),EZ(MSP),NAN,NP,MNL COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA DO 10 1=1,8 DO 10 .1=1,8 10 SXK(I,J)=0.0 DO 100 I=1,NP DX=XL/(NP-1) IF ( I.Eq.1 .0R. I.EQ.NP ) DX=0.5*XL/(NP-1) XI=(I-1)*XL/(NP-1) BX1(1.1)=PTK/2.*(M*PI/XL)*DCOS(M*PrXI/XL) BX1(1,3)=-DCOS(M - PI*XI/XL) BX1(i,5)=-PTK/2.*(M . PI/XL)*DCOS(M*PrXI/XL) BX1(1,7)=DCOS(M*PI*XI/XL) KX1 = (I - 1)*6 + 1 C BX1T(1,1)=PTK/2,*(N*PI/XL)*DCOS(N*PrXI/XL)*EX(IE) 1^*ICH(IE,KX1) BXIT(3,1)=-DCOS(N- PrXI/XL)*EX(IE) 1^*ICH(IE,KX1) BX1T(5,1)=-PTK/2.'(N*PI/XL)'DCOS(N'PrXI/XL)*EX(IE) 1^*ICH(IE,KX1) BX1T(7,1)=DCOS(N'PP, XI/XL)*EX(IE) 1^'ICH(IE.KX1) BX2(1,1)=-PTIC/2.*(M'PI/XL)*DCoS(M - PrXI/XL) BX2(1,3)=-DCOS(WPI'XI/XL) BX2(1,5)=PTK/2.*(WPI/XL)*DCOS(M'PrXI/XL) BX2(1,7)=DCOS(M*PrX1/XL) KX2 (I .^+ 2 C BX2T(1,1)=-PTK/2.'(1,PTI/XL) - DCOS(N*PI*XI/XL)*EX(IE) 1^*ICH(IE,KX2) 119 Appendix D. PTB SOURCE CODE AND I/O FILE^ BX2T(3,1)=-DCOS(N*PI*XI/XL)*EX(IE) 1^*ICH(IE,KX2) BX2T(5,1)=PTK/2.*(N*PI/XL)*DCOS(N*PI*X1/XL)*EX(IE) 1^*ICH(IE.KX2) BX2T(7,1)=DCOS(N*PrXI/XL)*EX(IE) 1^.ICH(IE,KX2) CALL DGMULT (BX1T,BX1,XK,8,1,8,8,1,8) DO 20 11=1,8 DC) 20 3=1,8 20 SXK(I1,J)=SXK(I1,3)-I.XK(I1,3) CALL DGMULT (BX2T,BX2,XK,8,1,8,8,1,8) DO 30 11=1,8 DO 30 3=1,8 30 SXK(II,J)=SXK(11,J)+XK(I1,J) 100 CONTINUE RETURN END C C THE CONTRIBUTION OF THE SPRING IN Y DIC. CALLED BY STIF12 C SUBROUTINE BY (N,M,IE) IMPLICIT DOUBLE PRECISION (A - 11, 0 - Z) PARAMETER (MIT 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD 12, MBC 20, MPF = 20) PARAMETER (MSP = 9, MF() = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG. 1^PDX, PDY, PDV, PDG COMMON /B13/ GSTIF(MFT.MSZ). ASTIF(MFO,MOZ), 1^XF(MFT,MEQ),FF(MFT,MEQ ),FLFF(MFT,MBQ ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ),SVEC2(MFT,MEQ ),XFX(MFT,MEQ ),EK(8,8) COMMON /B14/ XK(8,8),YK(8,8),ZK(8,8), BX1(1,8),I3X2(1,8),BY1(1,8),BY2(1,8), 2^BZ1(1,8),BZ2(1,8),BX1T(8,1)33X2T(8,1),BY1T(8,1), 3^BY2T(8,1),BZ1T(8,1),BZ2T(8,1).SXK(8,8), 4^SYK(8,8),SZK(8,8),EK12(12,12) COMMON /B15/ EX(MSP),EY(MSP ),EZ(MSP),NAN,NP,MNL COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA DO 10 1=1,8 DO 10 J=1,8 10 SYK(I,J)=0.0 DO 100 I=1,NP DX=XL/(NP-1) IF ( I.EQ.1 .OR. I.EQ.NP ) DX=0.5'XL/(NP-1) XI=(I-1).XL/(NP-1) 120 Appendix D. PTB SOURCE CODE AND I/O FILE ^ BY1(1.2)=PTK/(2.*SPJT)*DSIN(M*PI*XI/XL) BY1(1,4)=-DSIN(M*PI*XI/XL) BY1(1,6)=-PTK/(2.*SPJT)*DSIN(M'PrXI/XL) BY1(1,8)=DSIN(M*PI"XI/XL) KY1 = (I - 1) *6 + 3 IF(ICH(IE,KY1) .EQ. 1) EKY1^1.00*EY(IE) IF(ICH(IE.KY1) .EQ. -1) EKY1 = 1.00 1.00E+10 IF(ICH(IE,KY1) .EQ. 0 ) EKY1 = 0.00 C BY1T(2,1)=PTK/(2."SPJT)*DSIN(N'PI*XI/XL) 1^"EKY1 BY1T(4,1 )=-DSIN(N'PI*XI/ XL) 1^*EKY1 BY1T(6,1)=-PTK/(2."SPJT)*DSIN(N'PrXI/XL) 1^*EKY1 B Y1 T(8,1 )=DSIN(NWI'XI/ XL) 1^*EKY1 BY2(1,2)=-PTK/(2.*SPIT)*DSIN(M*PI"XI/XL) BY2(1,4)=-DSIN(M*PI*XI/XL) BY2(1,6)=PTK/(2."SPIT)*DSIN(M"PI'XI/XL) BY2(1,8)=DSIN(M*PI*XI/XL) KY2 = (I - 1) *6 +4 IF(ICH(IE,KY2) .EQ. 1) EKY2 = 1.00*EY(IE) IF(ICH(IE,KY2) .EQ. -1) EKY2 = 1.00 1.00E+10 IF(ICH(IE,KY2) .EQ. 0 ) EKY2 = 0.00 C BY2T(2,1)=-PTK/(2.*SPJT)*DSIN(WPI*XI/XL) 1^*EKY2 BY2T(4,1)=-DSIN(N*PrXI/XL) 1^*EKY2 BY2T(6,1)=PTK/(2.'SPJT)*DSIN(WPI*XI/XL) 1^*EKY2 BY2T(8,1)=DSIN(N*PrXI/XL) 1^*EKY2 CALL DGMULT (BY1T.BY1,YK,8,1,8,8,1,8) DO 20 11=1,8 DC) 20 J=1,8 20 SYK(I1,J)=SYK(I1,J)+YK(I1,J) CALL DGMULT (BY2T,BY2,YK,8,1,8,8,1,8) DO 30 11=1,8 DO 30 J=1.8 30 SYK(I1,J)=SYK(II,J)+YK(I1,J) 100 CONTINUE RETURN END c C THE CONTRIBUTION OF THE SPRING IN Z DIC. CALLED BY STIF12 121 Appendix D. PTB SOURCE CODE AND I/O FILE ^ C SUBROUTINE BZ (N,M,IE) IMPLICIT DOUBLE PRECISION (A - H, 0 - Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ 1296) PARAMETER (MFT = 4, MLD = 12. MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NET, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B07/ PTK, PEX. PEY, PG, PVXY. PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PDV, PDG COMMON /B13/ GSTIF(MFT,MSZ), ASTIE(MEO,MOZ), 1^XF(MFT,MEQ),17P(MFT,MEQ ),RFF(MFT,MEQ),SVEC(MPT,MEQ ), 3^SVEC1(MFT,MEQ),SVEC2(MET,MEQ).XFX(MFT,MEQ),EK(8,8) COMMON /B14/ XK(8,8),YK(8,8),ZK(8,8), 1^BX1(1.8),BX2(1,8).BY1(1,8),BY2(1,8), 2^BZ1(1,8),BZ2(1.8),DX1T(8,1),BX2T(8,1),BY1T(8.1), 3^BY2T(8,1),BZ1T(8,1),HZ2T(8,1),SXK(8.8), 4^SYK(8,8),SZK(8,8),EK12(12,12) COMMON /B15/ EX(MSP ),EY(MSP),EZ(MSP ),NAN,NP,MNL COMMON /B17/ ICH(MSP,MNP).ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA DO 10 1=1,8 DO 10 3=1,8 10 SZK(I,J)=0.0 DO 100 I=1,NP DX=XL/(NP-1) IF ( I.EQ.1 .0R. I.EQ.NP ) DX=0.5*XL/(NP-1) XI=(I-1)*XL/(NP-1) BZ1(1,1)=-DSIN(M*PI*XI/XL) BZ1(1,5)=DSIN(M*PI'XI/XL) KZ1 = (I -1) * 6 + 5 C BZ1T(1,1)=-DSIN(N*PI*XI/XL)*EZ(IE) 1^*ICH(IE,KZ1) BZIT(5.1)=DSIN(N*PI*XI/XL)*EZ(IE) 1^*ICH(IE,KZ1) CALL DGMULT (BZ1T,BZ1,ZK,8,1,8,8,1,8) DO 20 11=1,8 DO 20 J=1,8 20 SZK(I1,J)=SZK(I1,J)+ZK(11,J) 100 CONTINUE C DO 200 I=1,NP DX=XL/(NP-1) IF ( I.EQ.1 .OR. I.EQ.NP ) DX=0.5*XL/(NP-1) XI=(I-1)*XL/(NP-1) 122 Appendix D. PTB SOURCE CODE AND I/O FILE ^ B Z 2 ( 1 ) D SIN (M"P XI / X L ) B Z 2 (1 ,5 )= D SIN( M*P I'XI/ XL) KZ2^(I -1) * 6 -4- 6 C B Z2T(1,1)=-DSIN(N*PrXI/XL)*E Z(IE ) 1^*1CH(IE,KZ2) B Z2T(5,1)=DSIN(N*PI*XI/XL)*EZ(IE) 1^*ICH(IE,KZ2) CALL DGMULT (BZ2T,HZ2,ZK,8,1,8,8,1,8) DO 21 11=1,8 DO 21 J=1,8 21 SZK(I1,..1)=SZK(I1,J)+ZK(11,J) 200 CONTINUE C RETURN END C *************************************************************** C THIS SUBROUTINE EXPANDS EK(8,8) TO EK12(12,12) CALLED BY STIF12 C SUBROUTINE EXP IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER. (MFT = 4, MLD = 12, MBC = 20, MI'? = 20) PARAMETER (MSP = 9, MFO = 6. MNP = 378) COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ),FF(MFT,MBQ ),I1FF(MFT,MEQ ),SVEC(MFT,MEQ), 3^SVEC1(MFT,MEQ),SVEC2(MFT,MEQ ),XFX(MFT,MBQ ),EK(8,8) COMMON /B14/ XK(8,8),YK(8,8),ZK(8,8), 1^BX1(1,8).BX2(1.8),BY1(1,8),BY2(1,8), 2^BZ1(1,8).BZ2(1,8).BX1T(8,1),BX2T(8,1),I3Y1T(8,1), 3^BY2T(8,1),BZ1T(8.1),BZ2T(8,1),SXK(8,8), SYK(8,8),SZK(8,8),EK12(12,12) DO 10 1=1,12 DO 10 J=1,12 10 EK12(I,J)=0.0 DO 100 J=1,3 DO 20 1=1,3 20 EK12(I,J)=EK(I,J) EK12(5,J)=EK(4,J) DC) 40 1=7,9 15=1-2 90 EK12(I,J)=EK(I5,J) EK12(11,J)=EK(8,J) 100 CONTINUE DC) 120 1=1.3 120 EK12(I,5)=EK(I,9) DO 130 1=5 123 Appendix D. PTB SOURCE CODE AND I/O FILE^ 130 EK12(5,5)=EK(4,4) DO 140 1=7,9 15=1-2 140 EK12(I.5)=EK(15,4) EK12(11,5)=EK(8,4) DO 300 .1=7,9 .75=J-2 DO 220 1=1,3 220 EK12(I,J)=EK(I,J5) EK12(5,J)=EK(4,..75) DO 240 1=7,9 15=1-2 240 EK12(I,J)=EK(I5,J5) EK12(11,J)=EK(8,J5) 300 CONTINUE DO 320 1=1,3 320 EK12(I,11)=EK(I,8) EK12(5,11)=EK(4,8) DO 340 1=7,9 15=1-2 340 EK12(I,11)=EK(I5,8) EK12(11,11)=EK(8,8) RETURN c END C TO MULTIPLY A REC. M BY N MATRIX BY ANOTHER RECT. C N BY L MATRIX. c SUBROUTINE DGMULT(A. B, C, M, N, L. NA. NB, NC) IMPLICIT DOUBLE PRECISION (A - H, O -Z) REAL*8 A(NA.1), B(ND,1), C(NC,1) DO 20 .1 = 1, L DO 20 I = 1, M C(I,J)^0.D0 DO 10 K = 1, N 10^C(I,J) = C(I,J)^A(I,K) B(K,J) 20 CONTINUE RETURN END c C SUBROUTINE LOAD COMPUTES THE GLOBAL LOAD VECTOR. c SUBROUTINE LOAD(III,IIII) IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, M13C = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL 124 Appendix D. PTB SOURCE CODE AND I/O FILE^ COMMON /B03/ ETA(6), GWT(6) COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1 COMMON /B06/ XL, SPJT COMMON /B09/ X1LD(MLD), X2LD(MLD), YlLD(MLD), Y2LD(MLD), 1^PWLD(MLD),PPWLD(MLD). JTLD(MLD), NWLD COMMON /B10/ XPTF(MPF), PTF(MPF),NPTF COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ).FF(MFT,MEQ ),RFF(MFT.MEQ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ ),SVEC2(MFT,MEQ ),XFX(MFT.MEQ),EK(8,8) COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP).SL(12), 1^DDX,DDZ,FX,FZ,TOLL.NNA DIMENSION W(19) DO 20 IK = 1, NM, NSTEP IKO = IK IF (NSTEP .EQ. 2) IKO (IK +1) / 2 DO 10 I = 1, NEQ XF(IKO,I) = 0.0D0 10 CONTINUE 20 CONTINUE C CONTRIBUTION FROM TRANSVERSE LOADING IF (IIII .EQ. 2) THEN IF (NWLD^0) THEN Do 70 IL = 1, NWLD IE = JTLD(IL) 1.1 = (IE - 1)^19 DO 60 IK = 1, NM, NSTEP IK0 = IK IF (NSTEP .EQ. 2) IKO = (IK +1) / 2 PINL = rI IK / XL DO 30 I = 1, 19 W(I) = 0.0D0 30^CONTINUE XIX = 2.0D0 Y1LD(IL) / SPJT XIY 2.0D0 Y2LD(IL) / SPJT DO 40 I = 1, 6 XI = XIX + (XIY - XIX) * (1.0D0 + ETA(I)) / 2.0D0 XI2 = XI**2 XI3 = XI2 * XI XI4 = XI3 XI XI5 = XI4 XI W(1) = W(1) + (XI2 (5.0D0 XI3 4.0D0) - (XI4 / 2.0D0) + (3.0130 XI5 / 4.0D0)) GWT(I) W(2) W(9) = W(2) + (XI2 - X13 - XI4 + XI5) GWT(I) / 8.0D0 = W(9) + (XI - (2.0D0 * XI3) + XI5) * GWT(I) / 2.0D0 W(10) = W(10) + (1.0D0 (2.0D0 XI2) + XI4) GWT(I) W(14) = W(14) + (XI2 + (5.0D0 * XI3 / 4.0D0) 1^(XI4 / 2.0D0) (3.0D0 XI5 / 4.0D0)) GWT(I) 125 Appendix D. PTB SOURCE CODE AND I/O FILE ^ W(15) = W(15) + (-XI2 - X13 + XI4 + XIS) G^T(I) / 8.0D0 40^CONTINUE Y1 = YSLD(IL) Y2 Y2LD(IL) XI X1LD(IL) X2 = X2LD(IL) FACC = (1-0D0 / (2.0D0 PINL)) * (Y2LD(IL) Y1LD(IL)) (DCOS(PINL*X1LD(IL)) DCOS(PINL*X2LD(IL))) W(1) = W(1) FACC W(2) = W(2) * FACC W(9) = W(9) FACC W(10) = W(10) FACC W(14) = W(14) FACC W(15) = W(15) * FACC DO 50 I = 1, 19 = XF(IKO,I.1+I) + ((PWLD(IL) W(I))) 50^CONTINUE 60 CONTINUE 70 CONTINUE END IF END IF C CONTRIBUTION FROM POST TENSIONING FORCES IF (III .EQ. 1) THEN IF (NPTF .NE. 0) THEN DO 90 I = 1. NPTF DO 80 IK = 1, NM, NSTEP IKO = IK IF (NSTEP -EQ. 2) IKO = (IK +1) / 2 PINL = PI IK / XL V5 = DSIN(PINL*XPTF(I)) IJ = 5 XF(IKO,IJ) = XF(IKO,IJ) -f (PTF(I) V5) V18 = DSIN(PINL*XPTF(I)) IJ^((NJT - 1) * 19) + 18 XF(IKO,IJ) = XF(IKO,IJ) - (PTF(I) V18) 80 CONTINUE 90 CONTINUE END IF END IF RETURN END C C THIS SUBROUTINE CALCULATES THE EQUIVALENT LOADS CONTRIBUTED C FROM THE NON-LINEAR SPRINGS C SUBROUTINE SLOAD(IE,IK) 126 ^ Appendix D. PTB SOURCE CODE AND I/O FILE^ IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ.NA1 COMMON /B05/ XL, MIT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PDV, PDG COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,ME Q ).FF(MFT,ME Q ),RFF(MFT.ME Q ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ ),S VE C2(MFT,ME Q ),XFX(MFT,MEQ),EK(8,8) COMMON /B15/ EX(MSP ),EY(MSP ),EZ(MSP ),NAN,NP,MNL COMMON /B11/ ICH(MSP,MNP ).ICH1(MSP.MNP ).SL(12), 1^DDX.DDZ,FX,FZ,TOLL,NNA DO 10 I = 1,12 10^SL(I) = 0.0 DX = XL/(NP - 1) DO 100 IP = 1,NP IF (IP .EQ. 1 .OR. IP .EQ. NP) Dx=0.5-xL/(Np-1) XI = (IP - 1)*XL/(NP - 1) KX1 = (IP - 1)*6 + 1 KX2 = (IP - 1).6 + 2 KZ1 = (IP - 1)-6 + KZ2 = (IP - 1)*6 + 6 SL(1) = SL(1) 1^PTK/2.*(IK"PI/XL)'DCOS(IK*PI*XI/XL)*FX'ICH1(IE,KX1) 2^- PTK/2.*(IK*PI/XL)*DCOS(IK'PrXI/XL)*FX*ICH1(IE,KX2) 3^- DSIN(IIC*PI*XI/XL)*FVICH1(IE,KZ1) ^3 ^DSIN(IIC*PrXI/XL)*FZ*ICH1(IE,KZ2) SL(3) = SL(3) 1^DCOS(IK*PI*XI/XL)*FX*ICH1(IE,KX1) 2^DCOS(IK*PI*XI/XL)*FX*ICH1(IE,KX2) SL(7) = SL(7) 1^• PTK/2.'(IK*PI/XL)*DCOS(IK'PI*XI/XL)*FX*ICH1(IE,KX1) 2^PTK/2.'(IK*PI/XL)*DCOS(IK*PrXI/XL)xFX'ICH1(IE,KX2) 3^DSIN(IK'PI'XI/XL)*FZ'ICH1(IE,KZ1) 3^DSIN(IIC*PI*XI/XL)*FVICH1(IE,KZ2) SL(9) = SL(9) 1^DCOS(IICTI'XI/XL)*FX'ICH1(IE,KX1) 2^DCOS(IK*PVXI/XL)*FX.ICH1(IE,KX2) 100 CONTINUE RETURN END 127 Appendix D. PTB SOURCE CODE AND I/O FILE ^ C SUBROUTINE GBC APPLIES THE BOUNDARY CONDITIONS TO THE STIFFNESS MATRIX C GSTIF(IKO,JK) AND LOAD VECTOR XF(IKO,JK). C SUBROUTINE GBC (IKO) IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON 1B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1 COMMON /B11/ IBC(MBC,2), NBC COMMON /1313/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), XF(MFT,ME(4),FF(MFT,MEQ),RFF(MFT,MEQ),SVEC(MFT,MEQ), 3^SVECI(MFT,MEQ),SVEC2(MFTNEQ),XFX(MFT,MEQ),EK(8,8) IF (NBC.NE.0) THEN DO 30 I = 1, NBC NE = IBC(I,1) IDOF = IBC(I,2) M = (NE - 1) 19 IDOF C =1 MM = M - 1 IF (M .GE. 19) .11 = M 18 DO 10 .1 = .11, MM JK = (LHB - 1)^- 1) M GSTIF(IKO,JK) = 0.0 10 CONTINUE MI = M 1 Ms = M 4- 18 IF (M8 .GT. NEQ) MS = NEQ Do 20 .1 = MI, M8 JK = (LHB - 1) (M 1) J GSTIF(IKO,JK) = 0.0 20 CONTINUE JK = (LHB - 1) * (M 1) M GSTIF(IKO,JK) = 1.0 XF(IKO,M) = 0.0 30 CONTINUE END IF RETURN END C C SUBROUTINE ABC APPLIES THE BOUNDARY CONDITIONS TO THE STIFFNESS MATRIX C ASTIF(IKO,JK) . C SUBROUTINE ABC (IKO,IE) IMPLICIT DOUBLE PRECISION (A - H. o Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1 128 Appendix D. PTB SOURCE CODE AND I/O FILE^ COMMON /B11/ IBC(MBC,2), NBC COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ),FF(MFT,MEQ ).RFF(MFT.MEC4 ),SVEC(MFT.MEg ), 3^SVE Cl(MFT,ME ),SVEC2(MFT,MEQ ).XFX(MFT.MEQ ),EK(8.8) COMMON /B14/ XK(8,8),YK(8,8),ZK(8,8), 1^DX1(1,8),BX2(1,8),BY1(1,8),BY2(1,8), 2^B Z1(1.8 ),B Z2(1,4 ),B X1 T(8,1 ).B X2T(4 ),B Y1 T(8,1 ). 3^BY2T(8,1),BZ1T(4,1),BZ2T(4,1),SXK(8,4), 4^SYK(8,8),SZK(8,8),EK12(12,12) IF (NBC.NE.0) THEN DO 20 I = 1, NBC NE = IBC(I,1) IDOF^IBC(I,2) C IF (NE .EQ. IE) THEN IF (IDOF .GE. 14) THEN J= IDOF - 13 DO 10 K = 1,12 10 EK12(J,K) = 0.0 K = IDOF - 13 DO 14 J^1,12 14 EK12(J.K) = 0.0 END IF END IF IF ( NE .EQ. IE--1 ) THEN IF ( IDOF .LE. 6 ) THEN .I.1 = IDOF + 6 DO 16 KK = 1,12 16 EK12(JJ,KK) = 0.0 KK = IDOF + 6 DO 18 JJ = 1,12 18 EK12(JJ,KK) = 0.0 END IF END IF 20 CONTINUE END IF RETURN END c C SUBROUTINE DCOMP DECOMPOSES THE STIFFNESS DIAGONAL SUB-MATRICES 129 ^ Appendix D. PTB SOURCE CODE AND I/O FILE^ C USING THE CHOLESKY METHOD. C SUBROUTINE DCOMP (IKO) IMPLICIT DOUBLE PRECISION (A - H , O - Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 2o) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B04/ NM, NSTEP, NFT, NJT.NST, NEQ, LHB, NSZ,NA1 COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFo,M0Z), XF(MFT,MEQ).FF(MFT,MEQ ),RFF(MFT,MEq ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ ),SVEC2(MFT.MEq ),XFX(MFT,MEQ ).EK(8.8) C THE SYSTEM MATRIX GSTIF IS STORED COLUMNWISE. KB = LHB 1 TEMP = GSTIF(IKO,1) TEMP DSQRT(TEMP) GSTIF(IK0,1) = TEMP DO 10 I = 2, LHB GSTIF(IKO,I) = GSTIF(IKO,I) / TEMP 10 CONTINUE DO so .1 = 2, NEQ Ji = .1 - 1 IJD = LHB J KB SUM =_- GSTIF(IKO,IJD) Ko = 1 IF (.1 .GT. LHB) K0 = J • KD DO 20 K = KO. J1 .1K = KB K .1 - KB TEMP = GSTIF(IKO,JK) SUM = SUM • TEMP -* 2 20 CONTINUE GSTIF(IKO,IJD) DSQRT(SUM) DO 50 I = 1, KB II = .1 4- I K0= 1 IF (II .GT. LED) KC) II - KB SUM = GSTIF(IKO,IJD + I) IF (I .EQ. KB) GO TO 40 DC) 30 K = KC), .11 JK = KB K J - KB IK = KB K II - KB TEMP = GSTIF(IK0,JK) SUM = SUM - GSTIF(IKoJK) TEMP 30 CONTINUE 40^GSTIF(IKO,IJD + I) = SUM / GSTIF(IKO,IJD) 50 CONTINUE 60 CONTINUE RETURN END c 130 Appendix D. PTB SOURCE CODE AND I/O FILE^ C SUBROUTINE SOLVE SOLVES THE SYSTEM OF EQUATIONS USING THE C DECOMPOSED MATRIX OBTAINED FROM SUBROUTINE DCOMP. C SUBROUTINE SOLVE (IKO) IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ. LHB, NSZ,NA1 COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ),FF(MFT,MEQ),RFF(MFT,MEQ ),SVEC(MFT,MEQ ), 3^SVEC1(MFT.MEQ ),SVEC2(MFT,MEQ ),XFX(MFT,MEQ ),EK(8,8) C FORWARD SUBSTITUTION. KD = LHB - 1 TEMP = GSTIF(IKO,1) XF(IKO,1) = XF(IKO,1) / TEMP DO 20 I = 2, NEQ II = I - 1 KO = 1 IF (I .GT. LHB) KO = I - KB SUM = XF(IKO,I) II = LHB I - KB DO 10 K = KO, Ii IK = KB K + I - KB TEMP = GSTIF(IKO,IK) SUM = SUM - (TEMP XF(IKO,K)) 10 CONTINUE XF(IKO,I) = SUM / GSTIF(IKO,II) 20 CONTINUE C BACKWARD SUBSTITUTION. Ni = NEQ - 1 LB LHB NEQ - KB TEMP = GSTIF(IKO,LB) XF(IKO,NEQ) = XF(IKO,NEQ) / TEMP DO 40 I = 1, Ni Ii = NEQ . I + 1 NI = NEQ I KO = NEQ IF (I .GT. KB) KO = NI + KB SUM = XF(IKO,NI) II = LHB * NI • KB DO 30 K = Il, KO IK = KB NI K - KB TEMP = GSTIF(IKO,IK) SUM = SUM - (TEMP XF(IKO,K)) 30 CONTINUE XF(IKO,NI) = SUM / GSTIF(IKO,II) 40 CONTINUE 131 Appendix D. PTB SOURCE CODE AND I/O FILE ^ RETURN END C C ITERATING METHOD IS USED IN THIS SUBROUTINE TO SOLVE C THE COUPLING EQUATIONS. C SUBROUTINE ITERATE IMPLICIT DOUBLE PRECISION (A - H, 0 -Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 2o) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT, NST.NEQ, LHB, NSZ.NA1 COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFo,MOZ), 1^XF(MFT,MEQ ),FF(MFT,ME Q ),RFF(MFT,MEQ ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ).SVEC2(MFT,MEQ),XFX(MFT,MEQ),EK(8,8) COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12), 1^DDX,DDZ.FX.FZ,TOLL.NNA DIMENSION X(MEQ),XX(MEQ) IF (NM .EQ. 1) GO TO 123o C^ 1090 NFLAG = 1 DO 1220 IK = 1, NM, NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM = (IK + 1) / 2 DO 1100 I = 1. NEQ 1100 XX(I) = FF(IKM,I) DO 1170 IN = 1, NM, NSTEP IF (IN .EQ. IK) GO TO 1170 INM = IN IF (NSTEP .EQ. 2) INM = (IN + 1) / 2 IF (IN .GT. IK) GO TO 1110 IKO = (INM 1) * NFT + IKM INM (INM + 1) / 2 GO TO 1120 1110 IKO = (IKM 1)* NFT + INM IKM (IKM + 1) / 2 C 1120^DO 100 IE = 1, NST IJ1 = ( IE - 1 ) 19 + 14 1.12 = ( IE - 1 ) 19 + 14 + 11 = ( IE - 1 ) * 144 DO Bo I = 1.11,IJ2 TEMP = 0.0 CHI = I - (IE - 1 ) *19 .13 DO 60 = DI, 1.12 CHII = J - (IE - 1 ) *19 -13 TEMPI = SVEC(INM,J) CHICHI = IJJ + (CHI-1)*12 + CHII 132 Appendix D. PTB SOURCE CODE AND I/O FILE ^ 60^TEMP = TEMP + ASTIF(IKO,CHICHI) * TEMPI 80^XX(I) XX(I) - TEMP 100^CONTINUE C 1170 CONTINUE C ^ DO 1190 I = 1, NEQ 1190 XF(IKM,I) = XX(I) CALL SOLVE(IKM) DO 1195 I = 1,NEQ 1195 X(I) = XF(IKM,I) C TEMP = 0.0 TEMPI. = 0.0 DO 1200 I = 1, NEQ TEMP = TEMP + (X(I)**2) TEMPI = TEMPI (SVEC(IKM,I)**2) 1200 CONTINUE TEMP = DSQRT(TEMP) TEMPI = DSQRT(TEMP1) EPS = DABS((TEMP TEMPI)/TEMPI) C GIKM = 0.0 AIKO = 0.0 DO 1204 NANA = 1,NSZ 1204 GIKM = GIKM (GSTIF(IKM,NANA))**2 DO 1206 MAMA = 1,NA1 1206 AIKO = AIKO (ASTIF(IKO,MAMA))**2 • WRITE(",'")'EPS = SOLVE',EPS,'IKM^',IKM • WRITE(*,*)'GIKM = ',GIKM,' IKM = ',IKM • WRITE(',*)'AIKO = ',AIKO,' IKO = ',IKO C^ NAC = 1 IF (EPS .GT. TOL) NAC = 0 NFLAG = NFLAG NAC DO 1210 I = 1, NEQ 1210 SVEC(IKM,I) = X(I) 1220 CONTINUE IF (NFLAG .EQ. 1) GO TO 1230 GO TO 1090 1230 CONTINUE RETURN END C *************************************************************** C THIS SUBROUTINE CHECKS THE DEFORMATIONS OF THE SPRING ELEMENTS C TO SEE WHICH DEFORMATION EXCEEDS THE LIMITION. C .31K7/4*}[Mfitae”..malc.*#,Icyzac...*.*.M.....leat#*71[..1.1.1.1...1■WW*.WMC...kalcaKaRVeatalCa.F. 133 Appendix D. PTB SOURCE CODE AND I/O FILE^ SUBROUTINE CHECK(LMNLMN) C^ LMNLMN 1 WRITE = 0 DON'T WRITE IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI. PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, PDV, PDG COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ),FF(MFT,MEQ ),RFF(MFT.MEQ ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ),SVEC2(MFT,MEQ),XFX(MPT,MEQ).EK(8,8) COMMON /B15/ EX(MSP),EY(MSP),EZ(MSP),NAN,NP,MNL COMMON /B17/ ICH(MSP,MNP ),ICH1(MSP,MNP ),SL(12). 1^DDX,DDZ,FX.FZ,TOLL,NNA COMMON /B18/ SIGMAY,WDX1(MSP),WDX2(MSP), 1^WDY1(MSP ),WDY2(MSP ),WDZ1(MSP ),WDZ2(MSP ) DIMENSION S(12) C DO 5 IK =1,NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM =(1.+IK)/2 DO 5 I = 1,NEQ 5 SVEC(IKM,I) = SVEC(IKM,I) + SVEC2(IKM,I) DO 170 IE = 1,NST WDX1(IE) = 0.D0 WDX2(IE) = 0.D0 WDY1(IE) = 0.D0 WDY2(IE) 0.D0 WDZ1(IE) = 0.D0 WDZ2(IE) = 0.D0 DO 150 IP = 1,NP XI = (IP - 1)`XL/(NP - 1) DX1 = 0.0 DX2 = 0.0 DY1 = 0.0 DY2 = 0.0 DZ1 = 0.0 DZ2 = 0.0 DO 130 IK = 1,NM,NSTEP IKO = IK IF (NSTEP .EQ. 2) IKO = (IK + 1)/2 DO 20 I =14.19 .1 = (IE 1)"19 +I 134 ^ 0 9 +^- .11) = ZZX +^- dI) = TZX +^- dI) = ZAX C^9.(T • dI) = TAX Z^- dI) = ZXX I ÷^- dI) = TXX 1RNIS,K0 011 ^ 0 zz + zza = zza tz + TZCI = TZCI CA z.ka = zAa TA. + TA.411 = TA.(11 CX zxa = zxa IX + TXCI = TXCE (L)S.('IX/IX.I.I.XI)NISCIA- I (T)S.(7X/IX.Id.XI)NISCI - =ZZ (4)s . (gx / ix.tdotthosa + t (6)9.(gx/tx.td.m)sooa+ (z ) s.(gxlix.H .xt )so 0 a.(qx/ta,ou). .z /xxa + (c)s.(gx/xx.H.m)so0a(T)s.(qx/tx.moll)so3a*hx/td.)a)..zhna—zx (6)s.(qx/tx.L.L0n)so0a+ e (2,)s.(qx/tx.m.,u)so0a.(rix/H.3u).•z/ma,d- z (e)s.(qx/tx.mou)sooa- T (I)S.(•IX/IX.Id.XI)S00C1*('IX/Id.XI).'Z/X1,d=TX C (9)s.(rix/tx-moit)msa.(uas*•z)/mt,d+ z (9)s.(qx/tx,aa*xt)tcsa(z)s.(qx/tx,ad.m)Nusa.(1,ras.'z)/xJ..1-=zA (it)s.eix/tx,da*mhsasa+ e (9)s.(qx/tx.idoll)tosa.(.1,rds..z)/xLa- z (9)9.(qx/tx.ra.3a)msa- I ^nNISNO0 Of (r`03a)oans = (x)s 9+I=X I + 61. =r 9'1 = I OF OCI HRNISNO0 OZ (r`OXI)OHAS = (X)S CT - I = X 9ET TIM (VI (INV acrop 3 Janos Ead •a xrpuaddv - . Appendix D. PTB SOURCE CODE AND I/O FILE ^ IF (DY1 .LE. 0.0 .AND. ICH(IE,KY1) .NE. 0) THEN IF (ICH(IE.KX1) .EQ. 1) THEN IF (DADS(DX1) .GE. DDX) THEN ICH(IE.KX1) = 0 IF(DX1 .GT. 0.0) ICH1(IE,KX1) = 1 IF(DX1 .LT. 0.0) ICH1(IE,KX1) = -1 NNA NNA 1 END IF END IF ICH(IE,KY1) = -1 END IF C IF (DY1 .GT. 0.0) THEN ICH(IE,KX1) = 0 ICH1(IE,KX1) = 0 NNA = NNA 1 ICH(IE,KY1) = 0 NNA = NNA 1 END IF C IF (DY2 .LE. 0.0 .AND. ICH(IE,KY2) .NE. 0) THEN IF (ICH(IE,KX2) .EQ. 1) THEN IF (DABS(DX2) .GE. DDX) THEN ICH(IE,KX2) = 0 IF(DX2 .GT. 0.0) ICH1(IE,KX2) = 1 IF(DX2 .LT. 0.0) ICH1(IE,KX2) = -1 NNA = NNA 1 END IF END IF ICH(IE,KY2) = -1 END IF C IF (DY2 .GT. 0.0) THEN ICH(IE,KX2) = 0 ICH1(IE,KX2) = 0 NNA NNA 1 ICH(IE,KY2) = 0 NNA NNA 1 END IF C C IN Z - DIRECTION " Z1 " IF (DY1 .LE. 0.0 .AND. ICH(IE,KY1) .NE. 0) THEN IF (ICH(IE,KZ1) .EQ. 1) THEN IF (DABS(DZ1) .GE. DDZ) THEN 136 Appendix D. PTB SOURCE CODE AND I/O FILE^ ICH(IE,KZ1) = 0 IF(DZ1 .GT. 0.0) ICH1(IE,KZ1) = 1 IF(DZ1 .LT. 0.0) ICH1(IE,KZ1) = -1 NNA NNA 1 END IF END IF END IF C IN Z - DIRECTION ** Z2 " IF (DY1 .LE. 0.0 .AND. ICH(IE,KY1) .NE. 0) THEN IF (ICH(IE,KZ2) .EQ. 1) THEN IF (DABS(DZ2) .GE. DDZ) THEN ICH(IE,KZ2)^0 IF(DZ2 .GT. 0.0) ICH1(IE,KZ2)^1 IF(DZ2 .LT. 0.0) ICH1(IE,KZ2) = -1 NNA = NNA 1 END IF END IF END IF IF (DY1 .GT. 0.0) THEN ICH(IE,KZ1) = 0 ICH1(IE,KZ1) 0 NNA = NNA 1 END IF IF (DY2 .GT. 0.0) THEN ICH(IE,KZ2) = 0 ICH1(IE,KZ2) = 0 NNA = NNA 1 END IF IF (DABS(DX1) .GE. WDX1(IE)) WDX1(IE) = DABS(DX1) IF (DABS(DX2) .GE. WDX2(IE)) WDX2(IE) = DABS(DX2) IF (DABS(DZ1) .GE. WDZ1(IE)) WDZ1(IE) = DABS(DZ1) IF (DABS(DZ2) .GE. WDZ2(IE)) WDZ2(IE) = DABS(DZ2) IF ((DY1) .LE. WDY1(IE)) WDY1(IE) = (DY1) IF ((DY2) .LE. WDY2(IE)) WDY2(IE) = (DY2) 150 CONTINUE 170 CONTINUE DO 200 IK =1,NM,NSTEP IKM = IK IF (NSTEF .EQ. 2) IKM =(1+IK)/2 DO 200 I = 1,NEQ 200 SVEC(IKM,I) = SVEC(IKM,I) SVEC2(IKM,I) RETURN END 137 ^ Appendix D. PTB SOURCE CODE AND I/O FILE^ c seg. ae.**.*.**214....*^ ale***xe3k..mx*M..*.. C THIS SUBROUTINE CALCULATES THE IMBALANCE LOAD CAUSED C BY THE NON LINEAR DEFORMATION OF THE SPRING ELEMENTS C SUBROUTINE IMBALANCE IMPLICIT DOUBLE PRECISION (A - H, 0 Z) PARAMETER (MJT = 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /BOG/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, lux COMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG. 1^PDX, PDY. PDV, PDG COMMON /B11/ IBC(MBC,2), NBC COMMON /B12/ ESTIF(19,19) COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), XF(MFT,MEQ),FF(MFT,MEQ ).RFF(MFT,MEQ ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ ),SVEC2(MFT,MEQ ),XFX(MFT,MEQ),EK(8,8) COMMON /B14/ XK(8,8),YK(8,8),ZK(8,8), 1^BX1(1,8),BX2(1,8),BY1(1,8),BY2(1,8). 2^BZ1(1,8),B22(1,8),BX1T(8.1),13X2T(8,1),BY1T(8,1), 3^BY2T(8,1),BZ1T(8,1),BZ2T(8,1),SXK(8,8), 4^SYK(8,8),SZK(8,8),EK12(12,12) COMMON /B15/ EX(MSP),EY(MSP),EZ(MSP),NAN,NP,MNL COMMON /B17/ ICH(MSP.MNP).ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA DIMENSION SG(19),SG1(19) CALL CHECK(0) DO 5 IK =1,NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM =(1-1-IK)/2 DO 5 I = 1,NEQ 5 SVEC(IKM,I) = SVEC(IKM,I) SVEC2(IKM,I) DO 20 IK = 1,NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM = (IK + 1)/2 DO 10 I = 1,NEQ 10 XF(IKM,I) XFX(IKM,I) DO 20 IE = 1,NST CALL SLOAD(IE,IK) IJ = (IE^1)=19 -I- 13 DO 20 I = 1,12 20^XF(IKM,IJ-f-I) =^ SL(I) C DO 300 IK = 1,NM,NSTEP IKM = IK 138 Appendix D. PTB SOURCE CODE AND I/O FILE^ IF (NSTEP .EQ. 2) THEN IKM = (IK+1)12 END IF DO 300 IN = 1.NM,NSTEP IF (IN -EQ. IK) THEN FAC2 (IK"2) * PI2 (2.0D0 * XL) FAC4 = (IK"4) * PI4 / (2.0D0 (XL**3)) DO 200 IE = 1,NJT CALL COVCON(IK,IN) ESTIF(10,10) = ESTIF(10,10) + (EJT(IE) RIY FAC4) ESTIF(11,11) ESTIF(11,11) + (EJT(IE) AJT FAC2) ESTIF(12,12) = ESTIF(12,12) + (EJT(IE) *RIZ FAC4) ESTIF(13,13) = ESTIF(13,13) + (GJT(IE) RIT FAC2 / 1^ (SPJT"2)) Il = (IE •1)*19 +1 12 -= IX 19 J=0 DO 120 I = 11,12 J = 3+1 120 SG(J) = SVEC(IKM,I) CALL DGMATV (ESTIF,SG,SG1,19,19,19) =0 DC) 140 I = 11.12 J=3+1 140 XF(IKM,I) = XF(IKM,I) SG1(J) IF (IE .LE. NST) THEN CALL STIF12(IE,IK,IN) K=0 = (IE 1)*19 + 14 12 = Ii +11 DC) 160 I = 11,12 K = K+1 160 SG(K) = SVEC(IKM,I) CALL DGMATV (EK12,SG,SG1,12,12,12) K=0 DO 180 I = 11,12 K K+1 180 XF(IKM,I) = XF(IKM,I) SG1(K) END IF 200 CONTINUE END IF IF (IN .NE.IK) THEN INM = IN IF (NSTEP .EQ. 2) INM = (IN +1)/2 DO 290 IE = 1,NJT IF (IE .LE. NST) THEN CALL STIF12(IE,IK,IN) 139 Appendix D. PTB SOURCE CODE AND I/O FILE^ K=0 = (IE 1)*19 + 14 12 =^4-11 DO 260 I = 11,12 K = K+1 260 SG(K) = SVEC(IKM,I) CALL DGMATV (EK12,SG,SG1,12,12,12) K=0 DO 280 I = 11,12 K = K+1 280 XF(INM,I) = XF(INM,I) - SG1(K) END IF 290 CONTINUE END IF 300 CONTINUE C APPLIES THE BOUNDARY CONDITIONS TO THE LOAD VECTOR XF(IKO,M) IF (NBC.NE.0) THEN DO 330 I = 1, NBC NE = IBC(I,i) IDOF IBC(I,2) M = (NE - 1) * 19 4- IDOF C DO 330 IK = 1.NM,NSTEP IKM = IK IF (NSTEP .EQ. 2) IKM = (IK+1)/2 XF(IKM,M) = 0.0 330 CONTINUE END IF DO 365 IK =1,NM,NSTEP 1KM = IK IF (NSTEP .EQ. 2) IKM =(14-IK)/2 DC) 365 I = 1,NEQ 365 SVEC(IKM,I) = SVEC(IKM,I) SVEC2(IKM,I) RETURN END c C THIS SUBROUTINE IS TO MULTIPLY A RECTANGULAR M BY N C MATRIX BY A VECTOR OF LENGTH N . C SUBROUTINE DGMATV(A, V, W, M, N, NDIMA) IMPLICIT DOUBLE PRECISION (A - H, 0 Z) REAL'S A(NDIMA,1), W(1), V(1) DC) 10 I = 1, M W(I) = 0.D0 DC) 10 J = 1, N 10 W(I) = W(I)^A(I,J) V(J) RETURN 140 ^ Appendix D. PTB SOURCE CODE AND I/O FILE ^ END C C SUBROUTINE RESULT EVALUATES THE DEFLECTIONS AND STRESSES IN THE C JOIST AND PLATE COVER OF THE T-BEAM STRIP. C MC* MC W.....*...4.*...**D.1(.7..2***711.....‘,1,..*.W.11t MeaR W.M..31[4411,2411e.*.YItiM.AW SUBROUTINE RESULT IMPLICIT DOUBLE PRECISION (A - H,^Z) PARAMETER (MJT 10, MEQ = 190, MSZ = 3610, MOZ = 1296) PARAMETER (MFT = 4, MLD = 12, MBC = 20, MPF = 20) PARAMETER (MSP = 9, MFO = 6, MNP = 378) COMMON /B00/ PI, PI2, PI4, TOL COMMON /B03/ ETA(6), GWT(6) COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1 COMMON /B05/ XL, SPJT COMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RIT COMMON /B07/ PTK, PEX, FEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG, 1^PDX, PDY, FDV, PDG COMMON /B08/ CKPAL, CKPER, CKROT, SPC,XIN COMMON /B10/ XPTF(MPF), PTF(MPF),NPTF COMMON /B13/ GSTIF(MFT,MSZ), ASTIF(MFO,MOZ), 1^XF(MFT,MEQ ),FF(MFT,MEQ),RFF(MFT,MEQ),SVEC(MFT,MEQ ), 3^SVEC1(MFT,MEQ ),SVEC2(MFT,MEQ ),XFX(MFT,MEQ ),EK(8,8) COMMON /B15/ EX(MSP),EY(MSP),EZ(MSP),NAN,NP,MNL COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12), 1^DDX,DDZ,FX,FZ,TOLL,NNA COMMON /B18/ SIGMAY,WDX1(MSP ),WDX2(MSP), 1^WDY1(MSP),WDY2(MSP),WDZ1(MSP ),WDZ2(MSP) DIMENSION F(MEQ), WB(19), WS(19), ST(19), WP(38), S1T(19), 1^S1M(19),S1B(19),S2T(19),S2M(19),S2B(19) C WRITE(3,*) WRITE(3,1' WRITE(3,1' SUMMARY OUTPUT' WRITE(3,1' WRITE(3,1 WRITE(3,*)'WDY IS MAX. RELATIVE DISPL. IN Y-DIRECTION' WRITE(3,1'WDZ IS MAX. RELATIVE DISPL. IN Z-DIRECTION' WRITE(3,1 WRITE(3,1161)PTF(1) 1161 FORMAT(1X,' POST TENSION FORCE = E12.5) WRITE(3,*) WRITE(31' SIGMAY = ',SIGMAY WRITE(3,2) 2 FORMAT(1X ,'IE ',4X,'SIGMAY',9X ,'WDY1',9X ,'WDY2'.9X ,'WD Z1 ',11X , 1^'WDZ2') WRITE(3,*) DO 3 IE =1,NST 3^WRITE(3,4)IE,SIGMAY,WDY1(IE),WDY2(IE),WDZ1(IE),WDZ2(IE) 4^FORMAT(1 X ,I3 ,3X ,E10.4,3X ,E10 .4,3X ,E10.4,3X ,E10.4 ,3X ,E 10.4) 141 Appendix D. PTB SOURCE CODE AND I/O FILE^ W RITE(3,*) C IF (NAN .EQ.1) THEN WRITE(2,1 END IF IF (NAN .EQ. 0) THEN WRITE(2,*) WRITE(2,1' ELASTIC LINEAR ANALYSIS END IF W.TMAX = 0.0D0 SJMAX = 0.0D0 WPMAX = 0.0D0 SPMAX = 0.0D0 ALPHA = 1.2D0 C LOOP OVER EACH T-BEAM STRIP DO 170 IE = 1, NJT C C EVALUATING DEFLECTIONS AND STRESSES IN THE JOISTS. DC) 10 I = 1, 19 WH(I) = 0.0D0 WS(I) = 0.0D0 ST(I) 0.0D0 10 CONTINUE DO 90 IK = 1, NM, NSTEP PIN = PI * IK PINL = PIN / XL IKO IK IF (NSTEP .EQ. 2) IKO = (IK + 1) / 2 DO 20 I = 1, 19 J= ((IE - 1) * 19) + I F(I) = SVEC(IKO,J) 20 CONTINUE C FACS = EJT(IE) * ((HJT * F(10) * (PINL**2) / 2.0D0) 1^(F(11) PINL)) FAC1 = CKPAL / (2.0 * GJT(IE) * PINL BJT * SPC) FACW = EJT(IE) ((PINL HJT)**2) / (12.0D0 * GJT(IE)) FACW = FACW + (CKPAL (HJT + PTK)) / 1^(4.0D0 GJT(IE) BJT SPC) FACW = (FACW * F(10)) - (FAC1 * (F(7) - F(11))) FACW = ALPHA * FACW DO 30 I = 1, 19 XIL DSIN(I*PIN/20) WB(I) = WB(I) + (F(10) * XIL) WS(I) = WS(I) + (FACW XIL) ST(I) ST(I) + (FACS * XIL) 30 CONTINUE 90 CONTINUE C 142 Appendix D. PTB SOURCE CODE AND I/O FILE^ WJD = 0.0D0 WJS 0.0D0 WJT = 0.ODO ST.1 = 0.0130 DO 50 I = 1, 19 IF (DABS(WB(I)) .GE . DABS(WJB)) WJB = WB(I) IF (DABS(WS(I)) .GE. DABS(WJS)) WJS = WS(I) WT = WB(I) WS(I) IF (DABS(WT) .GE. DADS(WJT)) WJT WT IF (DABS(ST(I)) .GE. DABS(STJ)) STJ = ST(I) 50 CONTINUE C EVALUATING DEFLECTIONS AND STRESSES IN THE PLATE COVER. DO 60 I = 1, is S1T(I) = 0.0D0 S1M(I) = 0.0D0 C S1B(I)^0.0D0 S2T(I) = 0.0D0 S2M(I) = 0.ODO S2B(I) = 0.0D0 C 60 CONTINUE DO 62 I = 1,38 62 WP (I) = 0.ODO C EPC = PEY / (1.0D0 (PVXY * PVYX)) DO 130 IK = 1, NM, NSTEP PIN = PIIK PINL = PIN / XL IKO = IK IF (NSTEP .EQ. 2) IKO = (IK 1) / 2 DO 70 I = 1, 19 = ((IE - 1) 19) + I F(I) = SVEC(IKO,J) 70 CONTINUE C HT = PTK / 2.0D0 FAC1T = (F(6) / SPJT) (PVYX PINL * F(3)) 1^((-HT) ((-46.0D0 * F(1) )^(14.M° * F(14)) + 2^( 32.0D0 * F(10)) - (12.0D0 " F(2) ) 3^( 2.0D0 * F(15)) (16.0D0^F(9) )) / 4^(SPJT**2))^(PVYX * (-HT) F(1) * (PINL**2)) FAC1M = (F(6) / SPJT) (PVYX PINL * F(3)) 1^((0.0) * ((-46.0D0^F(1) )^(14.0D0 * F(14)) 2^( 32.0D0 F(10)) - (12.0D0 - F(2) ) 3^( 2.0D0 * F(15)) - (16.0130 * F(9) )) / 4^(SPJT”2))^(PVYX * (0.0) *F(1) * (PINL**2)) 143 ^ aana,Lism^oft ^ 0 (mx.(Wa) + (I)aM = (I)dM^bit (OZ/NId*(6T - I))NIS = 'IIX 41Z 4 OZ = I iZT Oa ('IIX. (Oa) + (I)dM = (I)dM (OZ/IslId*I)NIS = 'IIX 61'T = I ZZT OQ ^ 0 ariNumoo ('lix azona) + (i)azs = (Dais (ant vgzova) + (I)1lis = (DiAzzs ('IIX * Izova) + (i),Lzs = Wizs ('IIX amva) + Wats = (Dais vgmva) + (1)1-kas = Winas (gix^lova) + (i)Its = (I)Sts (oz/ma-Dmsa = 6T^= I OZT OCT Oda * aZOIrel = aZ0V3 Oda * L1tZOV3 = TAIZOV3 Oda I.Z0V3 = IZOVeI Oda * EITOArd = adGI * TAIT01,4 = mova oda^= STOV3 ((z..maa) - (004 - (Ix-) * - (((ta)a • (z)a + ((ot)a^ocro•90 -^- Kra) + ( (i)d^oao-8)) ((z..a.rds) / (Ix-))) + ((z)d^* XAAd) • (u.'s / (((g t)d^oasz•0) - ((Oa * (Ja•()) (0103 OQO•T) + (Wel * OCITT - ))) = azoird - ((z—rund) (ot)a (ro) XAAd) - (((g^- (z)a + ((ot)3^0Q0 '9T) - ((bad oac•g) + ( (Oa * °a•e)) ((z--ad ds) / (o•0))) + ((c)a^xAna) • (LI dS M6I)eI* OCI9•0) - ( MA 0a9Z • 0) - ((gT)a °a•l) + ((g)a^oas•t-))) = riqz0v3 - ((Z.WINId)* (OatI * 1,H XAAd) - (((7T)3 - (1)3 + ((oT)a * 0a0 • 9T) oao•g)) (( " )3 oac•e) + ( (Oa^ ((z—zras) / Ix)) (Wet^XAAd) - (Irds / (((603 * oagz•o) - ((Oa oagz•o) - ((Wei oaog•T) + ((g)a oag•T-))) = Izova 0 ((Z**rINId) (OA * (Lx) * XAAd)^((Z**.LfdS) (( (6)4 oacrgi) - ((Wet - (Ku•z ) - ( (z)a^oao•zt) - ((ot)a^oao•ze ) + ((g^oao•gt) + ( (Oa * oao•gs-)) (Ix)) ((1)3*^XAAd) (If<IS / Ma) = aTOV3 ttT T11,1 o/I GNV acroo aounos ELM 'a xnDua ddv Appendix D. PTB SOURCE CODE AND I/O FILE ^ WPC = 0.0D0 Si = 0.0D0 S2 = 0.0D0 S3 = 0.0D0 C S4 = 0.0D0 S5 = 0.0D0 S6 = 0.0D0 C DC) 140 I = 1, 19 IF (DABS(S1T(I)) .GT. DABS(S1)) S1 = DABS(S1T(I)) IF (DABS(S1M(I)) .GT. DABS(92)) S2 = DABS(S1M(I)) IF (DABS(S1B(I)) .GT. DABS(S3)) S3 = DABS(S1B(I)) IF (DABS(52T(I)) .GT. DADS(S4)) S4 = DABS(S2T(I)) IF (DABS(S2M(I)) .GT. DABS(S5)) S5 = DABS(S2M(I)) IF (DABS(S2B(I)) .GT. DABS(S6)) S6 = DABS(S2B(I)) 140 CONTINUE DO 144 I = 1,38 144 IF (DABS(WP(I)) .GT. DABS(WPC)) WPC = WP(I) C C PRINTING OUT RESULTS. C WRITE (2,155) IE 155 FORMAT (/, ' *" ELEMENT ', I2, ' **') WRITE (2,160) WJB, WJS, WJT. STJ, WPC, S1,S3,54,S6 160 FORMAT (' MAX. JOIST DEFLECTION (BENDING) = E12.5, /, 1^' MAX. JOIST DEFLECTION ( SHEAR) = E12.5, /, 2^' MAX. JOIST DEFLECTION ( TOTAL) ^E12.5, I, 3^' MAX. JOIST BENDING STRESS^= ',E12.5, /, 4^' MAX. COVER DEFLECTION ^= E12.5, /, 5^' MAX. COVER STRESS (NODE 1 TOP) = E12.5, /. 7^' MAX. COVER STRESS (NODE 1 BOTTOM)= E12.5, /, 8^' MAX. COVER STRESS (NODE 2 TOP) = E12.5, /, 9^' MAX. COVER STRESS (NODE 2 BOTTOM)= E12.5 ) C WRITE (3,155) IE WRITE (3,161) WJT, STJ, WPC 161 FORMAT (' MAX. JOIST DEFLECTION ( TOTAL) = E12.5, /, 3^' MAX. JOIST BENDING STRESS^= E12.5, /, 4^' MAX. COVER DEFLECTION^= E12.5, ) C IF (DABS(STJ) .GE. DABS(SJMAX)) SJMAX = STJ IF (DABS(WJT) .GE. DABS(WJMAX)) WJMAX = WJT IF (DABS(WPC) .GE. DABS(WPMAX)) WPMAX = WPC IF (DABS(S1) .GE. DABS(SPMAX)) SPMAX = 51 IF (DABS(S3) .GE. DABS(SPMAX)) SPMAX = S3 IF (DABS(54) .GE. DABS(SPMAX)) SPMAX = S4 IF (DABS(S6) .GE. DABS(SPMAX)) SPMAX = S6 145 Appendix D. PTB SOURCE CODE AND I/O FILE ^ 170 CONTINUE C SUMMARY OF THE ANALYSIS. WRITE (2,180) WJMAX, SJMAX, WPMAX, SPMAX 180 FORMAT (///, ix, 12('"'), ' SUMMARY OF FLOOR ANALYSIS ', 11(n, 1^//. ' MAX. JOIST DEFLECTION ^= E12.5, /, 2^' MAX. JOIST STRESS^= E12.5, /, 3^' MAX. COVER DEFLECTION ^E12.5, /, 4^' MAX. COVER STRESS BETWEEN JOISTS = E12.5 WRITE(2,90) 90 FORMAT (/,1X, 50(''')) RETURN END END 146 Appendix D. PTB SOURCE CODE AND I/O FILE ^ D.2 PTB.DAT "`' PTB^PTB^PTB " PTB 7 3 1^ 12000.0 1000.0^ 100.0 500.0^ FOURIER NO. JOIST SYMMETRIC XL, SPJT BJT,HJT 14000.00 19000.00 14000.00 17.0^ 100.0 .5 .65^ REG PTK, EMUX,EMUZ 19000.00 1400.00 1000.00 0.02 0.20^PEX,PEY,PG,PVXY,PVYX 10.0 0.0 1.0E06 1.0E06 1.0E06^SPC,XIN,CKPAL,CKPER,CKROT 3.0,2.0^ TAOSX,TAOSZ 1 21 1^ NAN, NP, MNL 16.9E+01 16.9E-4-01^ EX(I) 1.00E+10 1.00E+10^ EY(I) 17.4E+01 17.4E-f-01^ EZ(I) 4 1 4500.00 5000.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 15 1 7000.00 7500.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 15 2 4500.00 5000.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 15 2 7000.00 7500.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 15 10 2.0E5 78.53 10.0 2.0 600.0 1.00e+05 1800.0 1.00e+05 3000.0 1.00e+05 4200.0 1.00e+05 5400.0 1.00e+05 6600.0 1.00e+05 7800.0 1.00e-1-05 9000.0 1.00e+05 10200.0 1.00e+05 11400.0 1.00e+05 0 ES AS DIP RFAC 147 Appendix D. PTB SOURCE CODE AND I/O FILE^ D.3 PTB.OUT ""*****"*"**"*""*""" FLOOR ANALYSIS PROGRAM **"*****""**""*"***"" PROBLEM TITLE: " PTB PTB " PTB^PTB NUMBER OF FLOOR JOISTS^= 3 ^ NUMBER OF FOURIER TERMS USED 4 MAX. ORDER OF FOURIER TERM^--:-- 7 PROPERTIES AND DIMENSIONS OF JOISTS JOIST SPAN^= 0.12000E+05 JOIST SPACING = 0.10000E+04 JOIST WIDTH^= 0.10000E+03 JOIST DEPTH^= 0.50000E +03 JOIST NO. = 1 EJT = 0.19000E+05 GJT = 0.82353E+03 JOIST NO. = 2 EJT = 0.14000E+05 GJT = 0.82353E+03 JOIST NO. = 3 EJT = 0.14000E+05 GJT = 0.82353E+03 PROPERTIES AND DIMENSIONS OF PLATE COVER COVER THICKNESS = 0.10000E+03 KX = 0.11714E+10 KY = 0.11714E+09 KV = 0.23427E+08 KG = 0.83333E+08 DX = 0.14056E+07 DY = 0.14056E+06 DV = 0.28112E+05 DG ^0.10000E+06 PROPERTIES FOR CONNECTORS STIFFNESS PARALLEL TO JOIST ^= 0.10000E+07 STIFFNESS PERPENDICULAR TO JOIST = 0.10000E+07 ROTATIONAL STIFFNESS FLANGE/JOIST = 0.10000E+07 SPACING BETWEEN CONNECTORS = 0.10000E+02 APPLIED TRANSVERSE LOADING JOIST^X1^X2^Y1^Y2^LOAD 1^0.95000E+04 0.50000E+04 0.20000E+03 0.50000E+03 0.10000E+00 1^0.70000E+04 0.75000E+04 0.20000E+03 0.50000E+03 0.10000E+00 2^0.45000E+09 0.50000E+04 0.20000E+03 0.50000E+03 0.10000E+00 2^0.70000E+04 0.75000E+04 0.20000E+03 0.50000E+03 0.10000E+00 POST TENSIONING FORCES X-LOC^FORCE 0.60000E+03 0.10000E+06 0.18000E+04 0.10000E+06 0.30000E+04 0.10000E+06 0.42000E+04 0.10000E+06 148 Appendix D. PTB SOURCE CODE AND I/O FILE ^ 0.54000E+04 0.10000E+06 0.66000E+04 0.10000E+06 0.78000E+04 0.10000E+06 0.90000E+04 0.10000E+06 0.10200E+05 0.10000E+06 0.11400E+05 0.10000E+06 NO. OF BOUNDARY CONDITIONS = 0 THREE DIMENSIONAL SPRING BETWEEN TWO T- BEAMS 42 STIFF. OF THE SPRING IN X-DIC. JOINT NO. = 1 EX = 0.86224E+05 JOINT NO. = 2 EX = 0.86224E+05 STIFF. OF THE SPRING IN Y-DIC. JOINT NO. = 1 EY = 0.10000E+11 JOINT NO. = 2 EY = 0.10000E+11 STIFF. OF THE SPRING IN Z-DIC. JOINT NO. = 1 EZ = 0.88776E+05 JOINT NO.. 2 EZ = 0.88776E+05 FRICTION COEFFICIENT OF THE COVER (X-DIC.) = 0.50000E+00 FRICTION COEFFICIENT OF THE COVER (Z-DIC.) = 0.65000E+00 STIFFNESS OF THE POSTENSIONING CABLE = 0.20000E+06 THE AREA OF THE POSTENSIONING CABLE = 0.78530E+02 THE DISTANCE BETWEEN THE CABLE TEETH = 0.10000E+02 THE RELAX FACTOR OF THE CABLE = 0.20000E+01 THE TURN NC). REQUIRED = 0.41773E+01 ELEMENT 1 "* MAX. JOIST DEFLECTION (BENDING) = 0.12558E+02 MAX. JOIST DEFLECTION ( SHEAR ) = 0.12855E+01 MAX. JOIST DEFLECTION ( TOTAL ) = 0.13843E+02 ^ MAX. JOIST BENDING STRESS 0.56304E+01 ^ 0.13665E+02 MAX. COVER DEFLECTION MAX. COVER STRESS (NODE 1 TOP) = 0.10653E+01 MAX. COVER STRESS (NODE 1 BOTTOM)= 0.99706E+00 MAX. COVER STRESS (NODE 2 TOP) = 0.96243E+00 MAX. COVER STRESS (NODE 2 BOTTOM)= 0.87784E+00 ELEMENT 2 ** MAX. JOIST DEFLECTION (BENDING) = 0.13133E+02 MAX. JOIST DEFLECTION ( SHEAR ) = 0.15197E+01 MAX. JOIST DEFLECTION ( TOTAL ) = 0.14652E+02 ^ MAX. JOIST BENDING STRESS = 0.64713E+01 ^ MAX. COVER DEFLECTION = 0.13660E+02 MAX. COVER STRESS (NODE 1 TOP) = 0.17037E+01 MAX. COVER STRESS (NODE 1 BOTTOM). 0.80633E+00 MAX. COVER STRESS (NODE 2 TOP) = 0.11909E+01 MAX. COVER STRESS (NODE 2 BOTTOM)= 0.69862E+00 149 Appendix D. PTB SOURCE CODE AND I/O FILE^ "" ELEMENT 3 "" MAX. JOIST DEFLECTION (BENDING) = 0.80606E+01 MAX. JOIST DEFLECTION ( SHEAR ) = 0.81426E+00 MAX. JOIST DEFLECTION ( TOTAL ) = 0.88749E+01 MAX. JOIST BENDING STRESS^= 0.35843E+01 MAX. COVER DEFLECTION^= 0.11385E+02 MAX. COVER STRESS (NODE 1 TOP) = 0.15598E+01 MAX. COVER STRESS (NODE 1 BOTTOM)= 0.77819E+00 MAX. COVER STRESS (NODE 2 TOP) = 0.94627E+00 MAX. COVER STRESS (NODE 2 BOTTOM)_ 0.91086E+00 "*".."*"*"*" SUMMARY OF FLOOR ANALYSIS **"'"*"'"`" ^ MAX. JOIST DEFLECTION = 0.14652E+02 ^ MAX. JOIST STRESS = 0.64713E+01 ^ MAX. COVER DEFLECTION = 0.13665E+02 MAX. COVER STRESS BETWEEN JOISTS = 0.17037E+01 150 Appendix D. PTB SOURCE CODE AND I/O FILE^ D.4 PTBB.OUT *****"*"`"*"*"*"*"" FLOOR ANALYSIS PROGRAM """*""*"*`*****""" PROBLEM TITLE: " PTB^PTB "" PTB^PTB NUMBER OF FLOOR JOISTS^= 3 NUMBER OF FOURIER TERMS USED = 4 MAX. ORDER OF FOURIER TERM = 7 SUMMARY OUTPUT WDY IS MAX. RELATIVE DISPL. IN Y-DIRECTION WDZ IS MAX. RELATIVE DISPL. IN Z-DIRECTION POST TENSION FORCE = 0.10000E+06 SIGMAY = 0.83333333333333 IE SIGMAY^WDY1^WDY2^WDZ1^WDZ2 1^0.8333E+00 -.2473E-05 -.3233E-05 0.8878E-02 0.8878E-02 2^0.8333E+00 -.2487E-05 -.3222E-05 0.1018E-01 0.1018E-01 " ELEMENT 1 " MAX. JOIST DEFLECTION ( TOTAL) =^0.13843E+02 MAX. JOIST BENDING STRESS = 0.56304E+01 MAX. COVER DEFLECTION = 0.13665E+02 '" ELEMENT 2 " MAX. JOIST DEFLECTION ( TOTAL) =^0.14652E+02 MAX. JOIST BENDING STRESS = 0.64713E+01 MAX. COVER DEFLECTION = 0.13660E+02 ** ELEMENT 3 "' MAX. JOIST DEFLECTION ( TOTAL) =^0.88749E+01 MAX. JOIST BENDING STRESS = 0.35843E+01 MAX. COVER DEFLECTION = 0.11385E+02 151
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Non-linear analysis for transversely post-tensioned...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Non-linear analysis for transversely post-tensioned timber bridge Cai, Shen 1992
pdf
Page Metadata
Item Metadata
Title | Non-linear analysis for transversely post-tensioned timber bridge |
Creator |
Cai, Shen |
Date Issued | 1992 |
Description | Timber bridges have been very important in North America due to an abundant resource of the material and the relatively unsophisticated requirement for equipment and skilled labour. The transversely post-tensioned laminating method can improve the loading distribution capacity of the timber bridge structure and elongate the bridge service life. To develop a method for its analysis it is necessary to consider the non-linear behaviour of between-beams friction in the bridge structure. A method and corresponding computer program (PTB) has been developed for the non-linear analysis for the transversely post-tensioned timber bridge. The finite strip element method is used in the analysis. Wheel loadings are idealized as patches with static uniformly distributed loadings. The beam to beam friction parameters were obtained from tests. The relative movements between beams, in three-directions, stresses and deformations of the structure are obtained by the PTB program. As an application, abridge with post-tensioned T-beams made of Parallam is considered. |
Extent | 4934885 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-12-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050474 |
URI | http://hdl.handle.net/2429/3258 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1992_spring_cai_shen.pdf [ 4.71MB ]
- Metadata
- JSON: 831-1.0050474.json
- JSON-LD: 831-1.0050474-ld.json
- RDF/XML (Pretty): 831-1.0050474-rdf.xml
- RDF/JSON: 831-1.0050474-rdf.json
- Turtle: 831-1.0050474-turtle.txt
- N-Triples: 831-1.0050474-rdf-ntriples.txt
- Original Record: 831-1.0050474-source.json
- Full Text
- 831-1.0050474-fulltext.txt
- Citation
- 831-1.0050474.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050474/manifest