UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Non-linear analysis for transversely post-tensioned timber bridge Cai, Shen 1992

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

Item Metadata

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

NON-LINEAR ANALYSIS FOR TRANSVERSELY POST- TENSIONEDTIMBER BRIDGEByCai,ShenB.Eng. & M.Eng. Beijing Institute of Civil Engineering and ArchitectureA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESCIVIL ENGINEERINGWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1991© Cai, Shen, 1992In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.(Signature)Department of  Civil EngineeringThe University of British ColumbiaVancouver, CanadaDate ^DE-6 (2/88)ABSTRACTTimber bridges have been very important in North America due to an abundant resourceof the material and the relatively unsophisticated requirement for equipment and skilledlabour. The transversely post-tensioned laminating method can improve the loadingdistribution 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 behaviourof between-beams friction in the bridge structure.A method and corresponding computer program (PTB) has been developed for thenon-linear analysis for the transversely post-tensioned timber bridge. The finite stripelement method is used in the analysis. Wheel loadings are idealized as patches with staticuniformly distributed loadings. The beam to beam friction parameters were obtainedfrom tests. The relative movements between beams, in three-directions, stresses anddeformations of the structure are obtained by the PTB program. As an application, abridge with post-tensioned T-beams made of Parallam is considered.iiTable of ContentsABSTRACT^ iiList of Tables viList of Figures^ viiAcknowledgement1 INTRODUCTION 11.1 Background ^ 11.2 Literature Review 31.3 Objectives of the Thesis ^ 42 NON-LINEAR FINITE ELEMENT ANALYSIS 62.1 Structural Model ^ 62.1.1^Middle surface displacement functions for the flange of the T-beam 62.1.2^Displacement functions for the centroid of the T-beam web 92.1.3^Degrees of freedom for the joint element ^ 112.2 Strain Energy of the Structural System and Virtual Work ^ 132.2.1^Strain energy in the flange of the T-beam element 142.2.2^Strain energy in the web^ 182.2.3^Strain energy in the web-flange connector ^ 202.2.4^Virtual work done by non-conservative forces in the joint element 232.3 Global Equation ^ 26iii2.3.1^Stiffness contribution from one T-beam element  ^272.3.2^Stiffness contribution from one joint element  ^292.3.3^Consistent load vectors  ^302.3.4 System of equations corresponding to one T-beam (Element Matrix) 322.4 Global Equations and Solution Method ^  332.4.1^Characteristics of the global equation  ^332.4.2 The incremental procedure ^  342.4.3 Newton-Raphson method  352.4.4 Jacobi iteration method  ^393 FRICTION PARAMETERS FOR THE SPRING MODEL^403.1 The Friction Test  ^403.2 Spring Model ^  43ANALYTICAL RESULTS^ 494.1 Beam Test ^  494.2 T-beam Analysis 524.2.1^Analysis for post-tensioning force only  ^544.2.2^Analysis for the one patch of vertical loading  ^554.2.3^Wheel loading analysis  ^615 CONCLUSIONS^73Bibliography^ 74Appendices^ 76A SHAPE FUNCTIONS AND DERIVATIVES^76ivB THE PTB PROGRAM^ 79B.1 PTB Program Specifications ^  79B.2 Program Structure ^  80C PTB USER'S MANUAL 85D PTB SOURCE CODE AND I/O FILE^ 92D.1 PTB.FOR ^  93D.2 PTB.DAT  147D.3 PTB.OUT ^  148D.4 PTBB.OUT  151List of Tables3.1 Data of Friction Test Samples ^ 423.2 Statistical Results of the Friction Test ^ 434.1 Properties of the Cantilever Beam 504.2 Status of the Springs in the Cantilever Beams ^ 504.3 Properties of T-beam Structure ^ 534.4 Analytical Results: Post-tensioning Force Only ^ 55viList of Figures1.1 Analytical Structural Model ^ 52.1 The Structural Model 72.2 The T-beam Model and The Joint Model ^ 72.3 The Distribution of the Discrete Springs 122.4 The Constitutive Relationship of the Springs ^ 122.5 The Vechicle and Post-tensioning Load 302.6 The Global Equation System^ 342.7 The Incremental Method 362.8 Newton-Raphson Method ^ 363.1 The Friction Test Setup 413.2 Two Kinds of Friction Test ^ 413.3 Friction Test Recording Curves(1) 443.4 Friction Test Recording Curves(2) ^ 453.5 Friction Test Recording Curves(3) 463.6 Idealized Friction Test Curves ^ 473.7 Same Slip Assumption 484.1 Two Cantilever Beam Structure ^ 504.2 Cantilever Beam Structure Displacements ^ 514.3 Built-in Solid Beams ^ 514.4 Three T-beam Structure 52vii4.5 Three T-beam Structure with Post-Tensioning Force Only ^ 544.6 Displacements of T-beam Structure under Post-Tensioning Force Only^564.7 T-beam Structure One Patch Loading Case^  574.8 Bending Stress of T-beam Structure One Patch Loading Case ^ 584.9 Displacements of T-beam Structure One Patch Loading Case ^ 584.10 Relative Displ. of T-beam Structure One Patch Loading Case ^ 594.11 Upward Deflection of the T-beam ^  604.12 Displacements of T-beam under One Patch Loading ^ 604.13 T-beam Structure CASE 1 ^  624.14 Bending Stress in T-beam Structure CASE 1 ^  634.15 Displacements in T-beam Structure CASE 1  634.16 Relative Displacements in T-beam Structure CASE 1 ^ 644.17 Bending Stress under Moving Loading ^  644.18 Displacements under Moving Loading  654.19 T-beam Structure CASE 2 ^  654.20 Bending Stress in T-beam Structure CASE 2 ^  664.21 Displacements in T-beam Structure CASE 2  664.22 Relative Displacements in T-beam Structure CASE 2 ^ 674.23 The Effect of the Stiffness of the Spring  ^684.24 T-beam Structure CASE 3 ^  694.25 Bending Stress in T-beam Structure CASE 3 ^  694.26 Displacements in T-beam Structure CASE 3  704.27 Relative Displacements in T-beam Structure CASE 3 ^ 704.28 Effect of the Fourier Terms on Bending Stress  ^714.29 Effect of the Fourier Terms on Displacement ^  714.30 Effect of the Fourier Terms on Relative Displacement in Joint 1 ^ 72viii4.31 Effect of the Fourier Terms on Relative Displacement in Joint 2 ^ 72B.1 PTB Program Flow Chart ^  81ixAcknowledgementThe author takes this opportunity to gratefully acknowledge the guidance of his advisorDr. Ricardo 0. Foschi throughout this research and thesis preparation. His suggestionsand advice during our discussions proved invaluable. The author also expresses his sin-cere gratitude to Dr. D.L. Anderson and Dr. H.G.L. Prion for reading this thesis andsuggesting improvements.The financial support in the form of a Research Assistantship from the NationalEngineering Research Council of Canada is gratefully acknowledged. Also acknowledgedare contributions by MacMillan Bloedel Research Limited Centre such as the use of theirlaboratories in the friction testing of Para1] am, and for providing other data on thiscomposite product.Chapter 1INTRODUCTION1.1 BackgroundTimber has been very important in North American bridge construction due to its abun-dance and the relatively unsophisticated requirements for equipment and skilled labourin comparison to steel and concrete bridges. More than a thousand timber bridges withshort and medium spans have been built across Canada. Two kinds of bridges, usingeither nailed-laminated or glue-laminated construction have been most commonly de-signed.In the first half of this century, with the high quality and high strength of steel andreinforced concrete materials being used more and more, timber as a material for bridgeconstruction was progressively ignored. Structural engineers were reluctant to use timberin bridges, since the frequent failure of such bridges in the past convinced engineers thatthese structures could not be expected to last more than 40 years [1], leading to expensivereplacements.In 1973, the Ontario Ministry of Transportation and Communications began a lami-nated timber bridge testing program to evaluate the load carrying capacity of this kindof bridges. The results pointed out that the load carrying capacity and service life were adirect function of the 'tightness' of the structural system. 'Tightness' was defined as theability of the structure to prevent relative movement between components and preventthe entry of harmful agents between adjoining interfaces [2]. The tightness also reflects1Chapter 1. INTRODUCTION^ 2the loading distribution ability of the structure. In turn, this ability is an indication ofthe bridge load carrying capacity.The bridge is a structure subjected to vehicle wheel loading, which is relatively con-centrated. Slips between components of the deck (cover) would occur and form gaps.Water and incompressible materials (such as stones) would enter these gaps. The incom-pressible materials would force the gaps to remain open and make the deterioration ofthe 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 capacityof the bridge and, consequently, increasing its loading carrying capacity. The transverselypost-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 loadtest results obtained before and after the transverse post-tensioning indicated the effec-tiveness of the procedure in improving the load distribution and decreasing the deflectionin the bridge [2]. Since then, the transverse post-tensioning method has been acceptedin new bridge design and construction.Timber bridges can be constructed in the four seasons of the year. This is especiallyimportant in Canada because of its winter weather. Their economic advantage has alsobeen proven by cost statistics for short span bridges [3].A method is needed to analyze the transversely post-tensioned timber bridge structureand obtain the relative interface movement under load. Such analysis could be usedfor the study of the reliability of the bridge for different performance criteria. Thedevelopment of such an analysis is the main purpose of this research.Chapter 1. INTRODUCTION^ 31.2 Literature ReviewPrevious research has developed different models for the two or three-dimensional analysisof stiffened plates.Cheung [4] developed the finite strip approach to analyze thin rectangular plateswith two opposite edges simply supported . A trigonometric series function was usedin the approximation of the deflection of the plate between the two simply supportededges, while a finite element representation was used in the cross-direction. For thiscase, the finite strip element method required less computer memory and increases thecomputational efficiency when compared to standard finite elements techniques.Foschi [5] used a combination of Fourier series and finite elements to analyze a type oftimber 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 ofthe T-beam flange is equal to the spacing between the beams. The flange (cover) isconsidered 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 thecover perpendicular to the span. However, the consideration of gaps caused couplingbetween the different Fourier terms in the global equations, requiring more computermemory 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 behaviourbetween cover and beams. The method considered the floor as a system of crossing T-and rectangular beams, including the effects of slip between layers owing to fastenerdeformation.Taylor [8] used the computer program ORTHOP to analyze the transversely post-tensioned Hebert Creek timber bridge, which assumed no transverse flexural stiffness inChapter 1. INTRODUCTION^ 4the 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 tak-ing into account the effect of transverse shear deformation. In their model, the transversesection was discretized by one dimensional finite elements and the longitudinal behaviourof the structure was defined by Fourier series. The element was called the simple twonoded strip element with one single integrating point.Harik and Salamoun [10] also adopted the strip element to analyze the stiffenedorthotropic rectangular plates. They idealized the stiffened plate as a system of platedstrips 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 betweencomponents must be considered in a more general analysis.1.3 Objectives of the ThesisIn order to obtain a general theoretical analysis, the non-linear properties for the trans-versely post-tensioned timber bridge must be considered. On the basis of the programFAP by Foschi [5], a non- linear, transversely post-tensioned timber bridge analysis com-puter program, PTB(Post-tensioned Timber Bridge), has been developed. The structureis 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. Thethesis objective is the development of a method for the calculation of the relative dis-placement between two adjacent T-beams, under the vehicle wheel loading. The methodand the PTB program can be used to obtain the relationship between maximum relativeChapter 1. INTRODUCTION^ 5Figure 1.1: Analytical Structural Modeldisplacement and post-tensioning force, permitting the calculation of the force requiredto achieve a target maximum relative displacement.The Newton-Raphson method [14] and the incremental loading method [15] are usedto solve the non-linear problem. Since the non-linear behaviour of the friction springsproduces coupling between the Fourier terms, the Jacobi iteration method [5] has beenused to solve the coupled global equations. As an example, this Thesis considers apost-tensioned bridge with T-beams made with Parallam, a composite wood productmanufactured by MacMillan Bloedel Ltd. of Vancouver, B.C..Chapter 2NON-LINEAR FINITE ELEMENT ANALYSIS2.1 Structural ModelThe bridge structural model consists of two types of structural elements, ie. the T-beamand the joint between beams shown in Figure 2.1. The joint element, as mentionedin Chapter 1, is assumed to represent the non-linear friction behaviour between theinterfaces of the adjacent beams.It is assumed that all the T-beam elements remain linearly elastic, with small de-formations. The joint elements have non-elastic material properties. For the T-beamelement, its flange can be assumed fastened to the web to produce an assembly capableof composite action, behaving as a stiffened plate under loading [5]. The flange may alsobe assumed to be rigidly connected to the web. Orthotropic thin plate theory [11] hasbeen used in the analysis, to consider the different elastic properties of the flange in thedirections perpendicular and parallel to the beams.A semianalytical procedure [12] has been used here, assuming in the x-direction aFourier series for the displacement function, and a finite element approximation alongthe y-direction. The T-beam and the joint elements are illustrated in Figure 2.2 . Thedisplacements 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-beamThese displacements are represented as follows:6S flange HSYT-beam Element1Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 7Joint Element X T-beam ElementwebFigure 2.1: The Structural ModelJoint ElementFigure 2.2: The T-beam Model and The Joint ModelChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS1) in the z-directionNW(X,y) =^Fin (y)sin(nrx )n=12) in the x-directionN nrxu(x,y) = E F2n (y )cos(  L  )n=13) in the y-directionN rv(x,y) = E F3n(y)sin( n Lx )n=1where:N = the number of terms used in the Fourier seriesL = the span of the bridgeThese displacements satisfy simply supported conditions at x = 0 and x = L. Thefunctions Fin (y), F2n(Y) and F3n (y) can be expressed in terms of polynomials in thenon-dimensional variable e = 2y/s, where: s = the spacing of the T-beams. A 5t h orderpolynomial is used for Fin(Y), and 4 th order polynomials for F2n (y) and F3n (y), withdegrees of freedom associated with nodes 1, 2 and 3 in Figure 2.2. Thus, if {8 n } isthe vector of nodal displacements, Fin (e), F2,() and F3n (e) and their correspondingderivatives can be written as:F1n (S) = {M0(e)}T {sF2n(0 = { 114-3()}T { bn}F3n(e) = { ill5()}T { Sn}8(2.1)(2.2)(2.3)with corresponding derivatives:N8(x) = E e4nsin( nrx-)n=1(2.14)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 9dFin(e )1 = {MAW {an}ckd2 F1n(e) = {M2(e)}T { 8n}ck 2dF2n(e) = { 1114(e)}T { 8n}ckdF3n() = { 1116(e)}T {6n}c14.The vectors {M1 (01 to {M6 ()} are given in Appendix A.(2.10)2.1.2 Displacement functions for the centroid of the T-beam web(2.7)(2.8)(2.9)These are expressed as follows:1) in the z-direction2) in the x-direction3) in the y-direction4) rotation about x-axisN) :"--- E W4nsin ( n.trx )n=1(2.11)NU(x) = E u4ncos(n:x )^ (2.12)n=1NV(x) = E V4nsin(7)^ (2.13)n=1"W(The vector {8,j, associated with the nth term of the Fourier series, is shown inEquation 2.15. Each T-beam element has 19 degrees of freedom.Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 10U inn ln • SyinVln • SU2nV2nW2n • SW4nU4n"Kin94n • 3(2.15)IQ = <W3nW3n • 8U3n1U3n • Sv3nvan • sThe prime (') denotes the derivative with respect to y. For example, wi n • s is thederivative of w(x,y)with respect to y at node 1 associated with the nth item of Fourierseries. The rotation degree of freedom is multiplied by the beam spacing s to make allcomponents in {6n} dimensionally consistent.Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 112.1.3 Degrees of freedom for the joint elementSimilar to the T-beam element, the degrees of freedom for the joint element are associatedwith the each term of the Fourier series. =winW ln • 5UlnV inW2nw2n •£(2.16)2nV2nThe joint element consists of discrete three-dimensional springs at both the top andthe bottom of the flanges along the bridge span shown in Figure 2.3. The dots showthe positions of the three-dimensional springs which connect the flanges of two adjacentT-beam elements. The three-dimensional spring can act in the x-, the y- and the z-direction.The constitutive relationships for the x-, y- and z-direction springs are illustrated inFigure 2.4, where:fix = the friction coefficient of the T-beam flange in the x-direction;Az =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;Three-dimensional Springo^Ayy-directionFyE;Fxii,a,e t.paye tx-directionChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 12Figure 2.3: The Distribution of the Discrete SpringsFigure 2.4: The Constitutive Relationship of the SpringsChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 13Az = the elongation of the spring in the z-directionThe elongation of the spring is represented by the relative displacement between thetwo 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 springloses its stiffness in the corresponding direction, but maintains a constant force equal tothe friction force.The spring in the y-direction is used as a contact spring, to permit flange separationbut not overlapping. When the relative displacement is negative or zero, ie. the twoflanges tend to overlap each other, the stiffness of the y-direction spring is set veryhigh, to enforce no overlapping between them. Thus, when the relative displacement ispositive, which means that there is a separation between two adjacent T-beam flangesat the spring position, the spring is inactive or has no stiffness and no force in it. Whenthe spring in y-direction loses its stiffness, this spring will not make contribution to thestiffness 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-directionare not active.2.2 Strain Energy of the Structural System and Virtual WorkThe energy variational approach and the principle of virtual work are to establish theglobal equations for the structure. In this section, the strain energy of each compo-nent and the virtual work done by the conservative and the non-conservative forces arepresented in terms of nodal displacement degree of freedom.Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 142.2.1 Strain energy in the flange of the T-beam elementThe strain energy per unit area of flange can be obtained by applying small deflectionorthotropic plate theory.Ufu^Kx 8 2w 2 Ky 82w 2^a2w 82w2 ( axe + 2 ( 0y 2^Kv ( axe )( ay e02w Dx au 2 Dy av 2+2K00x0y ) 2 + 2(0x ) + 2(Oy)( 0u0v 1 DG au 07) 12ax ay ) + 2 `. ay OxFor the plate with the thickness d, the stiffnesses Kx ,Ky ,K„,KG, D.,Dy ,Dvare given as follows:1) the flexural stiffness in the x-directionKx =12(1 —2) the flexural stiffness in the y-directionEnd3 Ky =12(1 — vxyvyz )andKv = vzyKx3) the torsional stiffnessKG -Gd34) the axial stiffness in the x-directionDx =(1 — vxy vyx )Exd312Ex d(2.17)and DGChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^155) the axial stiffness in the y-directionEydDy =^(1 — vxy vyx )andD, = vxy Dz6) the shear stiffness in planeDG G•din 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.Uf^2 0L UfudxdyF (2.18)Substituting Equation 2.17 and derivatives of the w(x, y),u(x, y) and v(x,y) (shownin Appendix A) into Equation 2.18, the total strain energy in the flange of one T-beamelement is obtained in terms of the vectors {5,J. The integration makes use of theorthogonality of the trigonometric functions in the interval (0,L):r sin( nr: )sin( miarx^ )L m = n20 n nChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 16joL cos(nrL x )cos(m Lrx)m = n= { 20 Tit ny= —s, de = -2 dy, dy = -s2-deThus,8Uf = E Ufi^(2.19)i=1where the Ufi components are given as follows:Uf 3A^L= -Uf i^12. {1 Ufu ldX}ClYr 0N icn4 7.4 .9 1. 118L 3^-E^ {sn}T{m0(e)}{m0(0}T{sn}den=iLUf2 = 1 2. {1 Ufu2dX}dy2:FC ) 3 /1 {Sn} T {M2 (0}{ 1112 W}T {En}4n=i 2 S 2 -1I• L2.{fu3dX}dy= E --Kv( 77-1 ) 2 ( .9 ) 1- -121^{5.}T{M0(e)}{M2()} T{Sn}deL n=1(2.20)(2.21)(2.22)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 17{^U fu4dX}clyUf4^0Za J2KG (LN^nr s L t1)2(-2 )—2 --1{8n}T{M1()1{M1(e)}T{an}ck(2.23)Uf5^Ufu5dX}elyE n s L fiT ( '7 -Lr-)2() 2 J- 1 {871}T{M3V)}{M3(°}T{Sn}4(2.24)Uf6^1 L1{1 Ufu 6dX}cly0/1 fri}T{M6(0}{M6(0}T { 8.}4n=i 2 s 2 -1Uf7^a{ 0 Ufu7dx}dy2Dv(17r-)—L{8n}T{M3()}{M6(E)} T{8.}4n=1^L 2LUf 8^ 1= { i Ufusdx }dy1 02DG SL .1 1 {5n }T [XX]{8„}ck2 4(2.25)(2.26)(2.27)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 18Where:[xx] = ( .-9-2 ) 2 {m4 ()1{ 1114()}T + ( 7-7- ) 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 webThis is given in terms of the displacements U, V and W and the rotation O. Thus,Ely^d2W 2,^EIz^, d2 V 2^EA fL,dU 2^GIt 1,d8 , 2Uj =   dx+    dx+ dx+^dx (^)2.292 Jo dx 2)^2 Jo dx 2)^2 o dx )^2 o dxWhere:^Iy = (BH3 /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 alarge number, eg. E/G = 17.0.Substituting the derivatives of the W(x), U(x) V(x) and 8(x) (see Appendix A ) intoEquation 2.29 the strain energy of the web can be written as:Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 194uj^i=1^ (2.30)in whichEjlUj2U j3=EIy r( d2W ) 2 dx2 Jo dx 2n=1 m=1EIy v (winlor vi( Ln2 =1) L ) ( 2EIz I,( eV )2dx2 o dx 2EIy .‘1\r\ ■‘1‘r I L vnv, 4m ( nr )2 ( mr \ 2sin( nrx )sin( mrx2 y ^m=1Jo^L^L I^L^LEIz N LE(V4n ) 2 ( 17  )4 (2 )2 n=1 ^2EA iL ( dU )2dx2 Jo dxEA \1v, 2.",rIL u4nu4m(nr y mr \sin( nrx , s .n , mrx2 L-'d ‘---d o^L ( L )^k L ) I ( L^ )dxn=1 m=1EA N2A ^)n=1^EIy^fL^m7r^inrx ,s.n,mrx )dx^2 y^" L 2( L 12sin ( L ) ( L )a)dx(2.31)(2.32)(2.33)(2.34)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 202.2.3 Strain energy in the web-flange connectorThree kinds of deformations have been considered in the connectors between the flangeand the beam1) The deformation of the connector parallel to the beam:emu= ru2_ d ( aw2)1^H dW1[^2 Ox^[^2 dx2) The deformation of the connector perpendicular to the beam:Ay^[V2 d ( aw 2^ )]^[V + I 19]2 ay 23) The rotation deformation of the connector about the x-axisa2lv= ( ay ) 9Where:d = the thickness of the flangeH = the height of the webu2 = the displacement in the flange at point 2 in the x-directionv 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(2.35)(2.36)(2.37)The flange is assumed connected by uniformly spaced discrete connectors along thelongitudinal centerline of the beam. Thus, the total strain energy in the connectors ineach T-beam is:UN = Au) ?. + 14(64 +^ (2.38)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 21Where:NA = the number of the total connectors in one T-beam;(Au) i = the deformation in the x-direction for the ith connectors;(Av) i = the deformation in the y-direction for the i th connectors;(00) i = the rotation deformation for the ith connector about the x-axis;kz = the stiffness of a single connector corresponding to the x-direction;icy = the stiffness of a single connector corresponding to the y-direction;ke = the rotation stiffness of a single connector.Alternatively, an equivalent continuous connector model can be used to calculate thestrain energy of the connectors per T-beam.UN = IL 12ekAu\2 ky v ,2e^2e(A0)2]dxo where e = the connector spacing along the beam.(2.39)Substituting Equation 2.1 - Equation 2.3 , Equation 2.11 - Equation 2.14 and thecorresponding derivatives into Equation 2.35 - Equation 2.37,NAu^E [F2n (y = 0)cos( nrx ) d Fin (y = 0) cos(^ )rtrxL^2^L^Ln=1U4n cos ( nrx )— Tir nir cos( nrx )1-2- vv 4n L^L(2.40)UN =^ 2 2eE^[u2n —^— W4n(H d)]n=1 N L 1 2Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 22N d dFin (y = 0)^nrxAv = E [F3n(y = 0)   I74n^e4n]sin(^ )n=1 2^dy 2^LN^d ,^ nrxE [v2. — 2 w2n—^ — —2 94n]sin(  L  )n=1 (2.41)N nirxAO =^[w2n 04n]sin(—Tn=1Finally, introducing these three Equations into Equation 2.39,(2.42)d ,^H n—+ky [V2n —211)2n V4nn+k9 [w'2n — 04,] 2 } (2.43)In order to express the strain energy of the connector in terms of the nodal degreesof 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 }Tinto the Equation 2.43, the strain energy of the connector can be written as follows:Where:N L 1UN = E 2^ke2e{k{X} + ky {Y} + —{0}}n=1(2.44){X} = {8n } T [{e 7} — {e n } — {e io }(H + d)] •[{e 7} — {en } — {e io}(H + d)] T {8n }^(2.45)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 23{Y} = {an}T [{e8} fegl§ }{ e12  — 13 _ 2sdT[{e8} {69}-279 —1' 1— {613}- 51 {an} (2.46){O} = { 5n}T Re9} — {e13}1.Re91 — {en}i{sn}^(2.47)2.2.4 Virtual work done by non-conservative forces in the joint elementu(x i , zi )^u(x i ) — ziThe joint element consists of 2N8 'three dimensional' springs. The deformations of the ith`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-coordinateis 0.5d when it is at the bottom of the flange.Along the x-direction,with the spring elongation then being,5wC aw  .(2.48)Au(x i ,zi )^u2(xi,zi) —u 2 (x i) zi ( 5w1 )Oxxizi (2.49)Similarly, along the y-direction,v(x i , z i )^v(x i ) — zi OwCa (2.50)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 24with elongationAv(xi,zi)^v2(xi,zi)— vi(xi,zi)For the z-direction,( awl^(aw2—Vi (Xi ) + V2(Xi) +^z,ay) ay) xi(2.51)w (x i , z i )^w(x,)^ (2.52)and the spring elongation is:Aw(x i ,^=-- w 2 (x i ,zi ) —-=^w2(xi)^(2.53)Substituting the displacement functions of the T-beam flange and the correspondingderivatives in Au(x i ,zi ), Av(x i ,zi ) and Aw(x i ,zi ) above,N^ nir^fir^ nrxi)^Au(x i , zi) E[—u1n+u2r, +  w1n zi^w in ]cos ( (2.54)n=1nirxi)Av(x i ,^= E[—vin V2n ZiWln ZiW2n]Sin(^ (2.55)n=1NAw(x i , zi) = E[—w1n + w2n]sin(nrxi )n=1which can be written in terms of the nodal degree of freedom vector.(2.56)(2.57)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 25N^N= EfAnl = EAdfs'l^(2.58)n=1^n=1Where:{{Aun }{An} = {Avn }{Awn }1 (2. 59)zi (-yr )cos(Ti )^0^—sin(7-71)0^7-9isin(n)^0--cos ( mr: i )^0^00^—sin( mr i )^0- Zi (n7 )C08( 7"---ri)^0^sin(7-i-)0^—isin(Ti)^0cos ( 1"---1--7i )^0^00^sin(12T)^0The forces in one 'three-dimensional' spring will be:N{F} = [D]{A} = [D]^[Bn ]{8,J,}n=1Where:EI 0 0[D] = 0 E; 00 0 E:= the stiffness modulus of the spring in the x-directionEy = the stiffness modulus of the spring in the y-direction[Br,]T = (2.60)(2.61)(2.62 )Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 26E: = the stiffness modulus of the spring in the z-directionAssuming that the structural system is subjected to a virtual displacement, corre-sponding 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=1The internal virtual work in one 'three-dimensional' spring is, therefore,T47 = [{A}1 T{F}N= E({87.T)T[Bn]T{F}n=1N N= E E (feLY)T [Bn]T [D][B.]{8;17,}n=1772=1The internal virtual work in one joint element is then the sum(2.64)2N.= i=1^ (2.65)where 2N8 = the number of the total three dimensional spring in one joint element.2.3 Global EquationThe global system of equations for the structural system is[G]{U} {P}^ (2.66)in which:[G] = the global stiffness matrixChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 27{U} = the global nodal displacement vector{P} = the global consistent load vectorThe global stiffness matrix, the global nodal displacement vector and the global consistentload vector are assembled by the total contributions from each individual element usingthe conventional finite element method.2.3.1 Stiffness contribution from one T-beam elementThe stiffness contribution of one T-beam element can be derived by minimizing the totalstrain energy in the T-beam element. The total strain energy of the T-beam is thesummation of the strain energy of the flange, the strain energy of the web and the strainenergy of the connectors.U = Uf +Uj + UN8^4^3= E Ufi E uji + E uNi^ (2.67)i=1^i=1^i=1Take the first variation of the strain energy of the flange with respect to the {Sn }, weget the stiffness contribution from the T-beam.8(Un) = a(Uf )71 + 8 (Uj )7/ + 8 (UN)n8^4^3= E(sumn + E(suji )„ + E(SyNoni=i(2.68). n4r4s 1a(ifi)n^K 4L3 [/ if-M-0(0}{M0(0}Tda5.1(2.69)l8(Uf2). = - . ) 3L[fP^i{M2(e)}{M2(01T41{6.}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}a (Uf 5 )n = 221- (7) 2;L[fli {M3()}{M3 (e)}741{8n}8(Uf6)n = 2!-L{111{114-6()}{M6()}Tde]{an}8(Uf7) 71 = D, (74 [111{ m3 (4)}{ m6 (4)}Td4ll sri l(2.72)(2.73)(2.74)(2.75)8(Ufs). =DG( 2 )L[ li ( s2 )2{ m-4 ( ) 1{ 14 (0}T ig + 1 1 ( 7 )2{m5 ()}{1115 ( )}Tdenr 2^ ri 7r 2+ .1_ 1 (i);{M5(0}{M4(e)/T4 L i ( 1217);{M4(0}{M5(0}Tckif8n}(2.76)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 294E^nr ^I^Trir^W 4n( L^) 4 L + E 2^z V4„(—L )4L+EA u4„( nr ) 2L GIto4n(—nr ) 2L2^L^2^L (2.77)3E^2 —L flz[X][.X]T ky[Y][Y]T —ke [0][0]T}{8n}e1=1 (2.78)[X] = [fe 71— {en } — {e io}(H d)]^ (2.79)r^r^i H i[Y] [{e8} len/ le131 .0 (2.80)[0] = [{ e 9 } — { e 1 3 }]^ (2.81)2.3.2 Stiffness contribution from one joint elementBecause there are non-conservative forces in the joint element, the corresponding stiffnesscontribution is obtained by using the principle of virtual work.As shown in the previous section, the internal virtual work done by the spring forcesin the joint element is:2N. N N(Wirt 8 = E E E ({8T)T[Bn]T[D][Bm] le„}i=1 n=1 m=1Since {8;1 }* is an arbitrary nodal displacement the stiffness contribution from thejoint element corresponding to the n th and inth terms in Fourier series is:(2.82)P(x,Y)PTy2 y1Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 30PTXFigure 2.5: The Vechicle and Post-tensioning LoadY2N,BK[n,m]^[B„]T [D][Bm ]^ (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 andall springs are elastic, since such limit approaches the orthogonality situation among thetrigonometric shape functions.2.3.3 Consistent load vectorsTwo kinds of loading are shown in Figure2.5.Vehicle wheel loadUnder a virtual displacement w - (x,y), the potential of the vehicle wheel load p(x,y) isUi = 1 °2 I Y2 fulx,y)}T {p(x,y)}dxcly01^yiChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 31j(x2r2^ r=^{ {Mo()} 7' {S„}*}T {p(x,y)}sin(nLx )dxdyxi j Y1 n=1From the definition of the consistent load we have>{{Sn}*}T{Rn} =x2 Y2 {iv*(x,Y)}T{P(x,Y)}dxdYn=1^ Yi1V(2.84)- 1M 2 I Y2N°1IE{{84*}T{Mo(E)}{P(x,Y)}sin(nr X )dxdy(2.85)Then:{R.}= f 1 r{m0(0}{pcx,yilsin(n^ )dxdyx ^yi(2.86)where: {Rn } = consistent load corresponding to the n th term in the Fourier series.If p(x, y) = P, a constant uniformly distributed load over the patch, {R.} = P^ iziz2 1: 2 {Mo(e)}sin( n:x )dxd2Where:P s L^nirx2 nr {cos( Li) cos( n7 2 )} 166 {Mo()lck (2.87)^2y i^2y2el —;^=s sPost-tensioning loadThe individual concentrated post-tensioning load is assumed to be applied, respec-tively, at node 1 and node 3 of the first and last T-beams. Under an arbitrary displace-ment e(x,y), the potential of the post-tensioning load PT(x i ) isNTUPT = Elyi=1( ,O}T {PT(xi)} (2.88)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 32From the definition of the consistent load we have :NTE{{87,}1T {R,} = Et{8n}lT {M5()1 E PT (x i )sin( nrx^)^(2.89)n=1^n=1^i=1NT^nrx .{R„,} 5 = E sin(^')PT(x i )^ (2.90)Li=1NT^nrxi{R-.}18 — E sin(  L  )PT(x i )i=1where:PT(x i ) = the value of the ith 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.91)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:[AKe]{51 = {Re} (2.92)where:[AK (1,1)]^[0]^. [0][AK'] =[0]^[AK(2, 2)]^... [0](2.93)[0]^[0] [AK(N, N )](2.94){ 6-eN}Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 33{R7}={Re } (2.95){RN}Each submatrix AK [n, n] is a 19 x 19 symmetric square matrix. N is the number ofitems in the Fourier series. {Sn} and {R„e } ( n = 1, 2 ... N ) are the 19 x 1 vectors.2.4 Global Equations and Solution MethodThe characteristics of the global equations and the solution strategy will be discussedwith an example, using only two terms in the Fourier series and three T-beams ( NJT=3 ).2.4.1^Characteristics of the global equationThe form of the global equations is shown in Figure 2.6 and Equation 2.96:I[GK(1,1)] PK(2,1)][GK(2, 1)] PK(2,2)]_{U1} = {P1 }{U2} j^{P2}(2.96) AK 1 (1, 1) is the stiffness contribution from the T-beam element 1 associated withthe first term of the Fourier series. As mentioned in the previous section, it is a 19 x 19symmetric 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-thterms in Fourier series. It is a 8 x 8 symmetric matrix.Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 34n = 1^n = 2 n = 1n = 2P2Figure 2.6: The Global Equation SystemGK(n, n)(n = 1, 2) is the (NJT * 19) 2 symmetric banded matrix with a band widthof 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 nth term of theFourier series2.4.2 The incremental procedureThe incremental method [15) is used to solve the non-linear global equations. Theinitial values play a very important role in the non-linear analysis iteration method. TheChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 35incremental 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} = E3i=1 {AU2 }Where:{AUi} = the incremental displacement of each step.{APi } = the incremental loadThe initial value used in order to get {Ui} is the solution { Uj_ 1 } in the previous step.With in each step, the global stiffness matrix must be updated according to the currentdisplacement. The final solution is the summation of all displacement increments. Theprocess is illustrated in Figure 2.7.2.4.3 Newton-Raphson methodFrom the Figure 2.7 we could note that the more the number of steps, the more accuratethe solution is.In order to decrease the error within each incremental load step, the Newton-Raphsonmethod 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-Raphsonmethod converges. Thus,{AUi} = E{Auk }k=1Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 36A P 2APIU 1 A U2Figure 2.7: The Incremental MethodAAP14, up 2A U 1^AFigure 2.8: Newton-Raphson MethodChapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 37[G(U)]t{Auk} = {AQh}^ (2.97)I I { Auk} II 2 < cid 1 1{uk _1 } 112where:[G(U)] t = the global tangent stiffness matrix which must be updated to the currentdisplacements;{AC} = imbalance load vector;Eid = a tolerance;{uk - 1 } = Etil{Dui};Il{uk}112 = {Ei= 1 uZi } 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{R1-1 } = {p} — [K(ui- 1 )]{ui- 1 } 4 0^(2.100)We look for {uk } and make {Rk } equal zero{uk } = {uk_i} + {Auk}^ (2.101)Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 38{Rk} = {R(uk _ i + Auk)}= {R(uk )}^0 (2.102)Expand {Rk } in Taylor series in the neighbourhood of {uk } and only take the firsttwo terms.„^19{R} ,{R(uk _ i + Auk)} = {R(uk-1)1 (9{u}{uk-i} • {Auk} = {0}(2. 1 03)(NM-1ft,^{Auk}(9{ 14 k-1 {R(uk-i)}{p} — [K(uk_i)]{uk-il (2.104)where:(9{R} Ifuk^= ig(U)lta{u}[g(u)] t is the element tangent stiffness matrix corresponding t o the {uk _ 1 } and {uk _ 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 elementlevel.At the global level Assemble all element tangent stiffness matrices [g(u)] t and imbalance load vectors {q(u)}byconventional finite element method, we get the global tangent stiffness matrix [G(U)]tand the global imbalance load vector {.64k}.Chapter 2. NON-LINEAR FINITE ELEMENT ANALYSIS^ 392.4.4 Jacobi iteration methodWithin each iteration of the Newton Raphson method, the Jacobi iteration method isused to solve the global system of Equation 2.97. This avoids storing the entire globalmatrix, requiring less computer memory . If the diagonal matrices are dominant, therate of the convergence of this method is very high. The procedure can be written asfollows:N{ Al k};n = [CK(m,m)] t {{AQic}ni — E[GK(m,n)t{Auk}in 1}n.1with starting vector{Auk };„ = [G K (m, m)]t {AQ k}n(n m) (2.105)(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,...,Nwhere: E is the error tolerance (which could be .001 for example ).Theory presented in this Chapter has been incorporated into computer program PTBgiven in Appendix B.Chapter 3FRICTION PARAMETERS FOR THE SPRING MODEL3.1 The Friction TestIn order to determine the friction characteristics of Parallam, an experiment was con-ducted at MacMillan Bloedel Limited Research Centre. Four kinds of surface textureswere 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 sizeof the bottom specimen was 88.9mm x 88.9 x 38mm Both were cut from rough-sawnplanks. It should be pointed out that the surface of the preservative treated specimenswas 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 measuredbefore being soaked in water and weighted just after being taken out of water, and wipedclean of surface water. For each kind of surface texture, both the friction coefficient forthe sliding direction perpendicular and parallel to Parallam fibres were determined, asshown in Figure 3.2.40Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^41slidingbottom specimenweightclamp^ "`- top specimenFigure 3.1: The Friction Test Setupdirection slidingParalleldirection of slidingPerpendicularFigure 3.2: Two Kinds of Friction TestChapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^42Table 3.1: Data of Friction Test SamplesSpecimen Type Average Weight (gr) Number of SpecimensDry, Treated, Parallel 110.99 20Dry, Treated, Perpendicular 111.16 19Wet, Treated, Parallel 141.15 17Wet, Treated, Perpendicular 140.08 14Dry, Untreated, Parallel 92.36 15Dry, Untreated, Perpendicular 92.75 12Wet, Untreated, Parallel 135.18 19Wet, Untreated, Perpendicular 137.47 12Weight of the Clamp 109.91Applied Weight 5000.0Table 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 theParallam fibres. In the parallel direction test, the sliding direction was parallel to thefibres.The force FH required to produce sliding between the surfaces of the specimens wasproportional to the force applied normal to the plane of motion. The ratio is definedas the coefficient of friction.FH=^^ (3.1)FvIn all tests, an Instron 4210 testing machine was used to provide a 0.4in/min. constantsliding velocity. The horizontal force FH vs. sliding movement relationship through thetest was monitored and recorded, and typical curves are shown in Figure 3.3 to 3.5. Fromthese 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. Theresults and statistical data are presented in Table 3.2.Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^43Table 3.2: Statistical Results of the Friction TestClassified Specimens tt a of it Ct, of it E ( N/mm ) a- of E CC, of EDry, Treated, Parallel .655 .076 .115 16.902 2.895 0.171Dry, Treated, Perpendicular .831 .080 .096 17.378 1.922 .111Wet, Treated, Parallel .833 .038 .046 12.618 3.197 .254Wet, Treated, Perpendicular .878 .059 .068 12.353 3.863 .313Dry, Untreated, Parallel .373 .031 .083 11.183 1.751 .157Dry, Untreated, Perpendicular .419 .041 .098 10.891 1.992 .183Wet, Untreated, Parallel .804 .052 .065 14.900 2.523 .170Wet, Untreated, Perpendicular .811 .048 .059 12.339 1.943 .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 ModelAs previously mentioned, a 'three-dimensional spring' model is assumed in the non-linearanalysis and its constitutive relations are shown in Figure 2.4 of Chapter 2. The frictioncoefficients Ax 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 'sameslip' 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,fixcryetAlien =^E;(3.2)Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^44Figure 3.3: Friction Test Recording Curves(1)F(lb)Dry'Untreate::Parallel40.5Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^45Dry TreatedPenendicularF(1, )2) DEtn).0 5 1.0Figure 3.4: Friction Test Recording Curves(2)Wet TreatedParallel,D(inChapter 3. FRICTION PARA1VIETERS FOR THE SPRING MODEL^46^Ff ItsFmet Treated_Perpendicular'0.5^1.0^ 0.5 1.010 ^fi1-Figure 3.5: Friction Test Recording Curves(3)Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^47Figure 3.6: Idealized Friction Test Curves= AzEo- su etIn the parallel direction of the friction test:FHAix = ElIn the perpendicular direction of the friction test:FH=z zThe 'same slip' assumption implies that= Al. (3.6)Atim = Aix^ (3.7)(3.3)(3.4)(3.5)from which,Chapter 3. FRICTION PARAMETERS FOR THE SPRING MODEL^48 0=Figure 3.7: Same Slip AssumptionE t cr etx Y El =FT'(3.8)^-Et cr etz Y E: =Fir (3.9)E; and E.; are the stiffnesses used in the 'three dimensional' spring; Ext and E are thestiffnesses obtained from the friction test. The relations between the spring stiffness andthe stress ay is clearly shown in Figure 3.7. We can note that E; and E; are proportionalto the normal stress cry , however the maximum friction force in the spring cannot exceedthe shear strength of the material rs et, where Ta is the shear strength of the flange andet is effective area of the spring.Chapter 4ANALYTICAL RESULTS4.1 Beam TestIn order to test the spring model described in Chapter 2, a simple beam structuralproblem, shown in Figure 4.1, is now considered. The structure consists of two built-in beams connected by four springs. At the top, there are one y-direction and onez-direction spring. The other two are assumed in the y and z-directions at the bottom ofthe beams. The left beam is subjected to one vertical concentrated force at midspan. Theproperties and the dimensions of this structure are shown in the Table 4.1 The structurewas analyzed by a specific program using the same spring model theory, which is used inthe PTB program.When P > 0.0, this axial load produces tension in both beams. The results givenin 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 theTable 4.2In Table 4.2, F means the force in the spring is beyond the friction limit; E meansthat the deformation of the spring is within the elastic range, less than the friction limitAtim•From Figure 4.2 and Table4.2, we can note that increasing the post-tensioning force Pinfluenced the status of the springs, with an associated change in load distribution withinthe structure. Finally, when the P reached 100, 000(N), the two y-direction springs were49Chapter 4. ANALYTICAL RESULTS^ 504 L BLFigure 4.1: Two Cantilever Beam StructureTable 4.1: Properties of the Cantilever BeamBEAMwidth of the beam ( B ) 200 (mm)depth of the beam ( H )^. 400 (mm)span of the beam ( L ) 4000 (mm)elastic modulus of the beam (Eb) 14000 (MPa)elastic modulus of the y-dir. spring (E y ) 1.0E10 (MPa)elastic modulus of the z-dir. spring (Er ) 17.4 (MPa)friction coefficient (/./ z ) .65LOADINGvertical concentrated loading ( V ) -20,000 (N)Table 4.2: Status of the Springs in the Cantilever BeamsSTATUS OF THE SPRINGSPost Force (N) y-dir. at top z-dir. at top y-dir. at bottom z-dir. at bottom0.0 contact F open not active10 contact F open not active100 contact F open not active1,000 contact F open not active10,000 contact E open not active100,000 contact E contact E1000,000 contact E contact EChapter 4. ANALYTICAL RESULTS^ 51VERTICAL DISPL. vs. POST TENSIONING FORCE ." 4.002.00CD 0.00a)8-2.00 -co-4.00 -aCD -6.00 -'C^o^CD-810.8E+00NODE 2--B--NODE 31.0E+02^1.0E+04^1.0E+06^1.0E+08^1.0E+10Post-Tensioning Force ( N )Figure 4.2: Cantilever Beam Structure DisplacementsL/2L3Figure 4.3: Built-in Solid Beams1 Z y 2 3Joint 1^Joint 2Chapter 4. ANALYTICAL RESULTS^ 52Figure 4.4: Three T-beam Structureactive ( full contact between the beams) and the deformations of the two z-directionsprings ( in terms of relative displacement at the interface of the beams ) were all withinthe elastic deformation limit. When the P reached 100, 000(N) the displacement ofnode 2 was decreased by several times in comparison with that for P = 0.0(N). Thedisplacement of node 2 must then be nearly equal to the displacement of node 2 inFigure 4.3. The structure in Figure 4.3 is a solid beam and all its properties are the sameas those for the case of Figure 4.2. This was independently verified with a structuralanalysis of the beam in Figure 4.3.4.2 T-beam AnalysisThe configuration and the other properties of the structure are presented in Figure 4.4 andTable 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^ 53Table 4.3: Properties of T-beam StructureParameter Notation ValueT-beam span L 12000 ( mm )T-beam spacing s 1000 ( mm )web height II 500 ( mm )web width B 100 ( mm )the ratio of elastic modulusto shear modulus of the webREG 17.0flange height d 100( mm )friction coefficientin x-direction F. 0.50friction coefficientin z-direction fiz 0.65three dim. spring numberin one joint element 2N, 42spring elastic modulusin the x-direction Ex 16.9 (N/mm)spring elastic modulusin the y-direction EY 1.0E10 (N/mm)spring elastic modulusin the z-direction Etz 17.4 (N/mm)web elastic modulusin the x-direction EL eb 14000.0 (MPa)web elastic modulusin the y-direction EYweb 1400.0 (MPa)flange elastic modulus E flange 14000.0 (MPa)post-tensioning force spacing 1200 (mm)PTLPTChapter 4. ANALYTICAL RESULTS^ 54Figure 4.5: Three T-beam Structure with Post-Tensioning Force Only(deflection) of the structure, simply supported boundary conditions are satisfied at x = 0and x = L. Other edges are free. The analysis considered five loading cases and was doneby 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 wasto study the effect of the post-tensioning force in the loading distribution within thestructure. The program PTB can give either nonlinear or linear elastic results, the latterimplying that all springs keep their initial stiffness throughout the loading. Studyingthe relative displacement in z-direction between the flanges of adjacent the T-beamsconstitutes the main objective of this analysis.4.2.1^Analysis for post-tensioning force onlyThe structure and loading condition are shown Figure 4.5Three T-beams were subjected to ten pairs of post-tensioning forces along the span.Chapter 4. ANALYTICAL RESULTS^55Table 4.4: Analytical Results: Post-tensioning Force OnlyPTB PROGRAM SAME SLIP ( LINEAR)POST FORCE( N )MaximumDisplacement( T-beam 1)MaximumBEND. STRESS( T-beam 1)MaximumDisplacement(T-beam 2 )MaximumBEND. STRESS( T-beam 2 )( mm ) ( MPa ) ( mm ) ( MPa )1.00 -0.57952E-05 -0.14474E-05 -0.58064E-05 -0.11604E-0510,000 -0.57726E-01 -0.14681E-01 -0.58461E-01 -0.11126E-01PTB PROGRAM SAME SLIP (NON-LINEAR )POST FORCE( N)MaximumDisplacement( T-beam 1)MaximumBEND. STRESS( T-beam 1 )MaximumDisplacement(T-beam 2 )MaximumBEND. STRESS( T-beam 2 )( mm ) ( MPa ) ( mm ) ( MPa )1.00 -0.57952E-05 -0.14474E-05 -0.58064E-05 -0.11604E-0510,000 -0.57726E-01 -0.14681E-01 -0.58461E-01 -0.11126E-01From Table 4.4, we can note that the results from the non-linear analysis and those fromthe linear analysis are identical. This means all springs were in the elastic range at thedifferent levels of the post-tensioning forces.Figure 4.6 gives the relative deformations of the joint element 1 in z-direction alongthe span of the T-beam. The negative deformation of the joint element indicates thatthe 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 loadingThe load case is shown in Figure 4.7. This is a symmetric problem. Along the 12m longspan 10 pairs of post-tensioning force were acting. Each pair of the post-tensioning forcewas equal to PT. The structure carried one patch of vertical loading equal to 500 mmx 300 mm x 0.2 (N/mm 2 ) = 30,000 (N) at midspan of T-beam 2. For the symmetricChapter 4. ANALYTICAL RESULTS^ 565E-080-5E-08-1E-07-1.5E-07E• -2E-070-2.5E-07as-3E-07.S^-3.5E-07 0EUco0.cn0. 77^-0.0002Ti• 4.00042,000^4,000^6,000^8,000^10,000^12,000The x-coordinate along the span (mm)PT = 1.0(N)-0.0006-0.0008-0.001-0.0012 0^2,000^4,000^8,000^8,000The x-coordinate along the span (mm)PT = 10,000(N)Figure 4.6: Displacements of T-beam Structure under Post-Tensioning Force Only10,000 12,000Chapter 4. ANALYTICAL RESULTS^ 57Figure 4.7: T-beam Structure One Patch Loading Caseproblem , we only plot the stresses and the displacements of the T-beam 1 and T-beam2. 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 inbeams.When PT = 0.0, the total vertical loading is carried by beam 2 only. The maximumrelative displacement between beam 1 and beam 2 ( ie. the maximum deformation ofthe Joint element 1 in the z-direction ) was more than 20mm. With the increase in thepost-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 tothose from the linear elastic analysis. When PT reached 100,000 (N), the maximumrelative displacement, in z-direction, between the flanges of beam 1 and beam 2, wasnegligible when compared to the 20mm for PT = 0.0.Chapter 4. ANALYTICAL RESULTS^ 58coa.110ET-beam 1T-beam 2(D°3co6 -.5ou)4-ED. 2 ...O. •(.7)0 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07The Post-tensioning Force ( N )Figure 4.8: Bending Stress of T-beam Structure One Patch Loading Case• 2520EX15coco 105 5A 0lcp76'9-• iE+00LJT-beam 1T-beam 21E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07The Post-tensioning Force ( N )Figure 4.9: Displacements of T-beam Structure One Patch Loading CaseChapter 4. ANALYTICAL RESULTS^ 5920a)USI 150.2 .6 10>*fa'cl" a+oocc5Joint Element 11 E+06 1E+071E+01^1E+02^1E+03^1E+04^1E+05The Post-tensioning Force ( N )Figure 4.10: Relative Displ. of T -beam Structure One Patch Loading CaseIn the example shown in Figure 4.9, upward (negative) deflections of the T-beamscan be seen to be developing as the post-tensioning force becomes large. The reasonfor this behaviour is explained as follows: when the flange of the T-beam subjected topost-tensioning forces in the y-direction, the flange will expand in the x-direction dueto Poisson's effect. Since the flange is connected to the T-beam web by the connectorswith non-zero stiffnesses, shear forces develop to maintain compatibility. The directionof these forces as shown in Figure 4.11, is such that they cause an upward deflectionwhich, shown in Figure 4.9, becomes more important for post-tensioning forces greaterthan 1E + 05N.Figure 4.12 is the maximum displacements of the T-beams caused only by one patchloading. Those displacements were measured from the upward deflection caused by post-tensioning forces. See Figure 4.11. It is very obvious that the post-tensionging force canimprove the loading sharing capacity of the whole structure.T-beam 1T-beam 2Chapter 4. ANALYTICAL RESULTS^ 60PT*FlangeWebFlangeInternal ShearForceWeb4Figure 4.11: Upward Deflection. of the T - beam4^4^4 4^4^+Y PT---u) 25 ^szta)20 -Es.; m15a)2 -cwo ocn 5U_ ^C)1E+01^1E+02^1E+03^1E+04^1E+05^1E+06The Post-tensioning Force ( N )1 E+07Figure 4.12: Displacements of T-beam under One Patch LoadingChapter 4. ANALYTICAL RESULTS^ 614.2.3 Wheel loading analysisFour patches of vertical loading simulating a vehicle with four wheel were applied to thestructure. According to the different locations of the four wheels we have three cases tostudy.CASE 1The 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 inx-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 werethe 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 post-tensioning forces were increased, the contribution from beam 3 to the sharing of loadingincreased, producing a corresponding difference between the maximum bending stressesof 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 maximumrelative displacement between the beam 1 and beam 2 was very small due to their similarloading and boundary conditions. But the maximum absolute deformation of the joint2 in z-direction was very large. After the post-tensioning force got to 100,000 (N), thevalue of the maximum deformation of the joint element 2 in z-direction was negligible incomparison 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 wheelloading were then moved from x = 0.0 to x = 6000mm. The Figure 4.17 and Figure 4.18L PTPT/ Joint 1 Joint 2Chapter 4. ANALYTICAL RESULTS^ 621^2^3Figure 4.13: T-beam Structure CASE 1give us the non-linear analytical results from PTB program. When the vehicle movedto midspan, the displacement and the bending stress of the T-beam achieved the peakvalues.CASE 2The four patches of vehicle loading were acting at the edges of T-beam 1 and T-beam 2 asshown Figure 4.19. The deformations in joint elements were correspondingly more thanthose in CASE 1 at same level of the post-tensioning force. But the trend of the post-tensioning force effect was the same. When PT approached 10,000 (N), the tightness ofand load distribution in the structure was improved dramatically. The analytical resultsare shown from Figure 4.20 to Figure 4.22.The initial stiffnesses Ex and Ez of the springs have an influence on the magnitudeof the relative displacements. Figure 4.23 shows the effect on the maximum relative T-beam 1T-begm 2T-beam 3n insion toChapter 4. ANALYTICAL RESULTS^ 6302^.0 10 ^^8 ^E6° 6co4 -c—cowE) 2^ri 0 ^21E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07The Post-tensioning Force ( N )Figure 4.14: Bending Stress in T-beam Structure CASE 1L.s ens' ni^or -g- 25 ^cn^-^20a^E 15 -m 10 - 50tog (+00T-beam 1T-beam 2T-beam 31E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07The Post-tensioning Force ( N )Figure 4.15: Displacements in T-beam Structure CASE 1Chapter 4. ANALYTICAL RESULTS^ 6420 ^0 A a)15O1 0cnoa)> 5'o;crX 4+0021E+01^1E+02^1E+03^1E+04^1E+05The Post-tensioning Force ( N )1E+06^1E+07 Joint Element 1Joint  Elgment 2Figure 4.16: Relative Displacements in T-beam Structure CASE 1E ioa)I-"665 40)C 2'15CCOcd2oo^000 1000^2000^3000^4000^5000^6000Position of the Moving Loading(mm).... ... ...ndinCT-bepm 1T-Own 2T-b_eFn 3Figure 4.17: Bending Stress under Moving LoadingFigure 4.18: Displacements under. Moving LoadingPTLJoint 1^Joint 21^3E 25MSN.ro 2046 15Ea)MI 5ao0^ 1000^2000^3000Position of the Moving4000^5000Loading(mm)6000 T-beam 1T-beArri 2T-beam 3Chapter 4. ANALYTICAL RESULTS^ 65Figure 4.19: T-beam Structure CASE 2T-beam 1T-ber 2T-beam 3T-beam 1T-beam 2T-beam 3St4t6titibiliff'Chapter 4. ANALYTICAL RESULTS^ 66a.a10 ^a)U)8 ^6•c 42 2ei)^-"ti 00^a)co21E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07The Post-tensioning Force ( N )Figure 4.20: Bending Stress in T-beam Structure CASE 2^25 ^15a) 10co2 5 -U)^b 0^^g g ^2 11+00 1E+01^1E+02^1E+03^1E+04^1E+05^1E+06^1E+07The Post-tensioned Force ( N )Figure 4.21: Displacements in T-beam Structure CASE 21E+051E+01^1E+02^1E+03^1E+04The Post-tensioning Force ( N )•••••■■•Joint Element 1Joint Element 2t3^1E+06^1E+07Chapter 4. ANALYTICAL RESULTS^ 67Figure 4.22: Relative Displacements in T-beam Structure CASE 2vertical displacement if the horizontal spring stiffness Ex and the vertical spring stiffnessEz are changed from Ex = 16.9N/mm and Ez = 17.4N/mm to Ex = 169N/mm andEz = 174N/mm respectively. It is seen that these parameters should be known with someaccuracy if the maximum relative displacement is to be estimated at intermediate post-tensioning forces. However, the level of the post-tensioning force which makes the wholestructure work as a unit is less affected by the uncertainty in Ex and E. Determinationof Ex and Ex should include a larger experimental sample than used in this thesis.CASE 3In this case, the vehicle loadings were next the center lines of the T-beams as shownFigure 4.24. In addition to the maximum bending stress, maximum displacement in thebeam webs and the maximum deformation of the joint elements in the z-direction shownrespectively in Figure4.25 through Figure4.27, the program PTB was used to study theChapter 4. ANALYTICAL RESULTS^ 68EE 60c.) 50a) 40_c.• 3075._y 20 -• Ex = 169 N/mm76• 10 - Ez = 174 N/mma)CC1E+01^1E+02^1E+03^1E+04^1E+052^The Post-tensioning Force ( N )Ex= 16.9 N/mmE z = 17.4 N/mmJoint Element 1Softer SpringJoint Element 1Stiffer Spring1E+06^1E+07Figure 4.23: The Effect of the Stiffness of the Springconvergence 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, thegreater the coupling in the system of equations. Since the structure and the loadingcases were symmetric in x-direction, only odd number terms in the Fourier series wereincluded. From Figures 4.28 to 4.31, it can be seen that the convergence is quite goodwith only two Fourier terms (n = 1,3), being sufficient to obtain satisfactory answers.0-10—QCa)coco 6E 4C)C2Ccog Q2 it+oo 1E+01^1E+02^1E+03^1E+04The Post-tensioning Force ( N )1E+05^1E+06VS.T-beam 1—a-T-beam 2T-beam 3•••••...Chapter 4. ANALYTICAL RESULTS^69Figure 4.24: T-beam Structure CASE 3Figure 4.25: Bending Stress in T-beam Structure CASE 3T-beam 1—T•beem 2T-beam 325 ^E20,a^150co010a.U,0 ^1E+00 1E+01^1E+02^1E+03^1E+04The Post-tensioning Force ( N )1E+05^1E+06Joint Demerit 1—tt—Joint Element 2•Chapter 4. ANALYTICAL RESULTS^ 70E 35 ^O 300^0 25-• ._c020au, 15O 10CoT1 5ccii+ooFigure 4.26: Displacements in T-beam Structure CASE 31E+01^1E+02^1E+03^1E+04^1E+05^1E+06The Post-tensioning Force ( N )Figure 4.27: Relative Displacements in T-beam Structure CASE 3E 20asco 186".061a)14 -.0Ci 12.0 1Qi+cio 1E+01 1E+02^1E+03^1E+04 1 E+05 1E+06 1E+07Chapter 4. ANALYTICAL RESULTS^ 7108.5E 8a)co 7.5(/) 7ci)6.5co 65.5coa+oo„1E+01^1E+02^1E+03^1E+04^1E+05^1E+06The Post-tensioning Force ( N )1E+07Figure 4.28: Effect of the Fourier Terms on Bending StressThe Post-tensioning Force ( N )Figure 4.29: Effect of the Fourier Terms on Displacement1E+01 1E+04 1E+05^1E+06CC)E 2520515"e5a0 10Te.501E+0021E+02^1E+03The Post-tensioning Force ( N )Number of Terms713 051E+01^1E+02^1E+03^1E+04The Post-tensioning Force ( N )1 E+05 1 E+06Number of Terms7—a-135Chapter 4. ANALYTICAL RESULTS^ 72Figure 4.30: Effect of the Fourier Terms on Relative Displacement in Joint 1The Effect of the Fourier TermsFigure 4.31: Effect of the Fourier Terms on Relative Displacement in Joint 2Chapter 5CONCLUSIONSA non-linear analysis for the post-tensioned timber bridge has been developed and im-plemented in the computer program PTB. The analysis takes into account the frictionalcontact between the bridge beams. The relative movements between beams, the stressesand deformations of the structure can be obtained as a function of the post-tensioningforce 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 confirmthe 'same slip' assumption, and that reliability analyses be carried out to develop designcriteria. The structural analysis presented here can form the basis for such probabilisticinvestigation.73Bibliography[1] Csagoly, Paul F., and Taylor, Raymond J. 1980. " A Structural Wood System forHighway Bridges ", IABSE Periodica 4/1980 IABSE Proceedings pp. 35-80.[2] Research and Development Division Ontario Ministry of Transportation and Com-munications, 1979. " Transverse Post-Tensioning of Longitudinally Laminated Tim-ber Bridge Decks ", Research Report RR 220[3] Verna, J.R., Graham, J.F., Shannon, P.H. 1984. " Timber Bridges: Benefits andCosts ", J. Structural Engr., ASCE, 110(7): 1563-1571[4] Cheung, Y.K. 1968." Finite Strip Method in the Analysis of Elastic Plates withTwo Opposite Simply Supported Ends ", Proc., Inst. Civ. Engrs., 40(5):1-7[5] Foschi, Ricardo 0. 1982. " Structural Analysis of Wood Floor Systems ", StructuralDivision, Proceedings of ASCE, 108(ST7): 1557-1574[6] Thompson, E.G., Goodman, J.R., and Vanderbilt, M.D. 1975. " Finite ElementAnalysis 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: AProgram for the Analysis of Layered Wood Systems ", 7: 237-248[8] Taylor, Raymond J. and Csagoly, Paul F. 1980. " Rehabilitation of Wood HighwayBridges in Ontario Research Report Ministry of Transportation and Communi-cations74Bibliography^ 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 " Com-puters and Structures 407-426[10] Harik, Issam E. and Salamoun, Ghassan L. 1988. " The Analytical Strip Methodof 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 EngineeringDepartment, Colorado State Uni., Fort Collins, Colo., USA.[14] Cook, Robert D., Malkus,David S., Loesha,Michaele, 1989. " Concepts and Ap-plications of Finite Element Analysis " Third Edition University of Wisconsin-Madison[15] Dhatt, Gouri., Touzot, Gilbert., 1984. " The finite Element Method Displayed ",Chichester [ West Sussex ] New York: WileyAppendix ASHAPE FUNCTIONS AND DERIVATIVESPolynomial Shape Functions Vector {M0 } All components are zero except:mom= 42^-1M0(2 ) = —8 (e — e e4 e)M0( 9 ) = 2 — 2e + e)m0(10,1_42 +eM0(14) = E2^4E3^E411M0( 15) —8 ( —e e3^e5)Vector {M3} All components are zero except:M3 ( 3 ) = 1^ + 4e e376Appendix A. SHAPE FUNCTIONS AND DERIVATIVES^ 77M3(4 ) = 8^e + 63 - 6 4 )M3(7) = 1 — 26.2 + 4.41M3 (16) = —4(36 + 46 2 — 6 3 — 264 )M3(17)84^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) M1(k) —^; M2(k) — dMg(k);de dedM5(k).M3( k)  M6(k) = deM4(k) dWhere:^k = 1,2, ... ,19.Derivatives of the Displacement Functions aw(x,y)^N^ X = Etillo(4)1T {8,2 }( n7 r )COS( nr )az^n=1 L^LAppendix A. SHAPE FUNCTIONS AND DERIVATIVES(92 w 1x2 9)^E{iV10()}T{871}( 121-) 2sin( 717)n=1aW n(X' Y) =^{Mi(0}T{bn}( ) Sin ( n7rX(-/Y^n=1 s^LN52w5 (x ' Y) = E{M2()}T{Sn}(2) 2sin (nrx )s^L2^n=1192W(X , y)= NE{Mi()}T{Sn}(-2)( )c°s(nrx)OxOy ^L^LOu(x,y)^_^1r 1V1\lT fisnl(L )sn. (7717 )Ox n=1au(x,y) 014(0/T{871}(2)ws(nrxn=iav(x,y) N^nrx= E{m5(E)} {8n}( )cos(^ )Ox^n=1 L^Lav(x,y) —Oyd2 U (x)dx 2d2 W(x)dx 2NE{M6(0}T{8n}(-2^L)sin nrxsn=1NUn ( 117 ) 2 cos( n^ )n=1=^wn (7 )2sin(7 )n=1d2 V(x)dx 2d20(x)dx 2— 1%E7 Vn ( -7L-1- ) 2sin(-1-7)n=1 78Appendix BTHE PTB PROGRAMB.1 PTB Program SpecificationsThe non-linear analysis for the transverse post-tensioned bridge structure can be donewith the PTB computer program. The wheel loadings are idealized as patches of staticloading. PTB evaluates the response of the structure subjected to the transverse post-tensioning 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. PTBcan be run in an IBM 386 PC computer or a UNIX SUN work station. The programlimits 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)79Appendix 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' inthe 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 meansthat 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 expandedby changing the corresponding parameters in the PTB source code.B.2 Program StructureThe program PTB consists of one main program and a number of subroutines. Eachsubroutine performs its specific function. The flow chart of the Main Program of PTBis given below.PTB uses common blocks to convey data between the main program and subroutinesand between the subroutines as well.SUBROUTINE DATA This subroutine inputs data for the structure in free format through a data filePTB.DAT, which includes the size of the structure, the number of the terms in theCALCULATE IMBALANCE LOAD VECTORASSEMBLE TANGENTIAL GLOBAL STIFFNESS MATRIXCALCULATE INCREMENTAL DISPLACEMENTS OF THE STRUCTURE( STEP *)CHECK ALL THREE-DIMENSIONAL SPRINGSIF NO IF YESCHECK DISPLACEMENT CONVERGENCYAppendix B. THE PTB PROGRAM^ 81PTB PROGRAM FLOW CHART/ INPUT DATA /ASSEMBLE INCREMENTAL WHEEL LOADING VECTORAND GLOBAL STIFFNESS MATRIXBOUNDARY CONDITIONSSOLVE EQUATIONS CORRESPONDING n-th TERMOF FOURIER SERIES TO GET DISPLACEMENTS OF THE STRUCTUREUSING JACOBI ITERATION METHODCHECK ALL THREE-DIMENSIONAL SPRINGSIF YESIF NO/ OUTPUT RESULTS /Figure B.1: PTB Program Flow ChartAppendix B. THE PTB PROGRAM^ 82Fourier series, material properties, vehicle wheel loadings, transverse post-tensioningforces, number of three-dimensional springs in the joint elements, boundary conditionsetc. 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 six-point Gaussian quadrature to evaluate the contributions from the T-beam flanges to thestiffness 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 symmetricmatrix corresponding to each term of the Fourier series.SUBROUTINE STIF12 This subroutine creates the 12 x 12 stiffness matrix from one joint element with respectto each term of Fourier series. One joint element consists of 2 x N, three dimensionalsprings. N, has been input from the data file PTB.DAT.In this subroutine STIF12, subroutine BX, BY and BZ are called to evaluate thestiffness contributions from the x-, the y- and the z-direction springs according to thestatus of each spring respectively.SUBROUTINE CHECK This subroutine checks every spring in each joint element, according to the currentdisplacement and the constitutive relationship, to calculate the stiffness matrix and theinternal forces of the joint elements.SUBROUTINE LOAD This subroutine calculates the consistent load vector XF(n,IJ) according Equation 2.86for wheel loading and Equation 2.90 and Equation 2.91 for post-tensioning forces.In XF(n,IJ) where:Appendix B. THE PTB PROGRAM^ 83n = the n-th term in Fourier seriesIJ = 19 x NJTNJT = the number of the T-beam elements in the structureSUBROUTINE SLOAD According to the checking results obtained from SUBROUTINE CHECK, this sub-routine creates the internal force vector components from the 'yielded' springs.SUBROUTINE ASMBLYThis subroutine assembles the T-beam stiffness components, connector stiffness com-ponents and joint element stiffness components to form a global stiffness matrix by callingSUBROUTINE 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 bycolumn in the two dimensional matrix GSTIF(n,IJ). The components of the off-diagonalglobal stiffness matrix GK(n,m) are stored column by column in the two dimensionalmatrix ASTIF(nm,JK).SUBROUTINE GBC and SUBROUTINE ABC These two subroutines apply the boundary conditions to the stiffness matrices GS-TIF(n,JK), load vector XF(n,IJ) and stiffness ASTIF(nm,JK). To impose a specificboundary restrain to the structure, the corresponding columns and rows in the stiff-ness matrices are set to zero, and the corresponding components in the load vector isalso 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^ 84including the diagonal, is stored column by column in the GSTIF(n,JK). Cholesky decom-position 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 finalresults, according to Equation 2.106. Then the main program calls SUBROUTINE IT-ERATION to solve Equation 2.105. The final answer will be obtained when the iterationconverges.li{Auk}!,i^< Eil{Auk}i7n12 k =1,2, ... ,N and E is the error toleranceSUBROUTINE IMBALANCE This subroutine calculates the imbalance load for each iteration in Newton-Raphsonmethod. See Figure 2.8SUBROUTINE 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-directionand the z-directionIn all cases, the maximum response is obtained by comparing 19 values at equallyspaced locations along the longitudinal span.Appendix CPTB USER'S MANUALPTB: Post-Tensioned Timber Bridge Non-linear Analysis ProgramUser's ManualVersion 1.0by Cai ShenDepartment of Civil Engineering, UBCVancouver, B.C. Canada V6T 1W5November, 199185Appendix C. PTB USER'S MANUAL^ 86The computer program PTB Version 1.0, which is written in FORTRAN 77, is de-veloped 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 withthe input instructions given below. All data are to be entered in free format ( i.e. byproviding a space or comma between each data entry ). The capacity limitations of thePTB 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 deforma-tions and stresses in the structure and PTBB.OUT is a summary of the analysis resultsincluding the relative interface movements.Appendix C. PTB USER'S MANUAL^ 871. Enter: NTITLE NTITLE = problem title ( limited to 60 characters )2. Enter: NM, NJT, ISYMNM =maximum order of sine/cosine terms in the Fourier series. There can bemaximum of four terms. For x-direction symmetrical problem, NM can be upto 7 ( the four terms will be of 1, 3, 5 and 7 ). If NM = 5, there will be threeterms, 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)ISYM1^if the problem is x-direction symmetrical.0^if the problem is not x-direction symmetrical3. Enter: XL, SPJTXL = the span of the T-beamSPJT = T-beam spacing ( the width of the T-beam flange )4. Enter: BJT, HJT BJT = the width of the T-beam webHJT = the height of the T-beam webAppendix C. PTB USER'S MANUAL^ 885. Enter: EJT(I) EJT(I) = the elastic modulus of the i-th T-beam web6. Enter: REGREG = 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 flangeEMUX = the x-direction friction coefficient of the T-beam flangeEMUZ = the z-direction friction coefficient of the T-beam flange8. Enter: PEX, PEY, PG, PVXY, PVYXPEX = the elastic modulus of the T-beam flange in x-directionPEY = the elastic modulus of the T-beam flange in y-directionPG = the shear modulus of the T-beam flangePVXY = Poisson's ratio of T-beam flange, strain in x-direction while stress in y-directionPVYX = Poisson's ratio of T-beam flange, strain in y-direction while stress in x-direction9. Enter: SPC, XIN, CKPAL, CKPER, CKROT SPC = flange connector spacing along the T-beam spanXIN =distance between the end support and the first connector along the T-beamspanCKPAL, CKPER =Appendix C. PTB USER'S MANUAL^ 89connector load-slip moduli, respectively in the directions parallel and perpen-dicular to the T-beam ( unit = force/length )CKROT = connector rotation modulus10. Enter: TAOSX, TAOSZTAOSX = the shear strength of the T-beam flange in x-directionTAOSZ = the shear strength of the T-beam flange in z-direction11. Enter: NAN, N8 , MNL NAN =1^PTB will do non-linear analysis0^PTB will do linear analysis onlyN, =the number of the top three-dimensional springs in one joint element ( thetotal 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 thewheel 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 theI-th joint element.Appendix C. PTB USER'S MANUAL^ 9013. Enter: EY(I) EY(I) = the initial elastic modulus of the three dimensional spring in y-direction for theI-th joint element.14. Enter: EZ(I) EZ(I) = the initial elastic modulus of the three dimensional spring in z-direction for theI-th joint element.15. Enter: NWLD NWLD = the number of the wheel loadings; if NWLD = 0 goto step 1716. 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 appliedX1LD(I) =X2LD(I) =Y1LD(I) =Y2LD(I) =the coordinates of the I-th wheel loading patch, all coordinates are local T-beam coordinates ( I = 1, NWLD )17. Enter: NPTFNPTF =the number of the transversely post-tensioning forcesif NPTF = 0 goto step 20Appendix C. PTB USER'S MANUAL^ 9118. Enter: ES, AS, DIT, RFAC ES, AS = the elastic modulus and area of the post-tensioning cableDIT = the tooth distance of the threaded anchor head of post-tensioning tendonRFAC = factor to compensate the loss of post-tensioning force19. Enter: XPTF(I), PTF(I) XPTF(I) = the x-coordinate of the I-th pair of post-tensioning forcePTF(I) =the magnitude of the I-th pair of post-tensioning force( I = 1, NPTF )20 Enter: NBC NBC =the number of applied boundary conditionsif NBC = 0 skip step 2121. Enter: IBC(I,1), IBC(I,2) IBC(I,1) =the boundary conditions for I-th T-beamIBC(I,2) = the number of the constraint ( 1 to 19 ) ( I = 1, NBC )Appendix DPTB SOURCE CODE AND I/O FILEPTB.FORPTB.DATPTB.OUTPTBB.OUT92Appendix D. PTB SOURCE CODE AND I/O FILE^ 93D.1 PTB.FORCC^Post Tensioned Timber Bridge Non-Linear Analysis ProgramC Version 1.0 (Nov. 1991)^,kC^ Civil Engineering DepartmentC University of British ColumbiaC^ 2324 Main MallC Vancouver, Canada V6T-1W5CCCCC^C PROGRAM LIMITATIONSC--CC 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 )CC^ PROGRAM VARIABLES AND OPTIONSCC NTITLE = THE TITLE OF THE PROBLEMC NM^= MAXIMUM ORDER FOR THE SINE & COSINE SERIES. NM CAN BE UPC^TO 7 FOR SYMMETRIC PROBLEMS (ISYM = 1). NM CAN ONLY BE UPC TO 4 FOR. NON-SYMMETRIC PROBLEMS.C NJT^= NUMBER OF T-BEAMSC ISYM^= 1 FOR SYMMETRIC ABOUT X-DIC CASE; = 0 OTHERWISE.C XL^THE SPAN OF THE BRIDGEC SPJT^= THE FLANGE WIDTH OF THE T-BEAMC BJT^= THE WEB WIDTH OF THE T-BEAMC EMT^= THE WEB HEIGHT OF THE T-BEAMC EJT(I) = THE WEB ELASTIC MODULUS OF THE i-th T-BEAMC R.EG^= THE RATIO OF THE WEB ELASTIC MODULUS TO ITS SHEAR MODULUSC PTK^= THE FLANGE THICK OF THE T-BEAMC EMUX^= THE FRICTION COEFFICIENT OF THE FLANGE IN X-DIRECTIONC EMUZ^= THE FRICTION COEFFICIENT OF THE FLANGE IN Z-DIRECTIONC PEX^= THE ELASTIC MODULUS OF THE T-BEAM FLANGE IN X-DIRECTIONC PEZ^= THE ELASTIC MODULUS OF THE T-BEAM FLANGE IN Z-DIRECTIONC PG^= THE SHEAR MODULUS OF THE T-BEAM FLANGEC PVXY^= POISSON'S RATIO, STRAIN IN X-DIC WHILE STRESS IN Y-DICC PVYX^= POISSON'S RATIO, STRAIN IN Y-DIC WHILE STRESS IN X-DICC SPC^= NAIL SPACING ALONG THE LONGITUDINAL SPAN OF THE BRIDGEC XIN^= DISTANCE BETWEEN THE END SUPPORT AND THE FIRST NAILAppendix D. PTB SOURCE CODE AND I/O FILE^ 94C^ALONG THE LONGITUDINAL SPAN OF THE BRIDGEC CKPAL,CKPER= NAIL LOAD-SLIP MODULI, RESPECTIVELY, IN THE DIC.C^PARALLEL AND PERPENDICULAR TO THE WEB OF THE T-BEAMC CKROT = NAIL ROTATION MODULUSC TAOSX^= THE SHEAR STRENGTH OF THE FLANGE IN X-DIRECTIONC TAOSZ^= THE SHEAR STRENGTH OF THE FLANGE IN Z-DIRECTIONC NAN^= 1 - TO DO NON-LINEAR ANALYSISC^0 - TO DO LINEAR ANALYSIS ONLYC NP^= THE HALF NUMBER OF THREE DIMENSIONAL SPRINGS IN ONE JOINT ELE.C EX(I)^= THE INITIAL ELASTIC MODULUS OF THE X-DIRECTION SPRINGC EY(I)^= THE INITIAL ELASTIC MODULUS OF THE Y-DIRECTION SPRINGC EZ(I)^= THE INITIAL ELASTIC MODULUS OF THE Z-DIRECTION SPRINGC NWLD^= THE NUMBER OF THE WHEEL LOADINGSC JTLD(I) = THE NUMBER OF THE T-BEAM SUBJECTED BY THE WHEEL LOADINGCC X1LD(I)C X2LD(I) = THE COORDINATES OF THE i-th WHEEL LOADINGC Y1LD(I)C Y2LD(I)CC PPWLD(I) = MAGNITUDE OF THE WHEEL LOADINGC NPTF^= THE NUMBER OF THE TRANSVERSELY POST-TENSIONING FORCESC ES,AS^= THE MODULUS AND AREA OF THE POST-TENSIONING CABLEC DIT^= THE TOOTH DISTANCE OF THE CABLEC RFAC^= THEC XPTF(I) = THE X COORDINATE OF THE i-th POST-TENSIONING FORCEC PTF(I) = THE MAGNITUDE OF THE i-th POST-TENSIONING FORCEC 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).CIMPLICIT 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, TOLCOMMON /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,NA1COMMON /B05/ XL, SPJTCOMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RITCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY, PEW, PDGCOMMON /B08/ CKPAL, CKPER, CKROT, SPC,XINCOMMON /B09/ X1LD(MLD), X2LD(MLD), YiLD(MLD), Y2LD(MLD),1^PWLD(MLD).PPWLD(MLD), JTLD(MLD), NWLDCOMMON /B10/ XPTF(MPF), PTF(MPF),NPTFAppendix D. PTB SOURCE CODE AND I/O FILE^ 95COMMON /B11/ IBC(MBC,2), NBCCOMMON /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,MNLCOMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,F2,TOLL,NNACOMMON /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'*^ *',/T10,3'" PTB^ *',/T10.4'w^ *',/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**2PI4 PI"4TOL = .001CALL DATAWRITE(*,1)1 FORMAT(1X,'THE DATA HAVE BEEN INPUT')CALL GENMTXAppendix D. PTB SOURCE CODE AND I/O FILE^ 96NPT 6 * NPDC) 25 I = 1.NSTDO 25 J = 1,NPTICH(I,J) = 125 ICH1(I,J) = 0NNA 0DO 26 IK =1,NM.NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM =(1+IK)/ 2DO 26 I = 1,NEQSVEC2(IKM,I) = 0.D0SVEC(IKM,I) = 0.D0RFF(IKM,I) = 0.D026 FF(IKM,I) = 0.D0MMM = MNL + 1DO 350 LMN =1,MMMCIF ( LMN .EQ. 1 ) THENWRITE(*,*)WRITE(*,*)' APPLY POST-TENSIONING FORCE ONLY'END IFIF ( LMN .NE. 1 ) THENWRITE(*,*)WRITE(*,*)' APPLY WHEEL LOADING TO THE BRIDGE'END IFCDO 3 I = 1, NWLD3 PWLD(I) = (LMN-1)*PPWLD(I)/ (MMM 1)CALL LOAD(1.2)DO 30IK = 1,NM,NSTEPIKO =IKIF (NSTEP .EQ. 2) IK0 = (IK+1)/2DC) 30 I = 1,NEQ30 XFX(IKO,I) = XF(IKO,I)CALL ASMBLYWRITE(*,')8 FORMAT(BX,'THE GLOBAL STIFFNESS MATRIX HAS BEEN ASSMBLED')DO 9 I = 1, NWLD9 PWLD(I) = PPWLD(I)/ (MMM 1)WRITE(*,*)IF (LMN .EQ. 1) THENIF (NPTF .EQ.0) GOTO 350CALL LOAD(1,0)Appendix D. PTB SOURCE CODE AND I/O FILE^ 97END IFIF (LMN .GE. 2) THENIF (NWLD .EQ.0) GOTO 350CALL LOAD(0,2)END IFWRITE(*,io)10 FoRMAT(SX,'THE GLOBAL LOAD MATRIX HAS BEEN ASSMBLED')WRITE(*,*)DO 29 IK =1.NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM =(1-f-IK)/2DO 29 I = 1,NEQ29 XF(IKM,I) = XF(IKM,I) RFF(IKM,I)DO 60 IK = 1, NM, NSTEPDo 55 IN = 1,IK,NSTEPIKM = IKINM = INIF (NSTEP .EQ. 2) THENIKM = (IK + 1)/2INM = (IN + 1)/2END IFCIF (IN .NE. IK) THENIKO = (INM - 1) *NFT IKM^(INM 1)/2END IFIF ( IN .EQ. IK ) THENIKO = IKMCALL GBC(IKO)CALL DCOMP(IKO)Do 33 I = 1, NEQFF(IKO,I) = XF(IKM,I)33 CONTINUECALL SOLVE(IKO)DO 44 I = 1,NEQ44 SVEC(IKO,I) = XF(IKM,I)END IF55 CONTINUE60 CONTINUEC USING ITERATING METHOD TO SOLVE THE COUPLING EQUATIONCALL ITERATECIF ( NAN .EQ. o)G0To 315CCALL CHECK(0)WRITE(*,")^CHECK HAS BEEN FINISHED'WRITE( - ,*)Appendix D. PTB SOURCE CODE AND I/O FILE^ 98IF ( NNA .NE. 0 ) THENWRITE(',")'^NON-LINEAR ANALYSIS IS DOING 'WRITE(*,*)DO 240 IK = 1, NM, NSTEPIKM = IKIF (NSTEP .EQ. 2) THENIKM = (IK -I- 1)/2END IFDO 240 I = 1, NEQ240 SVEC1(IKM,I) = SVEC(IKM,I)C ^245 CALL IMBALANCE(1)CALL ASMBLYD() 249 IK = 1,NM,NSTEPDO 249 IN = 1,IK,NSTEFIKM = IKINM = INIF(NSTEP .EQ. 2) THENIKM = (IK-I-1)/2INM = (IN-f-1)/2END IFIF (IN .NE. IK) THENIKO = (INM - 1) 'NFT IKM INM" (INM 1)/2END IFIF ( IN .EQ. IK ) THENIKO = IKMCALL GBC(IKO)CALL DCOMP(IKO)DO 246 I =1,NEQFF(IKM,I) = XF(IKM,I)246 CONTINUECALL SOLVE(IKM)DO 297 I = 1,NEQ247 SVEC(IKM,I) = XF(IKM,I)END IF249 CONTINUECALL ITERATECNLF = 1DO 260 IK = 1, NM, NSTEPIKM = IKIF (NSTEP .EQ. 2) THENIKM = (IK + 1)/2END IFTP = 0.0TP1 = 0.0AMAX = 0.0DO 250 I = 1.NEQTP1 = TP1 + (SVEC1(IKM,I)"*2)+SVEC2(IKM,I)*x2TOLL = 0.01Appendix D. PTB SOURCE CODE AND I/O FILE^ 99250 TP = TP -1-(SVEC(IKM,I)"'2)WRITE(*,*)'TOLL = ',TOLLT1' = DSQRT(TP)TP1 = DSQRT(TP1)EP = TP/TP1WRITE(*,")'EP = NON LINEAR',EP,'IKM = ',IKMWRITE(*,*)'LMN = ',LMNWRITE(",*)'PWLD(I) = ',(LMN-1)*PPWLD(1)/ (MMM 1)NPA = 1IF ( EP .GT. TOLL) NPA = 0NLF = NLF NPADO 260 I = 1,NEQSVEC1(IKM,I) = SVEC1(IKM,I)^SVEC(IKM,I)260 CONTINUEDO 280 IK = 1, NM, NSTEPIKM = IKIF (NSTEP .EQ. 2) THENIKM = (IK -I- 1)/2END IFDO 280 I = 1, NEQ280 SVEC(IKM,I) = SVEC1(IKM,I)IF ( NLF .EQ. 1) GOTO 300WRITE(*,*)WRITE(*,*)'^NON-LINEAR ANALYSIS IS CONTINUING'WRITE(*,*)GOTO 245END IFCIF ( NNA .EQ. 0)THENWRITE(*,')'^NO NON-LINEAR SPRING OCCURS'WRITE(*,*)END IFC ^300 CONTINUECALL IMHALANCE(0)DO 310 I =1,NEQRFF(IKM,I) = XF(IKM,I)310 CONTINUE315 CONTINUEDO 320 IK =1,NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM =(1-1-IK)/2DO 320 I = 1,NEQ320 SVEC2(IKM,I) = SVEC(IKM,I)^SVEC2(IKM,I)350 CONTINUED() 380 IK =1,NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM =(1+IK)/2D() 380 I = 1,NEQ380 SVEC(IKM,I) = SVEC2(IKM,I)c^CALL CHECK(1)Appendix D. PTB SOURCE CODE AND I/O FILE^ 100CALL RESULTWRITE(",398)398 FORMAT(1X,'THE RESULTS HAVE BEEN OUTPUT')WRITE(*,*)CLOSE(1)CLOSE(2)CLOSE(3)STOPENDCC SUBROUTINE DATA READS INPUT DATA FROM DATA FILE.CSUBROUTINE DATAIMPLICIT 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,NA1COMMON /B05/ XL, SPJTCOMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RITCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY, PDV, PDGCOMMON /B08/ CKPAL, CKPEB., CKROT, SPC,XINCOMMON /B09/ X1LD(MLD), X2LD(MLD), Y1LD(MLD), Y2LD(MLD),1^PWLD(MLD),PPWLD(MLD), JTLD(MLD), NWLDCOMMON /B10/ XPTF(MPF), PTF(MPF),NPTFCOMMON /B11/ IBC(MBC,2), NBCCOMMON /B15/ EX(MSP ).EY(MSP ),EZ(MSP ),NAN,NP,MNLCOMMON /B17/ ICH(MSP.MNP),ICH1(MSP,MNP ),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNACOMMON /B18/ SIGMAY,WDX1(MSP ),WDX2(MSP),1^WDY1(MSP),WDY2(MSP),WDZ1(MSP),WDZ2(MSP)CHARACTER*72 NTITLEC INPUT DATA.READ (1,10) NTITLE10 FORMAT (A)READ (1,*) NM, NJT, ISYMNST=NJT-1READ (1,*) XL, SPJTREAD (1,*) BJT, HJTREAD (1,') (EJT(I),I=1,NJT), REGBETA = (1.0 - (0.63 BJT / HJT)) / 3.0AJT = BJT HJTRIY = BJT * (HJT*'3) / 12.0Appendix D. PTB SOURCE CODE AND I/O FILE^ 101RIZ = HJT (BJT"3) / 12.0RIT = BETA * HJT (BJT"3)DO 20 I = 1, NJTGJT(I) = EJT(I) / REG20 CONTINUEREAD (1,*) PTK,EMUX,EMUZREAD (1,") PEX, PEY, PG, PVXY, PVYXPKX = (PEX * PTK**3) / (12.0D0 (1.0D0 (PVXY*PVYX)))PKY = PKX PEY / PEXPKV = PVXY PKXPKG = PG * (PTK"3) / 12.0D0PDX = (PEX * PTK) / (1.0D0 - (PVXY*PVYX))PDY = PDX * PEY / PEXPDV = PVXY * PDXPDG = PG * PTKREAD (1,*) SPC, XIN,CKPAL, CKPER, CKROTCREAD(1,*) TAOSX,TAOSZREAD(1,*) NAN,NP,MNLIF ( NAN .EQ. 0)THENWRITE(*,*)' ** ONLY LINEAR ELASTIC ANALYSIS IS REQUIRED "'WRITE(*,*)END IFREAD(1,*) (EX(I), I=1,NST)READ(1,') (BY(I), I=1,NST)READ(1,*) (EZ(I), I=1,NST)c Jan. 20,1992c^ceexxx = ex(1)c^eeezaa^ez(1)c^ lcakREAD (1,") NWLDDO 150 I = 1, NWLDREAD (1,*) JTLD(I), X1LD(I), X2LD(I), Y1LD(I), Y2LD(I), PPWLD(I)150 CONTINUEREAD (1,*) NPTFIF (NPTF .NE. 0) THENREAD (1,')ES,AS,DIT,RFACDO 160 I = 1, NPTFREAD (1,1 XPTF(I), PTF(I)160 CONTINUEEND IFREAD (1,*) NBCIF (NBC .NE. 0) THEND() 170 I = 1, NBC170^READ (1,*) (IBC(I,J),J=1.2)180 CONTINUEAppendix D. PTB SOURCE CODE AND I/O FILE^ 102CDDXL = (5.0'9.8)*EMUX/EX(1)DDZL (5.0'9.8)*EMUZ/EZ(1)EXL = TAOSX *0.5"PTIC*XL/(NP-1) /DDXLEZL = TAOSZ '0.5"PTK'XL/(NP-1) /DDZLFXL = EXL DDXLFZL = EZL DDZLCIF (NPTF .NE. 0)THENAW = PTK*XL/NPTFFS = PTF(1)SIGMAY = FS/AWFX = SIGMAY EMUX *0.5"PTK*XL/(NP-1)FZ = SIGMAY EMUZ *0.5*PTK'XL/(NP-1)IF ( FX .GT. FXL ) FX = FXLIF ( FZ .GT. FZL ) FZ = FZLDO 190 I = 1,NSTEX(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) = EXLIF ( EZ(I) .GT. EZL ) EZ(I) = EELc Jan. 20.1992IF ( EX(I) .LT. EEEXXX ) EX(I) = EEEXXXIF ( EZ(I) .LT. EEEZZZ ) EZ(I) = EEEZZZc190 CONTINUEIF (EX(1) .NE. 0.0) DDX = FX/EX(1)IF (EZ(1) .NE. 0.0) DDZ = FZ/EZ(1)IF ( DDX .GT. DDXL ) DDX = DDXLIF ( DDZ .GT. DDZL ) DDZ = DDZLIF (EX(1) .EQ. 0.0) DDX = 0.0IF (EZ(1) .EQ. 0.0) DDZ = 0.0END IFCIF (NPTF .EQ. 0)THENEND IFAppendix D. PTB SOURCE CODE AND I/O FILE^ 103SIGMAY = 0.0FX = 0.0FZ = 0.0DO 191 I = 1,NSTEX(I) = 0.0EZ(I) = 0.0c Malcmagag.4.1(.****.W.K..IK*****WWW.71c.*.4.M.*..*MeMe.X.WA.c Jan. 20,1992IF ( EX(I) .LT. EEEXXX ) EX(I) EEEXXXIF ( EZ(I) .LT. EEEZZZ ) EZ(I) = EEEZZZ191 CONTINUEDDX = 0.0DDZ 0.0END IFC ^C TO CALCULATE THE NO. OF THE TURN OF THE POST TENSIONING CABLEC^IF (NPTF .NE. 0)THENTN = ( RFAC * FS * NJT SPJT ( 1/(ES'AS) 1/(PEY"AW) ))/DITEND IFC^C SET PARAMETERS FOR PROBLEM SIZE AND TYPE.C^ —NEg = NJT 19LHB = 19NSZ = 19 NEQC -NA]. = 144 * NSTIF (ISYM .EQ. 0) THENNSTEP = 1NFT = NMELSE IF (ISYM .EQ. 1) THENNSTEP = 2NFT = (NM + 1) / 2END IFC ECHOING INPUT DATA.WRITE(2,200) NTITLEWRITE(3,200) NTITLE200 FORMAT (/, lx, 26('''), ' FLOOR ANALYSIS PROGRAM *, 26('"'), //,' PROBLEM TITLE: ', A)WRITE(2,210) NJT, NFT, NMWRITE(3,210) NJT, NFT, NM210 FORMAT (/,' NUMBER OF FLOOR JOISTS1^' NUMBER OF FOURIER TERMS USEDAppendix D. PTB SOURCE CODE AND I/O FILE^ 1042^• MAX. ORDER OF FOURIER TERM^= 14 )WRITE(2,220)220 FORMAT (//,' PROPERTIES AND DIMENSIONS OF JOISTS')WRITE(2.230) XL, SPJT, BJT, HJT230 FORMAT (' JOIST SPAN =^E12.5, /,1 ' JOIST SPACING = E12.5, /,2 ' JOIST WIDTH = E12.5, /,3 JOIST DEPTH = E12.5 )Do 245 I = 1, NJTWRITE (2,240) I, EJT(I), GJT(I)240 FORMAT JOIST NO. =', 13, 2X, • EJT =', E12.5, 2X,1^GJT =', E12.5)245 CONTINUEWRITE(2,250)250 FORMAT (//,' PROPERTIES AND DIMENSIONS OF PLATE COVER')WRITE(2,260) PTK260 FORMAT (' COVER THICKNESS = E12.5)WRITE(2,270) PKX, PKY, PKV, PKG270 FORMAT (' KX = E12.5, 2X, ' KY = E12.5, /,1^' KV = E12.5, 2X, KG = E12.5 )WRITE(2,280) PDX, PDY, PDV, PDG280 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, CKROT350 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) SPC360 FORMAT SPACING BETWEEN CONNECTORS^E12.5, //)IF (NWLD .GT. 0) THENWRITE(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, NWLD390 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 IFIF (NPTF .GT. 0) THENWRITE (2,410)410 FoRMAT (' POST TENSIONING FORCES')WRITE (2,420)Appendix D. PTB SOURCE CODE AND I/O FILE^ 105420 FORMAT (' X-LOC', 6X, 'FORCE')DO 430 I = 1, NPTF430 WRITE (2,435) XPTF(I), PTF(I)935 FORMAT (2(2X,E12.5))END IFIF (NPTF .EQ. 0)THENWRITE(2,')WRITE(2,*)'NO POST TENSIONING FORCE'END IFWRITE (2,440) NBC440 FORMAT (//' NO. OF BOUNDARY CONDITIONS ', 14)IF (NBC .EQ. 0) GO TO 470DO 950 I = 1, NBC450 WRITE (2,960) (IBC(I,J)..1=1,2)460 FORMAT (' JOIST NO. = ', 13, 3X, ' B.C. CODE NO. = 13)970 CONTINUECNPP = NP*2WRITE(2,')WRITE(2,*)'THREE DIMENSIONAL SPRING BETWEEN TWO T. BEAMS',NPPWRITE(2,*)WRITE(2,1'STIFF. OF THE SPRING IN X-DIC.'DO 480 I = 1,NST980 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,NST500 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,NST520 WRITE (2,530) I,EZ(I)530 FORMAT ('^JOINT NO. = ', 13, 3X, ' EZ = E12.5)WRITE(2,540)EMUX590 FORMAT(1X,' FRICTION COEFFICIENT OF THE COVER (X-DIC.) =',E12.5)WRITE(2,543)EMUZ543 FORMAT(1X,' FRICTION COEFFICIENT OF THE COVER (Z-DIC.) =',E12.5)IF (NPTF .NE. 0)THENWRITE(2,700)ESWRITE(2,705)ASWRITE(2,710)DITWRITE(2,715)RFACWRITE(2,720)TN700 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)Appendix D. PTB SOURCE CODE AND I/O FILE^ 106END IFRETURNENDC SUBROUTINE GENTMX COMPUTES THE INTEGRALS REQUIRED IN THE DEVELOPMENTC OF THE ELEMENT STIFFNESS MATRIX NUMERICALLY USING SIX POINT GAUSSC QUADRATURE.C MCMCW*..,K*W**.747.(*.***..4..A(.1C...c*W7(..(*WWXDF11(.7..Y4.Wak...2K7(..**..,..f .*7K.Ma.”SUBROUTINE GENMTXIMPLICIT 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.93246951420315200ETA(2)^0.66120938646626500ETA(3)^0-23861918608331900ETA(4) -ETA(3)ETA(5) = -ETA(2)ETA(6) = -ETA(1)GWT(1) = 0.17132449237917000GWT(2) = 0.36076157304813900GWT(3) 0.46791393457269100GWT(4) = GWT(3)GWT(5) = GWT(2)GWT(6) = GWT(1)DO 10 I = 1, 6ETA2(I) = ETA(I) ETA(I)ETA3(I) = ETA2(I) ETA(I)ETA4(I) = ETA3(I) * ETA(I)ETAS(I) = ETA4(I) ETA(I)10 CONTINUEC MO AND M2 MATRICES.DO 20 I = 1, 6RM(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.0D0RM(I,15)^(-ETA2(I) ETA3(I)^ETA4(I)^ETA5(I)) / 8.000RM(I,9) = (ETA(I) - (2.000 * ETA3(I))^ETA5(I)) / 2.000Appendix D. PTB SOURCE CODE AND I/O FILE^ 107RM1(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.0D0RM1(I,15) = (-2.000 - (6.0130 * ETA(I)) + (12.000 ETA2(I))1^+ (20.000^ETA3(I))) / 8.0130RM1(I,9) = ((-12.0D0 ETA(I)) + (20.0130 ETA3(I))) / 2.00020 CONTINUECALL DMAT(1)DO 40 I = 1, 19DO 30 = 1, 19RMOO(I,J) = AA(I,J)30 CONTINUE40 CONTINUECALL DMAT(2)DO 60 I = I, 19DC) 50 J = 1, 19RM22(I,J) = AA(I,J)50 CONTINUE60 CONTINUECALL DMAT(3)D080 I = 1, 19DO 70 = 1, 19RMO2(I,J) = AA(I,J)70 CONTINUE80 CONTINUECALL DMAT(4)DO 100 I = 1, 19DO 90 .1 = 1, 19RM20(I,J) = AA(I,J)90 CONTINUE100 CONTINUECALL ZERO(RM, RM1)C M3 AND M6 MATRICES.DO 110 I = 1, 6RM(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.000RM(I,17) = (-ETA(I) - ETA2(I) ETA3(I) + ETA4(I)) / 8.000RM1(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)Appendix D. PTB SOURCE CODE AND I/O FILE^ 108RM1(I,6) = (-1.0D0^(2.01)0 ETA(I))^(3.0D0 ETA2(I))1^(4.0D0 ETA3(I))) / 8.01)0RM1(L19)^(-1.0D0 (2.0D0 ETA(I))^(3.0D0 ETA2(I))1^(4.0D0 * ETA3(I))) / 8.01,0110 CONTINUECALL DMAT(1)DO 130 I = 1, 19DO 120 J = 1,19RM33(I,J) = AA(I,J)120 CONTINUE130 CONTINUECALL DMAT(2)DO 150 I = 1, 19DO 140 J^1, 19RM66(I,J) = AA(I,J)140 CONTINUE150 CONTINUECALL DMAT(3)DO 170 I = 1, 19DO 160 J = 1, 19RM36(I,J) = AA(I,J)160 CONTINUE170 CONTINUECALL DMAT(4)DO 190 I = 1, 19DO 180 J = 1. 19RM63(I,J)^AA(I,J)180 CONTINUE190 CONTINUEC M5 AND M4 MATRICES.DO 200 I = 1.6RM(I,5) = RM(I,3)RM(I,3) = 0.0D0RM(I,18) = RM(I,16)RM(I,16)^0.0D0RM(I,8) = nm(I,T)Rm(i, ,r) = 0.0D0RM(I,6) = RM(I,4)RM(I,4) = 0.0D0RM(I,19) = RM(I,17)RM(I,17) = 0.0D0RM1(1,3) = RM1(I,5)RM1(I,5) = 0.0D0I1M1(1.16) = RMI(I,18)RMI(I,18) = 0.0D0= RM1(I,8)RM1(I,8) = 0.0D0RM1(I,4) = RM1(I,6)RM1(I,6) = 0.0D0RM1(I,17) = RM1(I,19)Appendix D. PTB SOURCE CODE AND I/O FILE^ 109itm1(I,19) = 0.000200 CONTINUECALL DMAT(2)DC) 220 I = 1, 19DC) 210 J = 1, 19RM44(I,J) = AA(I,J)210 CONTINUE220 CONTINUECALL DMAT(4)DC) 240 I = 1, 19DC) 230 .1 = 1, 19RM45(I,J) = AA(I,J)230 CONTINUE240 CONTINUECALL DMAT(3)DC) 260 I = 1, 19DC) 250 .1 = 1. 19RM54(I,J) = AA(I,J)250 CONTINUE260 CONTINUECALL DMAT(1)DO 280 I = 1, 19DO 270 .1 = 1, 19RM55(I,J) = AA(I,J)270 CONTINUE280 CONTINUECALL ZER.O(RM, RM1)C M1 MATRIX.DC) 290 I = 1, 6RM(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.000RM(I,15) = ((•2.000 * ETA(0) • (3.0D0 * ETA2(I))1^+ (4.000 * ETA3(I)) + (5.0D0 ETA4(I))) / 8.000RM(I,9) = (1.000 - (6.000 ETA2(I))1^+ (5.0D0 ETA4(I))) / 2.000290 CONTINUECALL DMAT(1)DO 310 I = 1, 19DO 300 .1 = 1, 19RM11(I,J) = AA(I,J)300 CONTINUE310 CONTINUERETURNENDAppendix D. PTB SOURCE CODE AND I/O FILE^ 110C SUBROUTINE ZERO CALLED BY GENMTX.CSUBROUTINE ZEROIMPLICIT DOUBLE PRECISION (A - H, 0 - Z)COMMON /B02/ AA(19,19), RM(6,19), RM1(6,19)DO 20 I = 1,6D 0 10 .1 = 1, 19RM(I,J)^0.0D0RM1(I,J) = 0.0D010 CONTINUE20 CONTINUERETURNENDC SUBROUTINE DMAT CALLED BY GENMTX.cSUBROUTINE 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, 19DO 10 .1 = 1, 19AA(I,J) = 0.0D010 CONTINUE20 CONTINUEDO 50 I = 1, 6DO 40 J = 1, 19Fl = RM(I,J)IF (II .EQ. 2 .OR. II .EQ. 4) F1 = R.M1(I,J)IF (F1 .NE. 0.0) THENDO 30 K = 1, 19F2 = RM1(I,K)IF (II .EQ. 1 .OR. II .EQ. 4) F2 = RM(I,K)IF (F2 .NE. 0.0) THENAA(J,K) = AA(J,K) + Fl * F2 GWT(I)END IF30^CONTINUEEND IF40 CONTINUE50 CONTINUERETURNENDC SUBROUTINE COVCON COMPUTES THE STIFFNESS MATRIX FOR THE PLATE COVERC AND CONNECTORS ASSOCIATED WITH A GIVEN T-BEAM STRIP.C^ mak M 1 11C.A.*.aRaleatalellt...”..ara, **.311■31C.*”....11(*.W.RSUBROUTINE COVCON (IN,IK)Appendix D. PTB SOURCE CODE AND I/O FILE^ 111IMPLICIT 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, TOLCOMMON /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,NA1COMMON /B05/ XL, SPJTCOMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RITCOMMON /B07/ PTK, PEX, PEy, PG, PVXy, PVYX, PKX, PKy, PKV, PKG,PDX, PDY, PDV, PDGCOMMON /B08/ CKPAL, CKPER, CKROT, SPC,XINCOMMON /B12/ ESTIF(19,19)COMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNADIMENSION AD(19,19)DO 20 I = 1, 19DO 10 J = 1, 19ESTIF(I,.I) = 0.010 CONTINUE20 CONTINUEC CONTRIBUTION FROM THE PLATE COVER.FAC = 0.0IF (IN .EQ. IK) FAC = PKX (IK**4) PI4 * Sp.TT / (4.0*XL* -3)DO 60 I = 1, 19DO 50 J = 1, 19AD(I,J) = ItMoo(I,J)50 CONTINUE60 CONTINUECALL ADD(AD, FAC)CFAC = 0.0IF (IN EQ. IK) FAC PKY * 4.0 XL / (SPJT**3)DO 100 I = 1, 19DO so = 1, 19AD(I,J) = RM22(I,J)90 CONTINUE100 CONTINUECALL ADD(AD, FAC)CFAC 0.0IF (IN Eq. IK) FAC = -PKV PI2 / (SPJT*XL)FAC1 FAC (IK"2)DO 140 I = 1, 19DO 130 .1 = 1, 19Appendix D. PTB SOURCE CODE AND I/O FILE^ 112AD(I,J) = KMO2(I,J)130 CONTINUE140 CONTINUECALL ADD(AD, FAC1)CFAC1 = FAC " (IN"2)DO 160 I = 1. 19DO 160 J = 1, 19AD(I,J) =150 CONTINUE160 CONTINUECALL ADD(AD, FAC1)CFAC = 0.0IF (IN .EQ. IK) FAC = PKG 4.0 " PI2 * (IK"2) / (SPJT"XL)DO 200 I = 1, 19DO 190 J = 1, 19AD(I,J)^RM11(I,J)190 CONTINUE200 CONTINUECALL ADD(AD, FAC)CPAC = 0.0IF (Ili .EQ. IK) FAC = PDX SPJT PI2 * (IK**2) / (4.0*XL)DO 240 I = 1. 19DO 230 J^1, 19AD(I,J) = RM33(I,J)230 CONTINUE240 CONTINUECALL ADD(AD, FAC)CFAC 0.0IF (IN .EQ. IK) FAC = PDY * XL / SPJTDO 280 I = 1, 19D0270 .1 = 1, 19AD(I,J) = KM66(I,J)270 CONTINUE280 CONTINUECALL ADD(AD, FAC)CFAC = 0.0IF (IN .EQ. IK) FAC = -PDV *PI 2.0FAC1 FAC " IKDO 320 I = 1, 19DO 310 J= 1, 19AD(I,J) = RM36(I,J)310 CONTINUE320 CONTINUECALL ADD(AD, FAC1)CFAC1 =PAC IND() 340 I = 1, 19Appendix D. PTB SOURCE CODE AND I/O FILE^ 113DO 330 J = 1, 19AD(I,J) = RM63(I,J)330 CONTINUE340 CONTINUECALL ADD(AD, FAC1)CFAC = 0.0IF (IN .EQ. IK) FAC = PDG. XL SP.1T / 4.0FAC1 = FAC 4.0 / (SPJT*'2)DO 380 I = 1, 19DO 370 J = 1,19AD(I,J) = RM44(I,J)370 CONTINUE380 CONTINUECALL ADD(AD, FAC1)CFAC1 = FAC 2.0 " IN PI / (SPJT*XL)DO 400 I = 1, 19DO 390 .1 = 1, 19AD(I,J) = RM45(I,J)390 CONTINUE400 CONTINUECALL ADD(AD, FAC1)CFAC1 = FAC * 2.0 " IK * PI / (SPJT"XL)DO 420 I = 1, 19DO 410 J = 1, 19AD(I,J) = 11M54(I,J)410 CONTINUE420 CONTINUECALL ADD(AD, FAC1)CFAC1 = FAC IK * IN PI2 / (XL*"2)DO 440 I = 1, 19DO 430 J = 1, 19AD(I,J) = RM55(I,J)430 CONTINUE440 CONTINUECALL ADD(AD, FACT)C CONTRIBUTION FROM THE COVER-TO-STIFFENER CONNECTORS.DC) 460 I = 1, 6DO 450 J = 1, 19AD(I,J) = 0.0450 CONTINUE460 CONTINUEAD(1,7) = 1.0AD(1,11) = -1.0AD(1,10) = -IK " PI (HJT t PTK) / (2.0"XL)AD(2,7) = 1.0AD(2,11) = -1.0Appendix D. PTB SOURCE CODE AND I/O FILE^ 114AD(2,10) = -IN * PI * (HIT + PTK) (2.0 *XL)AD(3,8) = 1.0AD(3,12) = -1.0AD(3,9) = -PTK / (2.0*SPJT)AD(3,13) = -HIT / (2.0*SPJT)AD(4,8) = 1.0AD(4,12) = -1.0AD(4,9) = -PTK / (2.0"SPIT)AD(4,13) = -HIT / (2.0*SPJT)AD(5,9) = 1.0 / SPITAD(5,13) = -1.0 / SPITAD(6,9) = 1.0 / SPITAD(6,13) = -1.0 / SPJTCIF (IN NE. IK) GO TO 510SPAR = XL / (2.0*SPC)SPER = XL / (2.0'SPC)SROT = SPERGO TO 520510 SPAR = 0.0SPER = 0.0SROT = 0.0520 SPAR = SPAR * CKPALSPER = SPER CKPERSROT = SROT CKROTCDO 540 I = 1, 19DC) 630 J = 1, 19ESTIF(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 CONTINUE540 CONTINUERETURNENDCC FUNCTION Z IS USED TO CALCULATE Ix( n,m,i)C alea....acakag 3r.FUNCTION Z(N, M, Z1, ZO, XL, NX, PI)IMPLICIT DOUBLE PRECISION (A - H, 0 -Z)IF (N .EQ. M) GO TO 20S1 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 10Z = XL "(D1- D2) /PIRETURN10 Z = XL * (Dl + D2 )/PIAppendix D. PTB SOURCE CODE AND I/O FILE^ 115RETURN20 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 30Z = XL * (DI - D2)/PIRETURN30 Z = XL * (DI + D2)/PIRETURNENDCC SUBROUTINE ADD, WHICH IS CALLED BY SUBROUTINE COVCON, BUILDS THEC STIFFNESS AND MASS MATRICES ASSOCIATED WITH THE PLATE COVERC AND THE CONNECTORS FOR A GIVEN T-BEAM STRIP.CSUBROUTINE ADD (AD, FAC)IMPLICIT DOUBLE PRECISION (A - H, 0 Z)COMMON /B12/ ESTIF(19,19)DIMENSION AD(19,19)IF (FAC .EQ. 0.0D0) RETURNDO 20 I^1, 19DC) 10 J^1, 19ESTIF(I,J) = ESTIF(I,J) + (FAC AD(I,J))10 CONTINUE20 CONTINUERETURNENDCC SUBROUTINE ASMBLY ASSEMBLES THE DIAGONAL STIFFNESS SUB-MATRICESC 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 ASMBLYIMPLICIT 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, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /B06/ EJT(MJT), GJT(MJT), BJT, OUT, AJT, RIY, RIZ, RITCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY, PDV, PDGCOMMON /B11/ IBC(MBC,2), NBCCOMMON /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),Appendix D. PTB SOURCE CODE AND I/O FILE^ 1161^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,MNLCOMMON /B17/ ICH(MSP.MNP),ICH1(MSP.MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL.NNADo 60 IK = 1, NM, NSTEPDO 55 IN = 1,IK,NSTEPIKM = IKINM = INIF (NSTEP .EQ. 2) THENIKM = (IK + 1)/2INM = (IN + 1)/2END IFCIF (IN .NE. IK) THENIKO = (INM 1) 'NFT IKM 'NM' (INM + 1)/2DO 5 I = 1, NA15 ASTIF(IKO,I)= 0.0DO 20 IE = 1,NJTIJ^(IE 1)*19IF (IE .LE. NST ) THENCALL STIF12(IE,IN.IK)C^FBC39.F NOV. 8,1991CIF ( NBC .NE. 0 ) THENCALL ABC(IKO,IE)END IFIJ = ( IE - 1 ) * 12Do 17 J= 1,12DO 17 K = 1,12JK = ( IJ J 1 ) * 12 KASTIF(IKO,JK) = ASTIF(IKO,JK) EK12(J,K)17 CONTINUECEND IF20 CONTINUEEND IFIF ( IN .EQ. IK ) THENIKO = IKMFAC2 = (IK**2) PI2 / (2.0D0 * XL)FAC4 = (IK**4) PI4 / (2.oD0 * (XL"3))Appendix D. PTB SOURCE CODE AND I/O FILE^ 117DO 25 I = 1, NSZGsTiF(IK0,i) = o.oDo25 CONTINUEC CONTRIBUTION FROM THE PLATE COVER AND CONNECTORS.CALL COVCON (IN,IK)Do 50 IE = 1, NJTIJ = (IE -^19DO 90 = 1, 19DO 30 K = 1,JK = LHB (IJ + K - 1) + J K + 1GSTIF(IKO,JK) = GSTIF(IKO,JK) + ESTIF(J,K)30^CONTINUE40 CONTINUEC CONTRIBUTION FROM THE JOISTS.JK = LHB (IJ + 9) + 1GSTIF(IKO,JK) = GSTIF(IKO,JK) + (EJT(IE) * RIY * FAC4)JK = LHB * (IJ + 10) + 1GSTIF(IKO,JK) = GSTIF(IKO,JK) + (EJT(IE) * AJT * FAC2)JK = LHB (IJ + 11) + 1GSTIF(IKO,JK) = GSTIF(IKO,JK) + (EJT(IE) " RIZ * FAC4)JK LHB^+ 12) + 1GSTIF(IKO,JK) = GSTIF(IKO,JK) + (GJT(IE) RIT FAC2 /(sPJT*'2))IF (IE .LE. NST ) THENCALL STIF12(IE,IN,IK)IJ=(IE-1)*19+13DO 45 .1=1,12DO 45 K=i,JJK=LHB -(IJ+K-1)+J-K+1GSTIF(IKO,JK) = GSTIF(IKO,JK) + EK12(J,K)45 CONTINUEEND IF50 CONTINUEEND IF55 CONTINUE60 CONTINUERETURNENDCC THIS SUBROUTINE GENERATES THE STIFNESS MATRICES OF THE JOINTC ELEMENT BETWEEN THE T-BEAM STRIP.CSUBROUTINE 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)Appendix D. PTB SOURCE CODE AND I/O FILE^ 118PARAMETER (MSP = 9, MFO = 6, MNP 378)COMMON /B00/ PI, PI2. PI4, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ.NA1COMMON /B05/ XL, SPJTCOMMON 0307/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY. PDV, PDGCOMMON /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,MNLCOMMON /1317/ ICH(MSP,MNP).ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNAM=IKN=INDO 20 J=1,8BX1(1,J)=0.0BX2(1,J)=0.0BY1(1,J)=0.0BY2(1,J)=0.0BZ1(1,J)=0.013Z2(1,3)=0.0BX1T(J,1)=0.0BX2T(J,1)=0.013Y1T(3,1)=0.0BY2T(J,1)=0.0BZ1T(J,1)=0.0BZ2T(.7,1)=0.020 CONTINUECALL BX(N,M,IE)CALL BY(N,M,IE)CALL BZ(N,M,IE)DO 60 1=1,8DO 60 J=1,8EK(I,J)=SXK(I,J)+SYK(I,J)+SZK(I,J)60 CONTINUECALL EXPRETURNENDCC THE CONTRIBUTION OF THE SPRING IN X DIC. CALLED BY STIF12C 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)Appendix D. PTB SOURCE CODE AND I/O FILE^ 119PARAMETER (MFT = 4, MLD = 12, MBC 20, MPF = 20)PARAMETER (MSP = 9, MFO = 6, MNP = 378)COMMON /1100/ PI, PI2, PI4, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY, PDV, PDGCOMMON /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,MNLCOMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNADO 10 1=1,8DO 10 .1=1,810 SXK(I,J)=0.0DO 100 I=1,NPDX=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 + 1CBX1T(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 .^+ 2CBX2T(1,1)=-PTK/2.'(1,PTI/XL)-DCOS(N*PI*XI/XL)*EX(IE)1^*ICH(IE,KX2)Appendix D. PTB SOURCE CODE AND I/O FILE^ 120BX2T(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,8DC) 20 3=1,820 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,8DO 30 3=1,830 SXK(II,J)=SXK(11,J)+XK(I1,J)100 CONTINUERETURNENDCC THE CONTRIBUTION OF THE SPRING IN Y DIC. CALLED BY STIF12CSUBROUTINE 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, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG.1^PDX, PDY, PDV, PDGCOMMON /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,MNLCOMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNADO 10 1=1,8DO 10 J=1,810 SYK(I,J)=0.0DO 100 I=1,NPDX=XL/(NP-1)IF ( I.EQ.1 .OR. I.EQ.NP ) DX=0.5'XL/(NP-1)XI=(I-1).XL/(NP-1)Appendix D. PTB SOURCE CODE AND I/O FILE^ 121BY1(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 + 3IF(ICH(IE,KY1) .EQ. 1) EKY1^1.00*EY(IE)IF(ICH(IE.KY1) .EQ. -1) EKY1 = 1.00 1.00E+10IF(ICH(IE,KY1) .EQ. 0 ) EKY1 = 0.00CBY1T(2,1)=PTK/(2."SPJT)*DSIN(N'PI*XI/XL)1^"EKY1BY1T(4,1 )=-DSIN(N'PI*XI/ XL)1^*EKY1BY1T(6,1)=-PTK/(2."SPJT)*DSIN(N'PrXI/XL)1^*EKY1B Y1 T(8,1 )=DSIN(NWI'XI/ XL)1^*EKY1BY2(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 +4IF(ICH(IE,KY2) .EQ. 1) EKY2 = 1.00*EY(IE)IF(ICH(IE,KY2) .EQ. -1) EKY2 = 1.00 1.00E+10IF(ICH(IE,KY2) .EQ. 0 ) EKY2 = 0.00CBY2T(2,1)=-PTK/(2.*SPJT)*DSIN(WPI*XI/XL)1^*EKY2BY2T(4,1)=-DSIN(N*PrXI/XL)1^*EKY2BY2T(6,1)=PTK/(2.'SPJT)*DSIN(WPI*XI/XL)1^*EKY2BY2T(8,1)=DSIN(N*PrXI/XL)1^*EKY2CALL DGMULT (BY1T.BY1,YK,8,1,8,8,1,8)DO 20 11=1,8DC) 20 J=1,820 SYK(I1,J)=SYK(I1,J)+YK(I1,J)CALL DGMULT (BY2T,BY2,YK,8,1,8,8,1,8)DO 30 11=1,8DO 30 J=1.830 SYK(I1,J)=SYK(II,J)+YK(I1,J)100 CONTINUERETURNENDcC THE CONTRIBUTION OF THE SPRING IN Z DIC. CALLED BY STIF12Appendix D. PTB SOURCE CODE AND I/O FILE^ 122CSUBROUTINE 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, TOLCOMMON /B04/ NM, NSTEP, NET, NJT,NST, NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /B07/ PTK, PEX. PEY, PG, PVXY. PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY, PDV, PDGCOMMON /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,MNLCOMMON /B17/ ICH(MSP,MNP).ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNADO 10 1=1,8DO 10 3=1,810 SZK(I,J)=0.0DO 100 I=1,NPDX=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 + 5CBZ1T(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,8DO 20 J=1,820 SZK(I1,J)=SZK(I1,J)+ZK(11,J)100 CONTINUECDO 200 I=1,NPDX=XL/(NP-1)IF ( I.EQ.1 .OR. I.EQ.NP ) DX=0.5*XL/(NP-1)XI=(I-1)*XL/(NP-1)Appendix D. PTB SOURCE CODE AND I/O FILE^ 123B 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- 6CB 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,8DO 21 J=1,821 SZK(I1,..1)=SZK(I1,J)+ZK(11,J)200 CONTINUECRETURNENDC ***************************************************************C THIS SUBROUTINE EXPANDS EK(8,8) TO EK12(12,12) CALLED BY STIF12CSUBROUTINE EXPIMPLICIT 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,12DO 10 J=1,1210 EK12(I,J)=0.0DO 100 J=1,3DO 20 1=1,320 EK12(I,J)=EK(I,J)EK12(5,J)=EK(4,J)DC) 40 1=7,915=1-290 EK12(I,J)=EK(I5,J)EK12(11,J)=EK(8,J)100 CONTINUEDC) 120 1=1.3120 EK12(I,5)=EK(I,9)DO 130 1=5Appendix D. PTB SOURCE CODE AND I/O FILE^ 124130 EK12(5,5)=EK(4,4)DO 140 1=7,915=1-2140 EK12(I.5)=EK(15,4)EK12(11,5)=EK(8,4)DO 300 .1=7,9.75=J-2DO 220 1=1,3220 EK12(I,J)=EK(I,J5)EK12(5,J)=EK(4,..75)DO 240 1=7,915=1-2240 EK12(I,J)=EK(I5,J5)EK12(11,J)=EK(8,J5)300 CONTINUEDO 320 1=1,3320 EK12(I,11)=EK(I,8)EK12(5,11)=EK(4,8)DO 340 1=7,915=1-2340 EK12(I,11)=EK(I5,8)EK12(11,11)=EK(8,8)RETURNENDcC TO MULTIPLY A REC. M BY N MATRIX BY ANOTHER RECT.C N BY L MATRIX.cSUBROUTINE 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, LDO 20 I = 1, MC(I,J)^0.D0DO 10 K = 1, N10^C(I,J) = C(I,J)^A(I,K) B(K,J)20 CONTINUERETURNENDcC SUBROUTINE LOAD COMPUTES THE GLOBAL LOAD VECTOR.cSUBROUTINE 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, TOLAppendix D. PTB SOURCE CODE AND I/O FILE^ 125COMMON /B03/ ETA(6), GWT(6)COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1COMMON /B06/ XL, SPJTCOMMON /B09/ X1LD(MLD), X2LD(MLD), YlLD(MLD), Y2LD(MLD),1^PWLD(MLD),PPWLD(MLD). JTLD(MLD), NWLDCOMMON /B10/ XPTF(MPF), PTF(MPF),NPTFCOMMON /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.NNADIMENSION W(19)DO 20 IK = 1, NM, NSTEPIKO = IKIF (NSTEP .EQ. 2) IKO (IK +1) / 2DO 10 I = 1, NEQXF(IKO,I) = 0.0D010 CONTINUE20 CONTINUEC CONTRIBUTION FROM TRANSVERSE LOADINGIF (IIII .EQ. 2) THENIF (NWLD^0) THENDo 70 IL = 1, NWLDIE = JTLD(IL)1.1 = (IE - 1)^19DO 60 IK = 1, NM, NSTEPIK0 = IKIF (NSTEP .EQ. 2) IKO = (IK +1) / 2PINL = rI IK / XLDO 30 I = 1, 19W(I) = 0.0D030^CONTINUEXIX = 2.0D0 Y1LD(IL) / SPJTXIY 2.0D0 Y2LD(IL) / SPJTDO 40 I = 1, 6XI = XIX + (XIY - XIX) * (1.0D0 + ETA(I)) / 2.0D0XI2 = XI**2XI3 = XI2 * XIXI4 = XI3 XIXI5 = XI4 XIW(1) = W(1) + (XI2 (5.0D0 XI3 4.0D0)- (XI4 / 2.0D0) + (3.0130 XI5 / 4.0D0)) GWT(I)W(2) = W(2) + (XI2 - X13 - XI4 + XI5) GWT(I) / 8.0D0W(9) = W(9) + (XI - (2.0D0 * XI3) + XI5) * GWT(I) /2.0D0W(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)Appendix D. PTB SOURCE CODE AND I/O FILE^ 126W(15) = W(15) + (-XI2 - X13 + XI4 + XIS) G^T(I) / 8.0D040^CONTINUEY1 = 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) FACCW(2) = W(2) * FACCW(9) = W(9) FACCW(10) = W(10) FACCW(14) = W(14) FACCW(15) = W(15) * FACCDO 50 I = 1, 19= XF(IKO,I.1+I) + ((PWLD(IL) W(I)))50^CONTINUE60 CONTINUE70 CONTINUEEND IFEND IFC CONTRIBUTION FROM POST TENSIONING FORCESIF (III .EQ. 1) THENIF (NPTF .NE. 0) THENDO 90 I = 1. NPTFDO 80 IK = 1, NM, NSTEPIKO = IKIF (NSTEP -EQ. 2) IKO = (IK +1) / 2PINL = PI IK / XLV5 = DSIN(PINL*XPTF(I))IJ = 5XF(IKO,IJ) = XF(IKO,IJ) -f (PTF(I) V5)V18 = DSIN(PINL*XPTF(I))IJ^((NJT - 1) * 19) + 18XF(IKO,IJ) = XF(IKO,IJ) - (PTF(I) V18)80 CONTINUE90 CONTINUEEND IFEND IFRETURNENDCC THIS SUBROUTINE CALCULATES THE EQUIVALENT LOADS CONTRIBUTEDC FROM THE NON-LINEAR SPRINGSCSUBROUTINE SLOAD(IE,IK)Appendix D. PTB SOURCE CODE AND I/O FILE^ 127IMPLICIT 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, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ.NA1COMMON /B05/ XL, MITCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,^1^PDX, PDY, PDV, PDGCOMMON /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,MNLCOMMON /B11/ ICH(MSP,MNP ).ICH1(MSP.MNP ).SL(12),1^DDX.DDZ,FX,FZ,TOLL,NNADO 10 I = 1,1210^SL(I) = 0.0DX = XL/(NP - 1)DO 100 IP = 1,NPIF (IP .EQ. 1 .OR. IP .EQ. NP) Dx=0.5-xL/(Np-1)XI = (IP - 1)*XL/(NP - 1)KX1 = (IP - 1)*6 + 1KX2 = (IP - 1).6 + 2KZ1 = (IP - 1)-6 +KZ2 = (IP - 1)*6 + 6SL(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 CONTINUERETURNENDAppendix D. PTB SOURCE CODE AND I/O FILE^ 128C SUBROUTINE GBC APPLIES THE BOUNDARY CONDITIONS TO THE STIFFNESS MATRIXC GSTIF(IKO,JK) AND LOAD VECTOR XF(IKO,JK).CSUBROUTINE 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,NA1COMMON /B11/ IBC(MBC,2), NBCCOMMON /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) THENDO 30 I = 1, NBCNE = IBC(I,1)IDOF = IBC(I,2)M = (NE - 1) 19 IDOFC= 1MM = M - 1IF (M .GE. 19) .11 = M 18DO 10 .1 = .11, MMJK = (LHB - 1)^- 1) MGSTIF(IKO,JK) = 0.010 CONTINUEMI = M 1Ms = M 4- 18IF (M8 .GT. NEQ) MS = NEQDo 20 .1 = MI, M8JK = (LHB - 1) (M 1) JGSTIF(IKO,JK) = 0.020 CONTINUEJK = (LHB - 1) * (M 1) MGSTIF(IKO,JK) = 1.0XF(IKO,M) = 0.030 CONTINUEEND IFRETURNENDCC SUBROUTINE ABC APPLIES THE BOUNDARY CONDITIONS TO THE STIFFNESS MATRIXC ASTIF(IKO,JK) .CSUBROUTINE 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,NA1Appendix D. PTB SOURCE CODE AND I/O FILE^ 129COMMON /B11/ IBC(MBC,2), NBCCOMMON /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) THENDO 20 I = 1, NBCNE = IBC(I,1)IDOF^IBC(I,2)CIF (NE .EQ. IE) THENIF (IDOF .GE. 14) THENJ= IDOF - 13DO 10 K = 1,1210 EK12(J,K) = 0.0K = IDOF - 13DO 14 J^1,1214 EK12(J.K) = 0.0END IFEND IFIF ( NE .EQ. IE--1 ) THENIF ( IDOF .LE. 6 ) THEN.I.1 = IDOF + 6DO 16 KK = 1,1216 EK12(JJ,KK) = 0.0KK = IDOF + 6DO 18 JJ = 1,1218 EK12(JJ,KK) = 0.0END IFEND IF20 CONTINUEEND IFRETURNENDcC SUBROUTINE DCOMP DECOMPOSES THE STIFFNESS DIAGONAL SUB-MATRICESAppendix D. PTB SOURCE CODE AND I/O FILE^ 130C USING THE CHOLESKY METHOD.CSUBROUTINE 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,NA1COMMON /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 1TEMP = GSTIF(IKO,1)TEMP DSQRT(TEMP)GSTIF(IK0,1) = TEMPDO 10 I = 2, LHBGSTIF(IKO,I) = GSTIF(IKO,I) / TEMP10 CONTINUEDO so .1 = 2, NEQJi = .1 - 1IJD = LHB J KBSUM =_- GSTIF(IKO,IJD)Ko = 1IF (.1 .GT. LHB) K0 = J • KDDO 20 K = KO. J1.1K = KB K .1 - KBTEMP = GSTIF(IKO,JK)SUM = SUM • TEMP -* 220 CONTINUEGSTIF(IKO,IJD) DSQRT(SUM)DO 50 I = 1, KBII = .1 4- IK0= 1IF (II .GT. LED) KC) II - KBSUM = GSTIF(IKO,IJD + I)IF (I .EQ. KB) GO TO 40DC) 30 K = KC), .11JK = KB K J - KBIK = KB K II - KBTEMP = GSTIF(IK0,JK)SUM = SUM - GSTIF(IKoJK) TEMP30 CONTINUE40^GSTIF(IKO,IJD + I) = SUM / GSTIF(IKO,IJD)50 CONTINUE60 CONTINUERETURNENDcAppendix D. PTB SOURCE CODE AND I/O FILE^ 131C SUBROUTINE SOLVE SOLVES THE SYSTEM OF EQUATIONS USING THEC DECOMPOSED MATRIX OBTAINED FROM SUBROUTINE DCOMP.CSUBROUTINE 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,NA1COMMON /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 - 1TEMP = GSTIF(IKO,1)XF(IKO,1) = XF(IKO,1) / TEMPDO 20 I = 2, NEQII = I - 1KO = 1IF (I .GT. LHB) KO = I - KBSUM = XF(IKO,I)II = LHB I - KBDO 10 K = KO, IiIK = KB K + I - KBTEMP = GSTIF(IKO,IK)SUM = SUM - (TEMP XF(IKO,K))10 CONTINUEXF(IKO,I) = SUM / GSTIF(IKO,II)20 CONTINUEC BACKWARD SUBSTITUTION.Ni = NEQ - 1LB LHB NEQ - KBTEMP = GSTIF(IKO,LB)XF(IKO,NEQ) = XF(IKO,NEQ) / TEMPDO 40 I = 1, NiIi = NEQ . I + 1NI = NEQ IKO = NEQIF (I .GT. KB) KO = NI + KBSUM = XF(IKO,NI)II = LHB * NI • KBDO 30 K = Il, KOIK = KB NI K - KBTEMP = GSTIF(IKO,IK)SUM = SUM - (TEMP XF(IKO,K))30 CONTINUEXF(IKO,NI) = SUM / GSTIF(IKO,II)40 CONTINUEAppendix D. PTB SOURCE CODE AND I/O FILE^ 132RETURNENDCC ITERATING METHOD IS USED IN THIS SUBROUTINE TO SOLVEC THE COUPLING EQUATIONS.CSUBROUTINE ITERATEIMPLICIT 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, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT, NST.NEQ, LHB, NSZ.NA1COMMON /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.NNADIMENSION X(MEQ),XX(MEQ)IF (NM .EQ. 1) GO TO 123oC ^1090 NFLAG = 1DO 1220 IK = 1, NM, NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM = (IK + 1) / 2DO 1100 I = 1. NEQ1100 XX(I) = FF(IKM,I)DO 1170 IN = 1, NM, NSTEPIF (IN .EQ. IK) GO TO 1170INM = INIF (NSTEP .EQ. 2) INM = (IN + 1) / 2IF (IN .GT. IK) GO TO 1110IKO = (INM 1) * NFT + IKM INM (INM + 1) / 2GO TO 11201110 IKO = (IKM 1)* NFT + INM IKM (IKM + 1) / 2C1120^DO 100 IE = 1, NSTIJ1 = ( IE - 1 ) 19 + 141.12 = ( IE - 1 ) 19 + 14 + 11= ( IE - 1 ) * 144DO Bo I = 1.11,IJ2TEMP = 0.0CHI = I - (IE - 1 ) *19 .13DO 60 = DI, 1.12CHII = J - (IE - 1 ) *19 -13TEMPI = SVEC(INM,J)CHICHI = IJJ + (CHI-1)*12 + CHIIAppendix D. PTB SOURCE CODE AND I/O FILE^ 13360^TEMP = TEMP + ASTIF(IKO,CHICHI) * TEMPI80^XX(I) XX(I) - TEMP100^CONTINUEC1170 CONTINUEC^DO 1190 I = 1, NEQ1190 XF(IKM,I) = XX(I)CALL SOLVE(IKM)DO 1195 I = 1,NEQ1195 X(I) = XF(IKM,I)CTEMP = 0.0TEMPI. = 0.0DO 1200 I = 1, NEQTEMP = TEMP + (X(I)**2)TEMPI = TEMPI (SVEC(IKM,I)**2)1200 CONTINUETEMP = DSQRT(TEMP)TEMPI = DSQRT(TEMP1)EPS = DABS((TEMP TEMPI)/TEMPI)CGIKM = 0.0AIKO = 0.0DO 1204 NANA = 1,NSZ1204 GIKM = GIKM (GSTIF(IKM,NANA))**2DO 1206 MAMA = 1,NA11206 AIKO = AIKO (ASTIF(IKO,MAMA))**2• WRITE(",'")'EPS = SOLVE',EPS,'IKM^',IKM• WRITE(*,*)'GIKM = ',GIKM,' IKM = ',IKM• WRITE(',*)'AIKO = ',AIKO,' IKO = ',IKOC ^NAC = 1IF (EPS .GT. TOL) NAC = 0NFLAG = NFLAG NACDO 1210 I = 1, NEQ1210 SVEC(IKM,I) = X(I)1220 CONTINUEIF (NFLAG .EQ. 1) GO TO 1230GO TO 10901230 CONTINUERETURNENDC ***************************************************************C THIS SUBROUTINE CHECKS THE DEFORMATIONS OF THE SPRING ELEMENTSC 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.Appendix D. PTB SOURCE CODE AND I/O FILE^ 134SUBROUTINE CHECK(LMNLMN)C^ LMNLMN 1 WRITE = 0 DON'T WRITEIMPLICIT 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, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,1^PDX, PDY, PDV, PDGCOMMON /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,MNLCOMMON /B17/ ICH(MSP,MNP ),ICH1(MSP,MNP ),SL(12).1^DDX,DDZ,FX.FZ,TOLL,NNACOMMON /B18/ SIGMAY,WDX1(MSP),WDX2(MSP),1^WDY1(MSP ),WDY2(MSP ),WDZ1(MSP ),WDZ2(MSP )DIMENSION S(12)CDO 5 IK =1,NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM =(1.+IK)/2DO 5 I = 1,NEQ5 SVEC(IKM,I) = SVEC(IKM,I) + SVEC2(IKM,I)DO 170 IE = 1,NSTWDX1(IE) = 0.D0WDX2(IE) = 0.D0WDY1(IE) = 0.D0WDY2(IE) 0.D0WDZ1(IE) = 0.D0WDZ2(IE) = 0.D0DO 150 IP = 1,NPXI = (IP - 1)`XL/(NP - 1)DX1 = 0.0DX2 = 0.0DY1 = 0.0DY2 = 0.0DZ1 = 0.0DZ2 = 0.0DO 130 IK = 1,NM,NSTEPIKO = IKIF (NSTEP .EQ. 2) IKO = (IK + 1)/2DO 20 I =14.19.1 = (IE 1)"19 +I09 +^- .11) = ZZX^+^- dI) = TZX+^- dI) = ZAXC^9.(T • dI) = TAXZ^- dI) = ZXXI ÷^- dI) = TXX1RNIS,K0 011^ 0zz + zza = zzatz + TZCI = TZCICA z.ka = zAaTA. + TA.411 = TA.(11CX zxa = zxaIX + 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 0a.(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=TXC(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)s9 + I = XI + 61. =r9'1 = I OF OCIHRNISNO0 OZ(r`OXI)OHAS = (X)SCT - I = XTIM (VI (INV acrop 3-Janos Ead. •a xrpuaddv9ETAppendix D. PTB SOURCE CODE AND I/O FILE^ 136IF (DY1 .LE. 0.0 .AND. ICH(IE,KY1) .NE. 0) THENIF (ICH(IE.KX1) .EQ. 1) THENIF (DADS(DX1) .GE. DDX) THENICH(IE.KX1) = 0IF(DX1 .GT. 0.0) ICH1(IE,KX1) = 1IF(DX1 .LT. 0.0) ICH1(IE,KX1) = -1NNA NNA 1END IFEND IFICH(IE,KY1) = -1END IFCIF (DY1 .GT. 0.0) THENICH(IE,KX1) = 0ICH1(IE,KX1) = 0NNA = NNA 1ICH(IE,KY1) = 0NNA = NNA 1END IFCIF (DY2 .LE. 0.0 .AND. ICH(IE,KY2) .NE. 0) THENIF (ICH(IE,KX2) .EQ. 1) THENIF (DABS(DX2) .GE. DDX) THENICH(IE,KX2) = 0IF(DX2 .GT. 0.0) ICH1(IE,KX2) = 1IF(DX2 .LT. 0.0) ICH1(IE,KX2) = -1NNA = NNA 1END IFEND IFICH(IE,KY2) = -1END IFCIF (DY2 .GT. 0.0) THENICH(IE,KX2) = 0ICH1(IE,KX2) = 0NNA NNA 1ICH(IE,KY2) = 0NNA NNA 1END IFCC IN Z - DIRECTION " Z1 "IF (DY1 .LE. 0.0 .AND. ICH(IE,KY1) .NE. 0) THENIF (ICH(IE,KZ1) .EQ. 1) THENIF (DABS(DZ1) .GE. DDZ) THENAppendix D. PTB SOURCE CODE AND I/O FILE^ 137ICH(IE,KZ1) = 0IF(DZ1 .GT. 0.0) ICH1(IE,KZ1) = 1IF(DZ1 .LT. 0.0) ICH1(IE,KZ1) = -1NNA NNA 1END IFEND IFEND IFC IN Z - DIRECTION ** Z2 "IF (DY1 .LE. 0.0 .AND. ICH(IE,KY1) .NE. 0) THENIF (ICH(IE,KZ2) .EQ. 1) THENIF (DABS(DZ2) .GE. DDZ) THENICH(IE,KZ2)^0IF(DZ2 .GT. 0.0) ICH1(IE,KZ2)^1IF(DZ2 .LT. 0.0) ICH1(IE,KZ2) = -1NNA = NNA 1END IFEND IFEND IFIF (DY1 .GT. 0.0) THENICH(IE,KZ1) = 0ICH1(IE,KZ1) 0NNA = NNA 1END IFIF (DY2 .GT. 0.0) THENICH(IE,KZ2) = 0ICH1(IE,KZ2) = 0NNA = NNA 1END IFIF (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 CONTINUE170 CONTINUEDO 200 IK =1,NM,NSTEPIKM = IKIF (NSTEF .EQ. 2) IKM =(1+IK)/2DO 200 I = 1,NEQ200 SVEC(IKM,I) = SVEC(IKM,I) SVEC2(IKM,I)RETURNENDAppendix D. PTB SOURCE CODE AND I/O FILE^ 138c seg. ae.**.*.**214....*^ ale***xe3k..mx*M..*..C THIS SUBROUTINE CALCULATES THE IMBALANCE LOAD CAUSEDC BY THE NON LINEAR DEFORMATION OF THE SPRING ELEMENTSCSUBROUTINE IMBALANCEIMPLICIT 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, TOLCOMMON /B04/ NM, NSTEP, NFT, NJT,NST, NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /BOG/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, luxCOMMON /B07/ PTK, PEX, PEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG.^1^PDX, PDY. PDV, PDGCOMMON /B11/ IBC(MBC,2), NBCCOMMON /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,MNLCOMMON /B17/ ICH(MSP.MNP).ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNADIMENSION SG(19),SG1(19)CALL CHECK(0)DO 5 IK =1,NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM =(1-1-IK)/2DO 5 I = 1,NEQ5 SVEC(IKM,I) = SVEC(IKM,I) SVEC2(IKM,I)DO 20 IK = 1,NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM = (IK + 1)/2DO 10 I = 1,NEQ10 XF(IKM,I) XFX(IKM,I)DO 20 IE = 1,NSTCALL SLOAD(IE,IK)IJ = (IE^1)=19 -I- 13DO 20 I = 1,1220^XF(IKM,IJ-f-I) =^ SL(I)CDO 300 IK = 1,NM,NSTEPIKM = IKAppendix D. PTB SOURCE CODE AND I/O FILE^ 139IF (NSTEP .EQ. 2) THENIKM = (IK+1)12END IFDO 300 IN = 1.NM,NSTEPIF (IN -EQ. IK) THENFAC2 (IK"2) * PI2 (2.0D0 * XL)FAC4 = (IK"4) * PI4 / (2.0D0 (XL**3))DO 200 IE = 1,NJTCALL 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 +112 -= IX 19J = 0DO 120 I = 11,12J = 3+1120 SG(J) = SVEC(IKM,I)CALL DGMATV (ESTIF,SG,SG1,19,19,19)= 0DC) 140 I = 11.12J=3+1140 XF(IKM,I) = XF(IKM,I) SG1(J)IF (IE .LE. NST) THENCALL STIF12(IE,IK,IN)K = 0= (IE 1)*19 + 1412 = Ii +11DC) 160 I = 11,12K = K+1160 SG(K) = SVEC(IKM,I)CALL DGMATV (EK12,SG,SG1,12,12,12)K = 0DO 180 I = 11,12K K+1180 XF(IKM,I) = XF(IKM,I) SG1(K)END IF200 CONTINUEEND IFIF (IN .NE.IK) THENINM = INIF (NSTEP .EQ. 2) INM = (IN +1)/2DO 290 IE = 1,NJTIF (IE .LE. NST) THENCALL STIF12(IE,IK,IN)Appendix D. PTB SOURCE CODE AND I/O FILE^ 140K = 0= (IE 1)*19 + 1412 =^4-11DO 260 I = 11,12K = K+1260 SG(K) = SVEC(IKM,I)CALL DGMATV (EK12,SG,SG1,12,12,12)K = 0DO 280 I = 11,12K = K+1280 XF(INM,I) = XF(INM,I) - SG1(K)END IF290 CONTINUEEND IF300 CONTINUEC APPLIES THE BOUNDARY CONDITIONS TO THE LOAD VECTOR XF(IKO,M)IF (NBC.NE.0) THENDO 330 I = 1, NBCNE = IBC(I,i)IDOF IBC(I,2)M = (NE - 1) * 19 4- IDOFCDO 330 IK = 1.NM,NSTEPIKM = IKIF (NSTEP .EQ. 2) IKM = (IK+1)/2XF(IKM,M) = 0.0330 CONTINUEEND IFDO 365 IK =1,NM,NSTEP1KM = IKIF (NSTEP .EQ. 2) IKM =(14-IK)/2DC) 365 I = 1,NEQ365 SVEC(IKM,I) = SVEC(IKM,I) SVEC2(IKM,I)RETURNENDcC THIS SUBROUTINE IS TO MULTIPLY A RECTANGULAR M BY NC MATRIX BY A VECTOR OF LENGTH N .CSUBROUTINE 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, MW(I) = 0.D0DC) 10 J = 1, N10 W(I) = W(I)^A(I,J) V(J)RETURNAppendix D. PTB SOURCE CODE AND I/O FILE^ 141ENDCC SUBROUTINE RESULT EVALUATES THE DEFLECTIONS AND STRESSES IN THEC 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.AWSUBROUTINE RESULTIMPLICIT 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, TOLCOMMON /B03/ ETA(6), GWT(6)COMMON /B04/ NM, NSTEP, NFT, NJT, NST,NEQ, LHB, NSZ,NA1COMMON /B05/ XL, SPJTCOMMON /B06/ EJT(MJT), GJT(MJT), BJT, HJT, AJT, RIY, RIZ, RITCOMMON /B07/ PTK, PEX, FEY, PG, PVXY, PVYX, PKX, PKY, PKV, PKG,^1^PDX, PDY, FDV, PDGCOMMON /B08/ CKPAL, CKPER, CKROT, SPC,XINCOMMON /B10/ XPTF(MPF), PTF(MPF),NPTFCOMMON /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,MNLCOMMON /B17/ ICH(MSP,MNP),ICH1(MSP,MNP),SL(12),1^DDX,DDZ,FX,FZ,TOLL,NNACOMMON /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)CWRITE(3,*)WRITE(3,1'WRITE(3,1' SUMMARY OUTPUT'WRITE(3,1'WRITE(3,1WRITE(3,*)'WDY IS MAX. RELATIVE DISPL. IN Y-DIRECTION'WRITE(3,1'WDZ IS MAX. RELATIVE DISPL. IN Z-DIRECTION'WRITE(3,1WRITE(3,1161)PTF(1)1161 FORMAT(1X,' POST TENSION FORCE = E12.5)WRITE(3,*)WRITE(31' SIGMAY = ',SIGMAYWRITE(3,2)2 FORMAT(1X ,'IE ',4X,'SIGMAY',9X ,'WDY1',9X ,'WDY2'.9X ,'WD Z1 ',11X ,1^'WDZ2')WRITE(3,*)DO 3 IE =1,NST3^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)Appendix D. PTB SOURCE CODE AND I/O FILE^ 142WRITE(3,*)CIF (NAN .EQ.1) THENWRITE(2,1END IFIF (NAN .EQ. 0) THENWRITE(2,*)WRITE(2,1' ELASTIC LINEAR ANALYSISEND IFW.TMAX = 0.0D0SJMAX = 0.0D0WPMAX = 0.0D0SPMAX = 0.0D0ALPHA = 1.2D0C LOOP OVER EACH T-BEAM STRIPDO 170 IE = 1, NJTCC EVALUATING DEFLECTIONS AND STRESSES IN THE JOISTS.DC) 10 I = 1, 19WH(I) = 0.0D0WS(I) = 0.0D0ST(I) 0.0D010 CONTINUEDO 90 IK = 1, NM, NSTEPPIN = PI * IKPINL = PIN / XLIKO IKIF (NSTEP .EQ. 2) IKO = (IK + 1) / 2DO 20 I = 1, 19J= ((IE - 1) * 19) + IF(I) = SVEC(IKO,J)20 CONTINUECFACS = 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 * FACWDO 30 I = 1, 19XIL DSIN(I*PIN/20)WB(I) = WB(I) + (F(10) * XIL)WS(I) = WS(I) + (FACW XIL)ST(I) ST(I) + (FACS * XIL)30 CONTINUE90 CONTINUECAppendix D. PTB SOURCE CODE AND I/O FILE^ 143WJD = 0.0D0WJS 0.0D0WJT = 0.ODOST.1 = 0.0130DO 50 I = 1, 19IF (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 WTIF (DABS(ST(I)) .GE. DABS(STJ)) STJ = ST(I)50 CONTINUEC EVALUATING DEFLECTIONS AND STRESSES IN THE PLATE COVER.DO 60 I = 1, isS1T(I) = 0.0D0S1M(I) = 0.0D0CS1B(I)^0.0D0S2T(I) = 0.0D0S2M(I) = 0.ODOS2B(I) = 0.0D0C60 CONTINUEDO 62 I = 1,3862 WP (I) = 0.ODOCEPC = PEY / (1.0D0 (PVXY * PVYX))DO 130 IK = 1, NM, NSTEPPIN = PIIKPINL = PIN / XLIKO = IKIF (NSTEP .EQ. 2) IKO = (IK 1) / 2DO 70 I = 1, 19= ((IE - 1) 19) + IF(I) = SVEC(IKO,J)70 CONTINUECHT = PTK / 2.0D0FAC1T = (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))^aana,Lism^oft^  0(mx.(Wa) + (I)aM = (I)dM^bit(OZ/NId*(6T -I))NIS = 'IIX41Z 4 OZ = I iZT Oa('IIX. (Oa) + (I)dM = (I)dM(OZ/IslId*I)NIS = 'IIX61'T = I ZZT OQ^ 0ariNumoo('lix azona) + (i)azs = (Dais(ant vgzova) + (I)1lis = (DiAzzs('IIX * Izova) + (i),Lzs = Wizs('IIX amva) + Wats = (Daisvgmva) + (1)1-kas = Winas(gix^lova) + (i)Its = (I)Sts(oz/ma-Dmsa =6T^= I OZT OCTOda * aZOIrel = aZ0V3Oda * L1tZOV3 = TAIZOV3Oda I.Z0V3 = IZOVeIOda * EITOArd =adGI * TAIT01,4 = movaoda^= 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)((" )3 oac•e) + ( (Oa^oao•g))((z—zras) / Ix)) (Wet^XAAd)- (Irds / (((603 * oagz•o) - ((Oa oagz•o)- ((Wei oaog•T) + ((g)a oag•T-))) = Izova0((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) = aTOV3T11,1 o/I GNV acroo aounos ELM 'a xnDua ddvttTAppendix D. PTB SOURCE CODE AND I/O FILE^ 145WPC = 0.0D0Si = 0.0D0S2 = 0.0D0S3 = 0.0D0CS4 = 0.0D0S5 = 0.0D0S6 = 0.0D0CDC) 140 I = 1, 19IF (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 CONTINUEDO 144 I = 1,38144 IF (DABS(WP(I)) .GT. DABS(WPC)) WPC = WP(I)CC PRINTING OUT RESULTS.CWRITE (2,155) IE155 FORMAT (/, ' *" ELEMENT ', I2, ' **')WRITE (2,160) WJB, WJS, WJT. STJ, WPC, S1,S3,54,S6160 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 )CWRITE (3,155) IEWRITE (3,161) WJT, STJ, WPC161 FORMAT (' MAX. JOIST DEFLECTION ( TOTAL) = E12.5, /,3^' MAX. JOIST BENDING STRESS^= E12.5, /,4^' MAX. COVER DEFLECTION^= E12.5, )CIF (DABS(STJ) .GE. DABS(SJMAX)) SJMAX = STJIF (DABS(WJT) .GE. DABS(WJMAX)) WJMAX = WJTIF (DABS(WPC) .GE. DABS(WPMAX)) WPMAX = WPCIF (DABS(S1) .GE. DABS(SPMAX)) SPMAX = 51IF (DABS(S3) .GE. DABS(SPMAX)) SPMAX = S3IF (DABS(54) .GE. DABS(SPMAX)) SPMAX = S4IF (DABS(S6) .GE. DABS(SPMAX)) SPMAX = S6Appendix D. PTB SOURCE CODE AND I/O FILE^ 146170 CONTINUEC SUMMARY OF THE ANALYSIS.WRITE (2,180) WJMAX, SJMAX, WPMAX, SPMAX180 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.5WRITE(2,90)90 FORMAT (/,1X, 50('''))RETURNENDENDAppendix D. PTB SOURCE CODE AND I/O FILE^ 147D.2 PTB.DAT"`' PTB^PTB^PTB " PTB7 3 1 FOURIER NO. JOIST SYMMETRIC12000.0 1000.0^ XL, SPJT100.0 500.0 BJT,HJT14000.00 19000.00 14000.0017.0^ REG100.0 .5 .65^ PTK, EMUX,EMUZ19000.00 1400.00 1000.00 0.02 0.20^PEX,PEY,PG,PVXY,PVYX10.0 0.0 1.0E06 1.0E06 1.0E06 SPC,XIN,CKPAL,CKPER,CKROT3.0,2.0^ TAOSX,TAOSZ1 21 1 NAN, NP, MNL16.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)41 4500.00 5000.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 151 7000.00 7500.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 152 4500.00 5000.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 152 7000.00 7500.0 200.0 500.0 1.00-01 DX1 = DY1 = DZO = 15102.0E5 78.53 10.0 2.0 ES AS DIP RFAC600.0 1.00e+051800.0 1.00e+053000.0 1.00e+054200.0 1.00e+055400.0 1.00e+056600.0 1.00e+057800.0 1.00e-1-059000.0 1.00e+0510200.0 1.00e+0511400.0 1.00e+050Appendix D. PTB SOURCE CODE AND I/O FILE^ 148D.3 PTB.OUT""*****"*"**"*""*""" FLOOR ANALYSIS PROGRAM **"*****""**""*"***""PROBLEM TITLE: " PTB PTB " PTB^PTBNUMBER OF FLOOR JOISTS^= 3NUMBER OF FOURIER TERMS USED^4MAX. ORDER OF FOURIER TERM^--:-- 7PROPERTIES AND DIMENSIONS OF JOISTSJOIST SPAN^= 0.12000E+05JOIST SPACING = 0.10000E+04JOIST WIDTH^= 0.10000E+03JOIST DEPTH^= 0.50000E +03JOIST NO. = 1 EJT = 0.19000E+05 GJT = 0.82353E+03JOIST NO. = 2 EJT = 0.14000E+05 GJT = 0.82353E+03JOIST NO. = 3 EJT = 0.14000E+05 GJT = 0.82353E+03PROPERTIES AND DIMENSIONS OF PLATE COVERCOVER THICKNESS = 0.10000E+03KX = 0.11714E+10 KY = 0.11714E+09KV = 0.23427E+08 KG = 0.83333E+08DX = 0.14056E+07 DY = 0.14056E+06DV = 0.28112E+05 DG^0.10000E+06PROPERTIES FOR CONNECTORSSTIFFNESS PARALLEL TO JOIST^= 0.10000E+07STIFFNESS PERPENDICULAR TO JOIST = 0.10000E+07ROTATIONAL STIFFNESS FLANGE/JOIST = 0.10000E+07SPACING BETWEEN CONNECTORS = 0.10000E+02APPLIED TRANSVERSE LOADINGJOIST^X1^X2^Y1^Y2^LOAD1^0.95000E+04 0.50000E+04 0.20000E+03 0.50000E+03 0.10000E+001^0.70000E+04 0.75000E+04 0.20000E+03 0.50000E+03 0.10000E+002^0.45000E+09 0.50000E+04 0.20000E+03 0.50000E+03 0.10000E+002^0.70000E+04 0.75000E+04 0.20000E+03 0.50000E+03 0.10000E+00POST TENSIONING FORCESX-LOC^FORCE0.60000E+03 0.10000E+060.18000E+04 0.10000E+060.30000E+04 0.10000E+060.42000E+04 0.10000E+06Appendix D. PTB SOURCE CODE AND I/O FILE^ 1490.54000E+04 0.10000E+060.66000E+04 0.10000E+060.78000E+04 0.10000E+060.90000E+04 0.10000E+060.10200E+05 0.10000E+060.11400E+05 0.10000E+06NO. OF BOUNDARY CONDITIONS = 0THREE DIMENSIONAL SPRING BETWEEN TWO T- BEAMS 42STIFF. OF THE SPRING IN X-DIC.JOINT NO. = 1 EX = 0.86224E+05JOINT NO. = 2 EX = 0.86224E+05STIFF. OF THE SPRING IN Y-DIC.JOINT NO. = 1 EY = 0.10000E+11JOINT NO. = 2 EY = 0.10000E+11STIFF. OF THE SPRING IN Z-DIC.JOINT NO. = 1 EZ = 0.88776E+05JOINT NO.. 2 EZ = 0.88776E+05FRICTION COEFFICIENT OF THE COVER (X-DIC.) = 0.50000E+00FRICTION COEFFICIENT OF THE COVER (Z-DIC.) = 0.65000E+00STIFFNESS OF THE POSTENSIONING CABLE = 0.20000E+06THE AREA OF THE POSTENSIONING CABLE = 0.78530E+02THE DISTANCE BETWEEN THE CABLE TEETH = 0.10000E+02THE RELAX FACTOR OF THE CABLE = 0.20000E+01THE TURN NC). REQUIRED = 0.41773E+01ELEMENT 1 "*MAX. JOIST DEFLECTION (BENDING) = 0.12558E+02MAX. JOIST DEFLECTION ( SHEAR ) = 0.12855E+01MAX. JOIST DEFLECTION ( TOTAL ) = 0.13843E+02MAX. JOIST BENDING STRESS^0.56304E+01MAX. COVER DEFLECTION 0.13665E+02MAX. COVER STRESS (NODE 1 TOP) = 0.10653E+01MAX. COVER STRESS (NODE 1 BOTTOM)= 0.99706E+00MAX. COVER STRESS (NODE 2 TOP) = 0.96243E+00MAX. COVER STRESS (NODE 2 BOTTOM)= 0.87784E+00ELEMENT 2 **MAX. JOIST DEFLECTION (BENDING) = 0.13133E+02MAX. JOIST DEFLECTION ( SHEAR ) = 0.15197E+01MAX. JOIST DEFLECTION ( TOTAL ) = 0.14652E+02MAX. JOIST BENDING STRESS^= 0.64713E+01MAX. COVER DEFLECTION^= 0.13660E+02MAX. COVER STRESS (NODE 1 TOP) = 0.17037E+01MAX. COVER STRESS (NODE 1 BOTTOM). 0.80633E+00MAX. COVER STRESS (NODE 2 TOP) = 0.11909E+01MAX. COVER STRESS (NODE 2 BOTTOM)= 0.69862E+00Appendix D. PTB SOURCE CODE AND I/O FILE^ 150"" ELEMENT 3 ""MAX. JOIST DEFLECTION (BENDING) = 0.80606E+01MAX. JOIST DEFLECTION ( SHEAR ) = 0.81426E+00MAX. JOIST DEFLECTION ( TOTAL ) = 0.88749E+01MAX. JOIST BENDING STRESS^= 0.35843E+01MAX. COVER DEFLECTION^= 0.11385E+02MAX. COVER STRESS (NODE 1 TOP) = 0.15598E+01MAX. COVER STRESS (NODE 1 BOTTOM)= 0.77819E+00MAX. COVER STRESS (NODE 2 TOP) = 0.94627E+00MAX. COVER STRESS (NODE 2 BOTTOM)_ 0.91086E+00"*".."*"*"*" SUMMARY OF FLOOR ANALYSIS **"'"*"'"`"MAX. JOIST DEFLECTION^= 0.14652E+02MAX. JOIST STRESS^= 0.64713E+01MAX. COVER DEFLECTION^= 0.13665E+02MAX. COVER STRESS BETWEEN JOISTS = 0.17037E+01Appendix D. PTB SOURCE CODE AND I/O FILE^ 151D.4 PTBB.OUT*****"*"`"*"*"*"*"" FLOOR ANALYSIS PROGRAM """*""*"*`*****"""PROBLEM TITLE: " PTB^PTB "" PTB^PTBNUMBER OF FLOOR JOISTS^= 3NUMBER OF FOURIER TERMS USED = 4MAX. ORDER OF FOURIER TERM = 7SUMMARY OUTPUTWDY IS MAX. RELATIVE DISPL. IN Y-DIRECTIONWDZ IS MAX. RELATIVE DISPL. IN Z-DIRECTIONPOST TENSION FORCE = 0.10000E+06SIGMAY = 0.83333333333333IE SIGMAY^WDY1^WDY2^WDZ1^WDZ21^0.8333E+00 -.2473E-05 -.3233E-05 0.8878E-02 0.8878E-022^0.8333E+00 -.2487E-05 -.3222E-05 0.1018E-01 0.1018E-01" ELEMENT 1 "MAX. JOIST DEFLECTION ( TOTAL) =^0.13843E+02MAX. JOIST BENDING STRESS = 0.56304E+01MAX. COVER DEFLECTION = 0.13665E+02'" ELEMENT 2 "MAX. JOIST DEFLECTION ( TOTAL) =^0.14652E+02MAX. JOIST BENDING STRESS = 0.64713E+01MAX. COVER DEFLECTION = 0.13660E+02** ELEMENT 3 "'MAX. JOIST DEFLECTION ( TOTAL) =^0.88749E+01MAX. JOIST BENDING STRESS = 0.35843E+01MAX. COVER DEFLECTION = 0.11385E+02

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050474/manifest

Comment

Related Items