Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The ultimate load capacity of square shear plates with circular perforations : (parameter study) Martin, Anthony George 1985

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
UBC_1985_A7 M37_7.pdf [ 7.53MB ]
[if-you-see-this-DO-NOT-CLICK]
Metadata
JSON: 1.0062574.json
JSON-LD: 1.0062574+ld.json
RDF/XML (Pretty): 1.0062574.xml
RDF/JSON: 1.0062574+rdf.json
Turtle: 1.0062574+rdf-turtle.txt
N-Triples: 1.0062574+rdf-ntriples.txt
Original Record: 1.0062574 +original-record.json
Full Text
1.0062574.txt
Citation
1.0062574.ris

Full Text

THE ULTIMATE LOAD CAPACITY OF SQUARE SHEAR PLATES WITH CIRCULAR PERFORATIONS (PARAMETER STUDY) by ANTHONY GEORGE MARTIN B.ApSc, University Of British Columbia, 1981 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard. THE UNIVERSITY OF BRITISH COLUMBIA SEPTEMBER 1985 ©Anthony George Martin, 1985 In presenting this thesis in partial fulfillment of the requirements for the advanced degree at the UNIVERSITY OF BRITISH COLUMBIA, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her Representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Anthony G. Martin, P. Eng. Department of CIVIL ENGINEERING THE UNIVERSITY OF BRITISH COLUMBIA 2070 Wesbrook Place, Vancouver, Canada. V6T-1W5 Date: September, 1985 ABSTRACT The incremental structural analysis program NISA83 was used to investigate various parameters affecting the ultimate capacity of square plates with circular perforations subjected to uniform shear stress. Both nonlinear material properties and nonlinear geometry were taken into account in determining the ultimate in-plane capacities and buckling capacities of perforated shear plates. The parameters investigated during this study were the hole size for a concentric location, and the hole location for a constant ratio of hole diameter to plate width of 0.2. In addition various doubler plates were studied to determine the most effective shape to restore a shear plate to its original ultimate in-plane capacity. For the first two parameters, the analysis was separated into three parts. The ultimate in-plane capacity, elastic buckling capacity and the ultimate elastic-plastic buckling capacity was determined for each combination of the two parameters. These were used to identify the importance of both elastic buckling and nonlinear material contribute to the reduced ultimate plate capacities. The results from plates with a concentrically located hole of varying size showed excellent correlation with other published experimental and analytical results for both the in-plane capacity and the 3-dimensional buckling capacities. Variation of the center location of a hole of a standard size provided some significant results. Little change was found in the ultimate in-plane capacity for all hole locations. On the i i other hand, the elastic buckling capacity was raised by 50% after moving the hole from the plate tension diagonal to the compression diagonal. Finally, from the ultimate elastic-plastic' buckling capacity results it was concluded that the concentric provides lower bound capacity for all other hole locations. The in-plane analysis of the optimum doubler plate size showed wide and thin plates to be more effective than narrow and thick plates. A doubler plate with the same thickness as the plate and twice the diameter of the hole is recommended to restore the perforated plate to its original in-plane capacity. In order to aid in the tedious task of checking the input data and to provide a convenient way of displaying the result, a full graphic post-processor was developed as part of this thesis. The program NISPLOT used color graphics available at the UBC Civil Engineering lab to process the output from NISA83. It was written in FORTRAN 77, utilizing subroutines from a commercial graphics package, DI3000, to obtain device independent graphics. NISPLOT generated plots of the nodes and element mesh for each data check. When a complete analysis was carried out by NISA83, nodes, element mesh, deflected shape, and color stress fill plots were generated. 9, fe.fs-TABLE OF CONTENTS Page ABSTRACT ii TABLE OF CONTENTS iv LIST OF TABLES viLIST OF FIGURES viii ACKNOWLEDGEMENTS x i NOMENCLATURE xiCHAPTER 1. INTRODUCTION 1 1 .1 Background1.2 Purpose and Scope 3 2. THEORETICAL BACKGROUND 5 2.1 The Ultimate Behavior of Shear Plates 5 3. COMPUTER PROGRAMS , 9 3.1 NISA83 9 3.1.1 General Background 9 3.1.2 Time Step Increment3.1.2.1 Convergence and Divergence 10 3.1.2.2 Constant Load Control 11 3.1.2.3 Constant Arclength Control 12 3.1.2.4 Iteration Technique 3 3.1.2.5 Applied Increment and Iteration Algorithm 14 3.1.3 Nonlinear Material 6 3.1.4 Element Library 18 3.1.4.1 2-Dimensional Plane Stress Element 9 3.1.4.2 3-Dimensional Plate Shell Element 21 3.2 NISPLOT 27 3.2.1 General Description 23.2.2 Stress Fill Routine 9 iv 3.2.3 Visible Surface Plotting 31 3.2.4 Metafiles 32 3.2.5 Flow Charts 5 4. PLATE ANALYSIS 39 4.1 Variation of Hole Size 40 4.1.1 Plate Geometry4.1.2 Finite Element Model 40 4.1.3 Results 44 4.1.3.1 In-plane Yielding 44.1.3.2 3-Dimensional Elastic Buckling ... 47 4.1.3.3 3-Dimensional Elastic-Plastic Buckling . 51 4.2 Variation of Hole Location 56 4.2.1 Plate Geometry 54.2.2 Finite Element Model 57 4.2.3 Results 60 4.2.3.1 In-plane Yielding 64.2.3.2 3-Dimensional Elastic Buckling ... 61 4.2.3.3 3-Dimensonal Elastic-Plastic Buckling 65 4.3 Optimume Doubler Plate 70 4.3.1 Plate Geometry 1 4.3.2 Finite Element Model 74.3.3 Results 73 4.3.3.1 In-plane Yielding 74.4 Convergence with Mesh Refinement 75 5. CONCLUSIONS 79 REFERENCES 83 APPENDIX A Derivation of Consistent Shear Load Vector .... 84 for the Bicubic Isoparametric Element APPENDIX B ASCE Suggested Design Guides for Beams with Web Holes 88 v APPENDIX C Modification of NISA80 at U.B.C 90 APPENDIX D Program Listings 97 D.1 NISPLOTD. 2 MESHGEN 112 APPENDIX E Communications Programs 118 E. 1 WORDSTAR Output on the MTS Xerox 9700 118 E.2 Transfer of a VAX/VMS File to the UBC/MTS System 119 vi LIST OF TABLES Table No. Title Page No. 2.1 Definition of Plate Slenderness Failure Modes 7 3.1 Interpolation Functions for Bilinear Plane Stress Element 21 3.2 Interpolation Functions for Bicubic Element 25 4.1 Failure Mode Classification 6A.1 Four Cubic Shape Functions along the Bicubic Element Boundary s=1, — 1 <r< 1 86 vi i LIST OF FIGURES Figure No. Title Page No. 2.1 Full Plate, Ideal Elastic-Plastic Behavior 7 2.2 Perforated PLate Elastic-Plastic 3ehavior3.1 Applied Increment and Iteration Algorithm 15 3.2 Initial Stress-Strain Diagram for Mild Steel 17 3.3 Elastic-Ideal-Plastic Material Model in Uniaxial Tension Test 18 3.4 Bilinear Plane Stress Isoparametric Element 20 3.5 Bicubic, Isoparametric, Degenerated, Plate Shell Element 24 3.6 Comparison of the Classical Concept and Degeneration Finite Element Formulation 26 3.7 Subdivision of the 16 Node Isoparametric Element into 16 Sub-Regions. Each Sub-Region is Filled with a Color According to the Stress Level at the Gauss Integration Point in the Sub-Region 30 3.8 Listing of the Metafile Source File "mtr.log" that Resulted in the Compound Picture in Figure [3.9] 33 3.9 Compound Picture using a Metafile 34 3.10 Flow Chart of the Frames Plotted by NISPLOT 35 3.11 Flow Chart of the Subroutine Interaction in the Program NISPLOT 36 4.1 Perforated Plate Showing 1/4 F.E. Model 41 4.2 1/4 Plate Model using 3x3 Element Mesh 43 4.3 In-plane Stress Distribution, 1/4 Model . .' 45 4.4 Comparison of the Ultimate In-plane Shear Capacities of Concentrically Perforated Plate as Calculated by the Finite Element Method and the ASCE Design Proposal given by Equation [4.2] 46 4.5 Variation of Elastic Buckling Coefficient with Concentric Hole Size 49 viii Figure No. Title Page No. 4.6 Elastic Buckling Mode,Concentric Hole 50 4.7 Decreasing Elastic-Plastic Buckling Capacity of a Perforated Plate with Increasing Applied Load 52 4.8 Load Deflection Curve of a Simply Supported Perforated Plate with Various Concentric Hole Sizes. 52 4.9 Variation of Ultimate Elastic-Plastic Buckling Capacity of Simply Supported Perforated Plate with Concentric Hole Size 54 4.10 E-P Buckling with von-Mises Stress 55 4.11 Plate Geometry and Loading used in the Analysis of the Variation of Hole Location Parameter 56 4.12 . Finite Element Model of Half the Plate 58 4.13 Finite Element Model of the Total Plate 59 4.14 Ultimate In-plane Capacity Results for Various Hole Locations Normalized to the Concentric Hole Ultimate In-plane Capacity 61 4.15 Yield at 90% Ultimate In-plane Load 62 4.16 Full Plate Model, von-Mises Stress 3 4.17 Elastic Buckling Capacity Factors for Various Hole Locations Normatized to the Concentric Hole Capacity 64 4.18 Profile of the Tension Diagonal 66 4.19 Profile of the Compression Diagonal 67 4.20 Ultimate Elastic-Plastic Capacity for Various Hole Locations Normalized to the Concentric Hole . Ultimate Elastic-Plastic Capacity 68 4.21 E-P Buckling at 90% Ultimate Capacity 69 4.22 Geometry and Loading of Perforated Plate with Doubler Plate 72 4.23 Finite Element Mesh of Perforated Plate with Doubler Plateix Figure No. Title Page No. 4.24 Effective Capacity Restoration Factor vs Doubler Plate Area for Various Doubler Plate Diameters 74 4.25 Spread of Yield Zones for Standard and Reinforced Perforated Plates 75 4.26 Convergence of the Elastic Buckling Load with Mesh Refinement for a Concentrically Holed Plate with D/b =0.2. t/6=0.0l 78 A.1 Bicubic Isoparametric Plate Shell Element with Uniform Shear Loading along One Edge 84 B.1 Ultimate In-plane Capacity of Shear Web with a Circular Perforation as Proposed in Reference [11] . 89 x ACKNOWLEDGEMENTS The author expresses his gratitude to his advisor Dr. S. F. Stiemer for his valuable advice during the research and thesis preparation. He would also like to thank Dr. P. Osterrieder for originally suggesting the topic and for his valuable help during the analysis. Much grateful appreciation is due to the National Research Council of Canada for their financial support throughout this work. xi NOMENCLATURE A list of symbols used repetitively in this thesis is given here. Conventional mathematical indices are not included. Symbol Description A equivalent hole length Ai cross sectional area of doubler plate {Dd~D)td Ah cross sectional area of circular perforation Dt b • width of plate B(un) operator matrix D diameter of circular perforation Dd diameter of doubler plate E Young's modulus (for steel, 200,000. MPa.) F(un) internal forces in configuration n H equivalent hole height J2 second invariant of the deviatoric stresses k elastic buckling coefficient K{un) incremental stiffness matrix using displacements un Kr(un) elastic stiffness matrix Kg{?'n) geometric stiffness matrix N(r,s) element shape functions P load vector AP increment in the load vector for one time step Q{un) out-of-balance load vector P - F(un) r,- radius from the center of the hole to node t (polar coordinates) R distance from center of plate to center of perforation S factored ratio of hole diameter to plate width 0.9D/6 xii Symbol Description S(u") stress matrix Sr.S..,S, deviatoric stresses t thickness of plate t, thickness of doubler plate elastic buckling stress of full plate elastic buckling stress of perforated plate exact elastic buckling stress of perforated plate ultimate in-plane capacity of perforated plate with doubler plate r„ ultimate elastic-plastic buckling stress of full plate r„ ultimate elastic-plastic buckling stress of perforated plate Ty yield shear stress of steel or in-plane ultimate stress of full plate Ty in-plane ultimate stress of perforated plate Txy,Tyz,Tzx shear stresses u exact equilibrium displacement vector u" approximate equilibrium displacement vector Aun increment in approximate equilibrium displacement vector a factor A load factor v Poisson's ratio (0.3) trm mean stress cr„ yield stress (for steel, 300. MPa.) °~zzi°'yy}0'zz normal stresses 1* doubler plate capacity restoration factor xiii 1 INTRODUCTION 1.1 Background Flat rectangular plates are very common structural component in modern buildings and offshore structures and are often subjected to a pure uniform shear stress. These plates are normally referred to as shear plates or shear webs. Cost effectiveness requires optimal use of space and weight of a structure. This often results in one or more holes being cut out of components such as shear webs of girder webs. These holes then allow utility pipes or ducting to pass through the web and reduce the overall structural weight. These shear plates with plain or reinforced circular holes are called perforated shear plates. Early elastic in-plane stress analysis of square perforated shear plates was done by WANG [1], who applied a finite Airy's stress function over the domain of the plate. His results show that the stress concentrations increased very rapidly with the hole size ratio D/b (hole diameter to plate width). He carried out an investigation into the stress distribution produced with the addition of reinforcement around the hole. He found that the circumferential stress is significantly decreased about the hole, but the doubler plate increases the radial shear stress from that of the unreinforced perforated plate. The stability of very slender perforated shear plates was first investigated by ROCKEY, ANDERSON and CHEUNG [2], They used the finite element method to establish the relationship between the elastic buckling load of the plates and the hole size ratio, Djb, for both clamped and simply supported boundary conditions. 1 They also investigated the relationship between the hole reinforcement shape and the elastic shear buckling coefficient, jfc. In addition to this analytical work they carried out an experimental study on perforated shear plates. Since their analytical model was only able to handle an elastic material model, they chose the plate slenderness ratio small enough so that neither the test plate nor the analytical plate would be subject to material yielding. The combined effects of the material plasticity and buckling stability were first dealt with by UENOYA and REDWOOD [3]. Using a 2-dimensional in-plane finite element program they calculated the elastic-plastic stress distribution due to shear loading. These stresses were then inserted into a Rayleigh-Ritz energy expression. By minimizing the potential energy of this system they were able to determine the elastic-plastic buckling capacity of the perforated shear plate. This method was applied to a one quarter model of a square shear plate for both simply supported and clamped boundaries. The method was later extended to rectangular beam webs with combined shear and moment loading for both circular and rectangular perforations. URENOYA and REDWOOD correlated their analytical results with an experimental study carried out concurrently, and with past experimental work. Their two part solution provided excellent results for perforated plates with small perforations. However, as the hole size ratio, D/b> became greater than 0.7, the restricted number of terms used in the Fourier series approximation of the deflected shape produced unexpected results. 2 For the clamped plate they found that by increasing the hole size the ultimate elastic buckling capacity is increased rather than decreased, as had been expected. The result from UENOYA and REDWOOD show that any circular perforation in a shear plate has a significant effect on the stability and ultimate load capacity of the plate. They also note that with increasing hole size there is a considerable increase in the amount of plasticity developed in the plate before the ultimate capacity is reached. 1.2 Purpose and Scope In the research of this thesis the Nonlinear Incremental Structural Analysis program, NISA83 [7], is used to carry out a parameter study on square shear plates with circular perforations. The analysis considers three separate loading capacities; in-plane yielding, elastic buckling and elastic-plastic buckling, for the two major parameters, hole size and hole location. The analysis of the third parameter, optimum doubler plate shape, only considers the ultimate in-plane capacities. Finally, an estimation of the accuracy of the elastic buckling analysis is obtained by refining the element model in order to establish the convergence rate and exact elastic buckling capacity. Included as part of this thesis is the development and implementation of a color graphics post processor called NISPLOT. The program was originally developed to check the input data. However, it was later extended to provide post processing for NISA83 output. NISPLOT is used throughout the research to provide simple plots of the model nodes and element meshes for data 3 checks. In addition, it is used to illustrate and summarize some of information provided by the finite element analysis. The program currently runs on a VAX 11/730 under the EUNICE operating system. It is written in FORTRAN 77, making use of the graphics package DI3000, by Precision Visuals [8], 2 THEORETICAL BACKGROUND 2.1 THE ULTIMATE BEHAVIOR OF SHEAR PLATES The theoretical behavior of square plates subjected to uniform shear loading along the boundaries can be separated into two types of behavior. If the square shear plate has no initial imperfections its ultimate capacity is governed by the lesser of the ultimate material shear stress or the elastic shear buckling stress of the plate. According to the von-Misps vield criteria the ultimate material shear stress is Ty = cry/\/3. ROCKEY [2] defines the elastic shear buckling stress by equation [2.1]. k = 9.34 for simply supported square plate By equating the ultimate material shear stress to the elastic buckling shear stress, a simple expression, equation [2.2], for defining the balance or transition point, in terms of the plate slenderness (t/b), is obtained. Theoretically, at this transition point the material will yield at exactly the same load level that the plate buckles. Equation [2.2] is also used to define the slenderness ratio above which a plate is considered stocky, or below which it is considered slender. The ultimate capacity of a stocky plate is governed by the ultimate material shear stress. A slender plate has a lower elastic buckling stress than material yield stress, so its ultimate capacity is governed by equation [2.1]. Illustration of the effect of slenderness on failure mode of ideal stocky and slender plates are shown in figure [2.1]. 5 CT y/3 12(1-0.32) E b cry _ 9.35TT2JE 0.2614 (2.2) In reality however, plates are never perfectly flat and the loading is never perfectly uniform. With the introduction of initial imperfections in the plate surface or eccentricities in the loading. The transition between yielding and buckling is no longer a single point. Now the transition from one failure mode to the other takes place over a range of slenderness values, depending on the degree of deviation from the perfect condition. However, the purpose here is not to examine the effects of imperfections on the unperforated plate. The analysis will assume an initially perfectly flat plate with uniform shear loading. The effects of imperfect conditions will be the subject of future studies. The effect of introducing of a circular hole in a plate is similar to that of imperfections. There is no longer a distinct point to distinguish between elastic buckling and material yielding. For intermediate slender plates, significant material yielding occurs around the hole. This causes a redistribution of the stresses and change in the load path. The net result of this is the reduction of the plate buckling capacity. In this transition region, referred to as zone II, the ultimate capacity of the plate is controlled by the ultimate elastic-plastic buckling stress. Unlike the elastic buckling stress in equation [2.1], no analytical solution existed for the zone II ultimate capacities. Therefore, solutions for the ultimate elastic-plastic 6 buckling capacities in zone II must be determined numerically. Fig. 2.1: Full Plate, Ideal Fig. 2.2: Perforated Plate, Elastic-Plastic Behavior Elastic-Plastic Behavior Table 2.1: Definition of Plate Slenderness Failure Modes ZONE FAILURE MODE I II III material yielding elastic-plastic buckling elastic buckling The objective of this research is to identify the ultimate plate capacities in each zone. By restricting the displacements or varying the material properties only one finite element model is required for each parameter change. The material capacity of zone I is determined by restricting the model to in-plane displacements only. The elastic buckling capacity of zone III is 7 evaluated by applying a small load increment to the model and then performing a bifurcation analysis on the resulting stiffness matrix. Finally the ultimate capacity in zone II is established by a full three dimensional elastic-plastic ultimate analysis. However, in order to observe the out-of-plane behavior a small lateral load is applied at or near the center of the plate. 8 3 COMPUTER PROGRAMS 3.1 NISA83 Nonlinear finite element programs have been available for some time and are easily accessible. The program, NISA83 used in this study, contains several special features which make the program well suited for the elastic-plastic instablity problem associated with perforated shear plates. 3.1.1 General Background NISA83 (Nonlinear Incremental Structural Analysis) is installed on the, University of British Columbia Civil Engineering Department's, VAX 11/730. The program is a development of the Instituts fur Baustatik, Universitat Stuttgart, under the direction of Dr. Hafner, Dr. Ramm and Dr. Sattele over the previous ten years. NISA83 is written in standard FORTRAN 77 and designed to be supported by a variety of main frame computers. A number modifications were necessary to get the program operating on the U.B.C. VAX. A comprehensive list of all changes made to NISA83 subroutines are included in Appendix C. 3.1.2 Time Step Increment In a non-linear analysis it is not always possible to determine the final equilibrium position of the structure by applying the ultimate load in one step. There maybe more than one final equilibrium configuration, or several load paths leading to the same equilibrium configuration. Like other non-linear programs, NISA83, follows the load deflection path by breaking the problem in to a series smaller increments. These increments are referred to as time steps 9 Each time step is considered as a separate problem. At the beginning of the time step the equilibrium displacements, forces and stresses are all known. Depending on what increment control is used the end of the time step is defined by either increasing the load or by incrementing the load-displacement vector by some constant. Within each time step several iterations of the solution may be required before convergence to the defined new equilibrium configuration is astablished. 3.1.2.1 Convergence And Divergence As NISA83 completes each iteration within a time step, checks are made for convergence or divergence of the solution. The NISA83 criteria used to define if the new equilibrium configuration has converged or, the solution has diverged are as follows: Convergence - if the change in the displacement vector between iterations is less than the user specified value (default RTOL = 0.001) e = iL_!i < RTOL IMI Divergence - if after 8 equilibrium iterations the out of balance load vector is still larger than incremental applied load vector IIAPII < - if after after a specified number of iterations the convergence requirement have not been satisfied (default 15 iterations) By using the restart option the analysis can be continued 10 from any last known equilibrium configuration. If at the end of an iteration the solution satisfies the convergence criteria NISA83 updates the restart file with the new equilibrium displacements and stresses. However, if a solution fails to converge by one of the criteria above, the program stops. Thus, whenever the program stops or complets its task the restart file contains the last known equilibrium configuration. By redefining any of the time step parameters and restarting NISA83 the analysis is continued from this last configuration. 3.1.2.2 Constant Load Control The constant load increment is the most common time step control used in finite element programs. The method allows the user to specify the load level used to define the end condition of each time step. The load level is held constant throughout the time step, while the solution iterates to the new configuration until convergence or divergence is established. The constant load increment method works well for problems where the nonlinearities remain small within each time step. Normally the load is specified in constant increments along the load path. As the applied load approaches the ultimate load, the increments are reduced in magnitude. If a load level is specified that is higher than the ultimate capacity of the problem the program is unable to converge to an equilibrium solution. For this case the ultimate load is defined as the average of the last known equilibrium load step and the load of the divergent step. Problems are encountered with this method when the load level approaches a bifurcation point, or the behavior becomes 11 highly nonlinear. In the these cases small changes in the load increment produce large deformations. These deformations in turn caused a further reduction in the system stiffness resulting in still higher deformations. Thus, convergence to the equilibrium solution close to the ultimate capacity becomes very difficult using this method. 3.1.2.3 Constant Arclength Control One of the special features of NISA83 is the Riks-Wempner Constant Arclength time step control. This method lets the user specify the initial arclength to be used as a reference. In each time step the applied load vector is assumed a variable along with the displacement field. The load level is scaled within each time step so the normal of the incremental displacement vector plus incremental load vector is equal to the defined reference arclength. This method is extremely powerful in its ability to follow the load deflection path. As explained before, when using the constant load time step, once the specified load level is higher than the ultimate load it is impossible for the program to converge to an equilibrium solution. However, with the constant arclength controlled method, and the load level as an unknown, the applied load increment is reduced until equilibrium is obtained. This allows the program to obtain an equilibrium solution even if the stiffness matrix has gone negative. This mean that in the buckling shear plate problem, the program is able to follow the load path to the plate ultimate buckling 12 capacity and then on into the post-buckling region. 3.1.2.4 Iteration Technique In the two time step control methods above, the required end conditions are specified. Both methods assume that with each iteration the solution will improve and eventually the convergence requirement will be satisfied. In NISA83 the user can select one or a combination of the Modified Newton-Raphson and the Standard Newton-Raphson iteration methods. The first method uses less CPU time but the second method converges more rapidly if the tangent stiffness matrix is changing rapidly with each iteration. The difference between the two iteration techniques is illustrated by steping thougth a time step iteration. Associated with each time step are three configurations. (1) is the last known equilibrium configuration and all the information on stresses, strains, and displacements is known. (2) is the next unknown equilibrium comfiguration on the load path. Only the end condition (load level or arclength) is known in this position. (n) is some position in between (1) and (2). The procedure to get from configuration (1) to (2) is as follows: Step 1 calculate the unbalanced forces for the time step. Q(ul) = P2 - (3.1) Step 2 formulate the tangent stiffness matrix in configuration (1) K(ul) = Ke(ul) + Kg(ul) (3.2) 13 Step 3 solve for the incremental displacements Au=[^(u1)]"1P2 (3.3) Step 4 find the total displacement vector for the new configuration (n) un = u1 + Au (3-4) Step 5 calculate the internal forces in configuration (n) F(un) = J [B{un)]T S{un) dv (3.5) Vol Step 6 check for convergence or divergence, section [3.1.2.1] Step 7 calculate the remaining unbalanced forces for this time step Q(un) = P2 - F{un) (3.6) Step 8 form the new tangent stiffnes matrix K(un) (Standard Newton-Raphson) or use the old stiffness matrix K^u1) (Modified Newton-Raphson) Step 9 solve for the next displacement increment and continue with step 4 to 9 untill convergence. 3.1.2.5 Applied Increment and Iteration Algorithm Figure [3.1] show how the various increment controls, iteration techniques and the restart option are combined to provide an efficient algorithm. This algorithm is applied to all the in-plane and elastic-plastic buckling analyses. The approach provides an efficient method of raising the load level to the beginning of the nonlinear portion of the load path. It also provides an accurate procedure to follow the load path in this region beyond into the post-buckling range. The constant load step control is combined with the modified 14 MODIFIED NEWTON-RAPHSON WITH CONSTANT LOAD 1/ Deflection RESTART Deflection Fig. 3.1: Applied Increment and Iteration Algorithm 15 Newton-Raphson iteration technique for the first part of the analysis. The load level is specified to approximately 80% of the estimated ultimate capacity. Since the material yielding is confined to small areas around the hole most problems will behave almost elastically up to this load. Then the program parameters are changed so that the constant load step control is combined with the Standard Newton-Raphson iteration procedure. Then NISA83 is restarted and the analysis continues from the last known equilibrium configuration. The method was found to be adequate for up to 95% of the ultimate load, using increments in the order of 2% of the current specified load level. After the load step control fails to converge, the constant arclength time step control is invoked. This time step control combined with the Standard Newton-Raphson iteration is then used to follow the load path for the rest of the analysis. This method is able to follow the true load-deflection path until the the plate buckles and beyond into the post-buckling region. 3.1.3 Nonlinear Material In order to identify the point at which elastic deformations stop and plastic deformations start, for various stress states, a yield criteria is used. The von-Mises yield criteria, used by NISA83 to identify plastic strains, is well accepted as a reasonable model of the elastic plastic behavior of steels. Mathematically the criteria can be expressed by equation [3.7] 16 *9/y/3 = J21/2 = y/l/2 (5| + SJ + S?) + r%, + rj, + r* (3.7) = C__ — (7, x xx : m m 5* = - °"i m NISA83 also allows the user to select the usual material parameters E, <ry , and v along with strain hardening modulus Eh • For mild steel if the strains are between 0.002 and 0.02, then the stress-strain curve is almost horizontal. Most strain hardening effects occurs at strain levels above 0.02, as shown in figure [3.2]. Since all of the plates considered were expected to have strains below this level no strain hardening is considered. Thus, the material model commonly known as "elastic ideal plastic" is used throughout the analysis. Under uniaxial tension of this material model give the stress-strain diagram shown in figure [3.3]. 400 -o 0.000 0.005 0.010 0.015 0.020 0.025 Strain mm/mm Initial Stress-Strain Diagram for Mild Steel Fig. 3.2: 17 o co CO 400-300 -200 -100 0.000 0.005 0.010 0.015 Strain mm/mm 0.020 0.025 Fig.3.3: Stress-Strain Diagram for Elastic Ideal Plastic Material Model in Uniaxial Tension The values of the parameters used to define the material properties are; cry= 300. MPa. E= 200,000. MPa. Eh= 0.0 MPa. v= 0.30 3.1.4 Element Library NISA83 contains an extensive element library. Elements range from truss, beam, and curved beam elements for longitudinal elements and, 4 and 8 node isoparametrics with plane stress or plane strain for 2-dimensional analysis to, a family of degenerate plate shell elements with 4,8,9 or 16 nodes for problems involving membrane forces and plate bending. A short description of the two elements used in this research is presented here. 18 3.1.4.1 2-Dimensional Plane Stress Element For the two dimensional analyses of doubler plate capacities a simple plane stress element is used. This 4 node isoparametric element is thoroughly documented as being well-conditioned and giving reliable results for plane stress problem [10], This element is used to discretize the 2-dimensional domain of the perforated plate with a doubler plate, and determine the ultimate in-plane plate yield capacity. The element is derived from the interpolation functions given in table [3.1]. As for all isoparametric elements the functions are used for mapping the element from the global x-y coordinates to the local r-s coordinates, and also for interpolation of the element displacement field, equation [3.8,3.9], In order to maintain the energy bound of the finite element formulation, 2x2 Gauss numerical integration is used with this element. This order of integration exactly evaluates the element stiffness integral provided the coordinate transformation matrix is a constant. Some higher order error is introduced if the coordinate transformation matrix is not constant. However, if the element mesh is refined the trasformation matrix approches a constant and the error becomes negligible. Thus, this order of numerical integration is sufficient to ensure convergence as the (3.8) (3.9) 19 mesh is refined and preseve the energy bound. A full derivation of the complete element stiffness matrix is given in reference [10]. Table 3.1: Interpolation Functions for Bilinear Plane Stress Element Element Node No. Coordi r nate s Interpc 1 +r )lation E 1 -r 'unction 1+s * 4.0 1-s 1 1 1 1 - 1 -2 -1 1 - 1 1 -3 -1 -1 - 1 - 1 4 1 -1 1 - - 1 3.1.4.2 3-Dimensional Plate Shell Element For all 3-dimensional and some 2-dimensional analysis, the higher order degenerated plate shell element is used. BATHE [12] found that the 16 node, bicubic isoparametric plate shell element provides reliable results for 3-dimensional analyses. However, inorder to eliminate variations in results due to modelling changes, this same element and mesh is also used for some in-plane ultimate capacity calculations. By applying the appropriate boundary conditions to the nodes, the 3-dimensional model is restricted to in-plane displacements and can be successfully used for the calculation of the in-plane yield capacities. The interpolation functions, called "Lagrangian cubic polynomials", used for the 16 node element formulation are given in table [3.2], The term "degenerated element" refers to the way element is formulated. The classical formulation of a finite element is based on plate or shell theory. This theory is developed from the general 3-dimensional field equation by making 21 some restricting assumptions (ie. small displacements and ignoring higher order terms) to reduce it to a 2-dimensional equation on the shell surface. A structural system is then discretized using the elements in this 2-dimensional domain. The degenerate formulation is more general than the classical approach. No behavior or displacement assumptions are used to reduce the 3-dimensional field equation to a 2-dimensional domain. Instead the formulation is based on the full 3-dimensional field equations. Therefore, any 3-dimensional structural system is discretized directly by a finite element model. A schematic comparison of the two formulations is shown in figure [3.6] Geometrically nonlinear problems are solved using large deflection theory which consider the second order strain displacement terms that are neglected in linear analyses. NISA83 allowed the user to apply either the Updated Lagrangian or Total Lagrangian large deflection formulation. Here the Updated Lagrangian formulation is used for all buckling analysis and a brief description follows. The Updated Lagrangian Formulation is based on the element local coordinate system tracking with the element incremental deformations. At the beginning of each time step the element coordinate system, called the "corotational system", is updated to the last known equilibrium configuration. All information about stresses, strains and the equilibrium displacements are known in this configuration and it is used as a reference throughout the increment. At the end of the time step, when a new 22 equilibrium configuration is established, the coordinate system, strains and stresses are updated or transformed to this new configuration, ready to start the next time increment. With this large deflection theory the element corotational system shares the same rigid body modes as the element, thus any stresses or strains calculated are with reference to the last known deformed equilibrium configuration of the element. A 4x4 Gauss numerical integration in the r-s plane and a 5 point Simpson integration through the thickness is used with this element. Again this order of integration exactly evaluates the element stiffness integral if the coordinate transformation matrix is a constant. The Simpson numerical integration in the thickness direction provides reliable results if the width to thickness ratio of the plate is greater than 10. 23 10 In local space (projection on to the plane t = 0) Fig. 3.5: Bicubic, Isoparametric, Degenerated, Plate Shell Element 24 Table 3.2: Interpotation Functions for Bicubic Element. Elem Node No. coc r rd. s Fact I nte 1+r rpola 3r+1 t ion 3r-1 Funct 1 -r ion * 1 +s 256 3s+1 3s-1 1 -s 1 1 1 1 1 1 1 - 1 1 1 -2 -1 1 - 1 1 1 1 1 1 -3 -1 1 - 1 1 1 - 1 1 1 4 1 1 1 1 1 - . - 1 1 1 5 1/3 9 1 1 - 1 1 1 1 -6 -1/3 -9 1 - 1 1 1 1 1 -7 -1 1/3 9 - . 1 . 1 1 1 1 - < 8 -1 "1/3 -9 - 1 1 1 1 - 1 9 "1/3 -1 -9 - 1 1 - 1 1 < 1 0 1/3 -1 9 -1 - 1 - 1 1 1 1 1 -1/3 -9 1 1 - - 1 12 1 1/3 9 < 1 1 - 1 --1 3 1/3 1/3 81 1 - 1 -1 -1 4 -1/3 1/3 -81 < - 1 1 1 -1 5 -1/3 "1/3 81 - 1 1 - 1 -16 1/3 -1/3 -81 • 1 - 1 < - 1 25 CLASSICAL CONCEPT DEGENERATION Assumptions Analytical reduction _3_dim—— 2 dim over thickness Discretiza tion 2 dim —point numerical j* over surface Dis place mentL*. solution Assumptions Discretiza tion 3 dim point numerical ^ over volume 3.6: Comparison of the Classical Concept and Degeneration Finite Element Formulation 26 3.2 NISPLOT 3.2.1 General Description As commonly known, output generated by finite element program can be very lengthy. The program NISPLOT was developed to use the color graphics capabilities of the Department of Civil Engineering, U.B.C, to process the output from NISA83 into a more concise form. NISPLOT was used for checking the input data before running an analysis. After an analysis it was used as a color graphic post-processor, producing element mesh, 3-dimensional deflected shapes, or element stresses (displayed in up to seven colors). NISPLOT is written in a very general and modular form so that any changes or subsequent modifications that are required can be carried out quickly and efficiently. The program is written in standard FORTRAN 77 programing language and makes calls to standard graphics subroutines in the commercial graphics software package, DI3000 by Precision Visuals [8], NISPLOT runs on the University of British Columbia VAX 11/730 under the EUNICE 4.1 environment, which is a UNIX emulator. It is nessesary to the EUNICE operating system as opposed to the native VAX-VMS system because of the storage restrictions on the VAX. The DI3000 subroutine library requires a large amount of disk storage and only one object library can be stored on the VAX 11/730 system disk. Since the Civil Engineering Department is currently in the process of converting to the EUNICE operating system, only the EUNICE library is stored on the disk. On the other hand, NISA83 only runs under the VAX-VMS system. However all input and output files are VMS to EUNICE 27 compat ible. Like any other program that makes calls to the graphics subroutines of DI3000, the object module is linked to the correct device driver to create the device dependent executable module. At the Civil Engineering Graphics lab the program can be linked to a Seiko D-SCAN GR-1104, Silicon Graphics Inc. IRIS-1200, and a Watanabye Instrument Corp. WX4636R plotter. The image could also linked to any other the device driver supported by DI3000. NISPLOT is sequentially structured and prompts the user for the input data files. The default NISA83 output geometry, displacement, and stress files are "plot.geo", "plot.dis", and "plot.str" respectively. The user can supply one, two or all three of the input files when requested, but the geometry file must always be one of the files. The program checks to see if the files specified exist and prompt the user if it does not. Unfortunately, the user must-be sure that the correct file type is given when prompted, as the program does not check the validity of a file type. When given all three files as input, NISPLOT will plot the following frames: 1) the mesh nodes and node numbers 2) the element mesh and element numbers with each element group in a new color. 3) a solid color stress fill of the 16 node isoparametric plate shell element 4) 3-dimensional plot of the deflected shape with 2 colors, one for the near and one for the far side, of the 16 node isoparametric plate shell element. 28 5) 3-dimensional deflected shape with a seven color solid fill of the 16 node isoparametric plate shell element The program routines for plotting the element nodes, element mesh and the title are universal routines and work for any of the elements contained in the NISA83 library. The other routines involving plotting the element stresses, deflected shape with solid fill, or deflected shape with color stress fill, are element dependent, and are only supported for the 16 node isoparametric plate shell element with 4x4 Gauss integration. 3.2.2 Stress Fill Routine This routine fills each sub-region of all the 16 node elements with a color that is consistent with the von-Mises stress level calculated at the 4x4 Gauss integration points. First all the element stress are searched to determine the maximum and minimum values of the von-Mises stresses. Then this range is divided into six levels and a color assigned to each level. Next the routine proceeds through all the elements as follows. Each 16 node isoparametric element is subdivided into 16 sub-regions. Each region contains one Gauss integration point. The boundaries of these 16 sub-regions are defined by 25 nodes as shown in figure [3.9]. Calculation of these points is done using the cubic interpolation functions given in table [3.2], By examining the von-Mises stress at the gauss integration point contained within a sub-region the routine selects the correct color and fills the sub-region area. 29 <5> 6 1—Q 6 19 7 O 7 20 8 O + o—t-3 I I + 1 4 1 3 + 17 16 -I-21 1 5 -4.-4--25 10 6" 1 6 I + + 11 o 9 Oi2 + 22 15 --T- 4 + <V1 ^ 23 14 + 1 12 1 3 4> 10 O Global Nodes # Subdivided Element Nodes (25 nodes, 16 regions) -f- Gauss Integration Points Fig 3.7: Subdivision of the 16 Node Isoparametric Element into 16 Sub-Regions. Each Sub-Region is Filled with a Color According to the Stress Level at the Gauss Integration Point in the Sub-Region. 30 3.2.3 Visible Surface Plotting The subroutine VISBLE is a simple routine that attempts to improve the quality of any three-dimensional plot by distinguishing between surfaces that are visible in the viewing plane and those that are not. By examining the normal of any given surface in the viewing plane the subroutine is able to detect if the top or bottom surface is visible to the viewer. The hidden surface subroutine VISBLE is called from both the stress fill routine and the solid fill routine just before the element sub-region is filled. The global x,y,z coordinates of three or more points that lie on a sub-region surface are passed to this subroutine. These points are then converted into u,v coordinates of the local viewing system. The subroutine then calculates the cross product of two vectors formed by point(1)-point(2) and point(3)-point(2). If the result is positive the top surface is visible on the viewing plane, and if it is negative the bottom surface is visible on the viewing plane. By calling the solid fill routine twice, once filling all the visible bottom surfaces with dark blue and then filling all the visible top surfaces with light blue, it is possible to create plots in which the two surfaces are easily distinguished. Similarly, by calling the stress fill routine twice and using a different fill pattern for top and bottom surfaces a clearer illustration results. The routine works well for the simple convex deflected shapes of the plate problem. For complex surfaces, however, a more sophisticated algorithm is required, but this is beyond the scope of this thesis. JENSSEN [6] provides a good reference for future work in this area. 31 3.2.4 Metafiles A metafile is a device independent plot file which is stored on the host and may be kept as a permanent file for future reference. The idea of the Metafiles is that they can be recalled and manipulated by anyone using the DI3000 metafile translator without having to run NISPLOT or having access to the NISA83 output files. When the user requests that a Metafile be created, NISPLOT creates and opens a new file called "NISAPLOT.MFL". All the subsequent frames drawn by NISPLOT on the current device are then stored in this file, as well as appearing on the viewing device. These plot frames, or pictures, are stored in the standard format of the-DI3000 graphics package. When NISPLOT finishes running it closes the file. If a metafile is requested, the plotting sequence is modified. NISPLOT plots the title only in the first frame. The reason for this will become apparent after further discussion of the manipulations that can be done with the translator. A detailed description of the translator can be found in the DI30.00 Users' Manual [8] under the chapter "Metafiles". Only a brief explanation of possible manipulations is outlined.. To recall pictures from the metafile, NISAPLOT.MFL, to a new device, the DI3000 metafile translator linked to the device is invoked (ie. mtrans.seiko). The metafile translator allows the user to define a viewing port and window and then request any of the pictures from the Metafile be drawn in these viewing attributes. By dividing the screen of the viewing device into two or four viewing ports and requesting that different picture be 32 drawn in each port, a compound picture is generated. However, since the title was only drawn on the first picture generated, it will not appear in any of the view ports. To place the title on the subsequent compound picture the user can define a new view port for the full device screen and then select picture one to be drawn. This causes the title to appear under the total compound picture instead of in each sub-picture. One last feature about the metafile translator is that the commands to define the metafile, view ports, windows and even the order of pictures drawn can be set up in a log file. When the Metafile translator is activated, the user assigns this file as a command source. The translator will then execute these commands in sequence. This source file can be set up to draw pictures from several Metafiles in sequence to display a series of results. SET MF 1 NISAPLOT.MFL SET W 1 (-1.0 1.0 -1.0 1.0) SET W 2 (-1.0 1.0 -0.95 0.95) SET V 1 (-1.0 0.0 0.05 1.0) SET V 2 ( 0.0 1.0 0.05 1.0) SET V 3 (-1.0 0.0 -0.9 0.05) SET V 4 ( 0.0 1.0 -0.9 0.05) SET V 5 (-1.0 1.0 -1.0 1.0) SET DEFAULT W 2 SET DEFAULT BOX ON DRAW P 2 V 1 NOEJECT DRAW P 3 V 2 NOEJECT DRAW P 4 V 3 NOEJECT DRAW P 5 V 4 NOEJECT SET BOX OFF DRAW P 1 V 5 W 1 Fig. 3.8: Listing of the Metafile Source File "mtr.log" that Resulted in the Compound Picture in Figure [3.9] 33 Fig. 3.9: Compond Picture using a METAFILE 3.2.5 Flow Charts Figure [3.10] shows a flow chart of NISPLOT with respect to the different frames plotted. Figure [3.11] shows a more detailed flow chart of NISPLOT showing how subroutines interact. ( START) DEFINE INPUT FILES PLOT AND NUMBER NODES PLOT AND NUMBER ELEMENTS IF (3dstr) NO PLOT ELEMENTS WITH SOLID COLOR STRESS FILL PLOT UNDEFLECTED SHAPE IN 3-D THEN DEFLECTED SHAPE WITH ONE COLOR SOLID FILL PLOT UNDEFLECTED SHAPE IN 3-D THEN DEFLECTED SHAPE WITH SOLID COLOR STRESS FILL Fig. 3.10 Flow Chart of the Frames Plotted By NISPLOT 35 ( START ) /SETUP/ T /READNO/ MAXMIN VIEW PLTNOD - NODNUM PLTHED /READEL/ • PLTELE r PLTHED DRAY ELENUM DATIN ELESTR SHAPE 1 * VISBLE PLTELE DRAY LEGEND J^d5plt),  /READNO/ VIEW r ADDDIS PLTELE ( PLTHED ADDDIS FILLEL > PLTELE DRAY ELENUM VISBLE ELESTR ( PLTELE -LEGEND (eof.ne. NO <END) Fig. 3.11: Flow Chart of the Subroutine Interaction in the Program NISPLOT 36 Descriptions of subroutines in figure [3.10] Subroutine Description SETUP asks the operator for the input file names containing the geometry, displacements and stresses, and opens these files READNO reads the x,y,z coordinates of the nodes MAXMIN determines the maximum and minimum global dimensions of the node coordinates VIEW sets the frame viewing attributes ie. 2-D or 3-D view, window, view port etc. PLTNOD marks the node point with "+" NODNUM numbers the node point marks PLTHED plots the heading at the bottom of the page READEL reads in the element node numbers PLTELE sets the line style and color attributes for each element and collects the node numbers for each line in one array, ready for DRAY DRAY draws the line given by PLTELE ELENUM numbers the elements at their mid point READST reads in the element stresses ELESTR separates the 16 node isoparametric element into 16 sub-regions with 25 nodes, then fills each of the sub-regions with the appropriate color according to the stress level at the Gauss point in the sub-region DATAIN initializes variables used in ELESTR SHAPE converts the 16 node isoparametric element to a 25 node element using the element shape functions for interpolation VISBLE checks to see if the sub-region to be plotted is visible in the viewing plane 37 Subroutine Description LEGEND writes the stress legend in the upper right-hand corner of the frame ADDDIS scales and adds the nodeal displacements to the original nodal coordinates FILLEL separates the 16 node isoparametric element into 9 sub-region and fills each sub-region with one of two colors depending on the response from VISBLE 38 4 PLATE ANALYSIS Numerical analysis was performed on a series of standard square plates with circular perforations to establish the effects of three parameters on the ultimate capacities of shear plates. For the first two parameters (hole size and hole location) the ultimate in-plane, elastic buckling, and ultimate elastic-plastic buckling capacities were determined. For the third parameter (the shape of a doubler plate) only the in-plane ultimate capacities was determined. The plate thickness was selected so that the thickness to width ratio (t/b ) was close to the ideal balanced slenderness ratio, as described in section [2.1]. The outside dimensions for the standard plate, used throughout the analysis, were 1000 mm x 1000 mm x 10 mm. This provided a plate slenderness of 1/100, which was very close to the balance ratio for the full plate of 1/98.7, given by equation [2.2]. This slenderness ensured that the resulting analyses dealt with the combined material and buckling failure modes. The investigation of the first parameter, hole size, was done on the standard square plate described above with a concentric hole. The diameter of the hole was varied from 0.156 to 0.90b. For the hole location parameter, the center of a 0.26 diameter hole was moved about the surface. Many different geometries and models were used during the study of the hole size and hole location. In order to determine the significance of nonlinear material or geometry on the capacity of the perforated plate, three ultimate capacities were calculated for each variation of the two 39 major parameters. For each model configuration, the ultimate in-plane yield capacity, elastic buckling capacity and the ultimate elastic-plastic buckling capacity were calculated. The in-plane yield capacity indicate the influence of material yielding around the hole on the plate capacity. The ultimate elastic buckling capacity was used to establish the change in the elastic out-of-plane stiffness due to the presence of the perforation. Finally, the two factors were considered together with a calculation of the ultimate elastic-plastic buckling capacity. 4.1 Variation of Hole Size The first parameter investigated for its effect on the ultimate plate capacity was the size of a centrally located perforation. Results for the ultimate in-plane, elastic buckling, and elastic-plastic buckling capacities were obtained for hole diameters from 0.156 to 0.96. 4.1.1 Plate Geometry The standard square plate, 6=1000 mm and r=l0 mm, was analyzed with concentric hole diameters. The ratio of hole diameter to plate thickness, D/b, varied from 0.15, 0.2, 0.3, 0.4,...to 0.9. All material parameters were held constant to the specifications in section [3.1.3], 4.1.2 Finite Element Model By observing symmetry of the plate geometry and loading, and by applying appropriate boundary conditions, only one quarter of the plate needed to be analysed. The two diagonals of the square plate (see figure [4.1]) are two axes of symmetry for the plate loading and geometry. For the symmetric buckling modes there will 40 Fig. 4.1: Perforated Plate Showing 1/4 F.E. Model be no rotation about these axes, and these two axes remain perpendicular to each other for all deformations due to this loading. As shown in figure [4.2] the 16 node plate shell element in a 3x3 grid was used to model a quarter of the plate. This same 9 element mesh was used for all the concentric hole analyses. The boundary conditions for one model were varied to suit the ultimate load calculation desired. For 3-dimensional buckl_ing, displacement boundary conditions were imposed along the plate outer edge and the two axes of symmetry. Then, by restricting the displacement field to in-plane movement only the same element mesh was used for the analysis of the perforated plate ultimate in-plane capacity. The quarter plate model made of a 90 deg. section of plate was separated into nine slices by 10 equally spaced radial lines. Nodes were then placed along these lines using a proportional spacing given by equation [4.1]. This provided a dense concentration of nodes around the perforation where the stress gradients were the highest and material yielding was most severe. (4.1) 42 4.2: 1/4 Plate Model using 3x3 Element Mesh The consistent load vector for a uniform shear stress was applied along the plate boundaries. Using such a distorted element mesh required the use of the exact consistent load vector for each element to correctly model the uniform stress traction along the plate boundary. A general formulation for the element consistent load vector, defined in terms of element global coordinates, was developed and is documented in Appendix A. The resulting equations [A.8] were used with the element global coordinates to calculate each element consistent load and then this was added to the total model consistent load vector. 4.1.3. Results The ultimate in-plane capacity, ultimate elastic buckling capacity and the ultimate elastic-plastic buckling capacity was calculated for each variation of the hole size. 4.1.3.1 In-plane Yielding The hole size had a direct effect on the ultimate in-plane capacity of the perforated plate. With increasing hole size there was found to be a corresponding decrease in ultimate capacity. The results of the finite element work and the ASCE [11] design proposal show an excellent correlation for the entire range of hole sizes. The design equation [4.2] is derived in appendix B and was based on the ASCE proposal assuming a square plate and a circular perforation. This equation along with the results of the finite element work have been illustrated in figure [4.4] and show almost a straight line correspondence between decreasing capacity and increasing hole size. (4.2) >/4/352 -25 + 1 44 1.2 1.1-Hole Size Ratio D/b Fig. 4.4: Comparison of the Ultimate In-plane Shear Capacities of Concentrically Perforated Plate as Calculated by the Finite Element Method and the ASCE Design Proposal given by Equation [4.2] 46 4.1.3.2 3-Dimensional Elastic Buckling The variation of the ultimate elastic buckling capacity with hole size was significantly different from that of the in-plane yield capacity. For the smaller hole sizes the decrease in capacity was almost proportional to the hole size. For the medium size holes the elastic buckling capacity was much lower than the straight line correlation found in the in-plane yield capacity. When the hole diameter became greater than 0.2 6 the buckling capacity became significabtly lower than a straight line correlation. The elastic buckling coefficients for both the clamped and simply supported plate boundaries are illustrated in figure [4.5], and compared with other theoretical results by UENOYA, REDWOOD [3] and MARCO [9]. The ultimate elastic buckling capacity is related to the elastic buckling coefficient by equation [4.3]. UENOYA and REDWOOD used a combination of in-plane finite element stress analysis and a Rayleigh-Ritz energy method to determine the elastic buckling coefficient. The perforated plate surface was discretized by the finite element method using constant stress triangle (CST) elements. An elastic analysis provided the stress distribution throughout the domain. The stresses were then substituted into the minimum potential energy expression, and an bifurcation analysis provided the buckling loads. The deflected shape in the energy expression was represented by the first eight terms of a Fourier series. The (4.3) 47 resulting values showed a good correlation with the classical solution for a full plate and agree with the current work for the smaller holes. However, as the hole size, D/b, becomes greater 0.4 their elastic analysis gives a much higher buckling load than the current finite element work or that of MARCO. The elastic buckling coefficients determined by MARCO were obtained using the same program as used in the current work. A full model of the perforated plate was used throughout the analysis. The model was made up of 16 bicubic, isoparametric, plate shell elements. Despite using fewer elements in the model, the results he obtained were slightly lower than the current work. One possible explanation for this is the use of a diferent initial load increament before the bukling analysis is done. If MARCO applied a higher initial load to the plate before the bifurcation analysis was done he would get lower buckling load. Also, if a lower order integration was in his analysis it would have softening the stiffness matrix, thus a slightly lower buckling load would result. A typical buckling mode for the quarter plate model is illustrated in figure [4.6], A along the compression diagonal, on the right-hand side, the displacements tend to be concentrated in the center of the plate. Along the tension diagonal, on the left-hand side, the displacements are more evenly distributed. This is because along the compression diagonal the plate is subjected to tension across the boundary. This tension tends to restrain the out-of-plane movement of the plate in this area. Thus, lower displacements occur in these areas. 48 16 Hole Size Ratio D/b Fig. 4.5: Variation of Elastic Buckling Coefficient with Concentric Hole Size 49 Fig. 4.6: Elastic Buckling Mode, Concentric Hole 4.1.3.3 3-Dimensional Elastic-Plastic Buckling The combined effects of nonlinear material and geometry were studied by determining the ultimate elastic-plastic buckling capacity of the perforated plate. The significance of the nonlinear material behavior is displayed in the typical load-deflection path shown in figure [4.7], The plate buckling load was determined at various stages along the load-deflection path. As each new equilibrium configuration was established a bifurcation analysis was performed to determine the buckling load. When the plate material first starts to yield, the stiffness matrix is effectivly softened and the buckling capacity is reduced. As the plate move out-of-plane the post buckling strength of the plate increase the buckling capacity. Finally, the ultimate elastic-plastic buckling capacity is reached when the material yeilding has lowered the buckling load to the same value as the current load level. The resulting load deflection paths of four perforated plates with different hole sizes are shown in figure [4.8]. While the material behaves elastically the lateral displacements are small. At approximately 90% of the ultimate load there is a very rapid change in the plate stiffness. The deflection increases very quickly with only a small increase in applied load. There is a long plateau as displacements continue to increase with only a small increase in applied loading. Finally, when the nonlinearities have lowered the buckling capacity sufficiently the determinant of the stiffness matrix becomes negative. This point defines the ultimate elastic-plastic bukling capacity of the plate 51 150 140-130-120 110-100-90-80-70 Stable 4 Unstable • Legend Load Path Buckling Load I 12 14 i 16 6 8 10 Lateral Deflection mm ig. 4.7: Decreasing Elastic-Plastic Buckling Capacity of a Perforated Plate with Increasing Applied Load. 18 200 Full Pate D/b=0.15 D/b=0.2 D/b=0.3 0 2 4 6 Lateral Deflection mm Fig. 4.8: Load Deflection Curve of a Simply Supported Perforated Plate with Various Concentric Hole Sizes 8 D/b=0.5 i 1 i 1 i 1i 10 12 14 16 i 18 -I—'—l—'—I—'—l— 20 22 24 26 28 52 Although the curves have similar characteristics, figure [4.8] clearly demonstrates that the length of the plateau changes with hole size. As the hole size increase the plateau becomes longer and at the same time the slope of the load-deflection path becomes steeper. This means is that for larger holes there will be more out-of-plane distortion in the plate before the ultimate load is reached. Also, a steeper slope indicates that the plate is stiffer and will still carry more load. Thus, plates with larger holes will be more ductile. Figure [4.8] also clearly shows that the ultimate elastic-plastic buckling capacity is substantially reduced with increasing hole size. The results of the present study along with numerical results by UENOYA for the ultimate elastic-plastic buckling capacities of a concentrically perforated plate are illustrated in figure [4.9]. The results of the present study were obtained from load-deflection curves similar to those in figure [4.8]. The ultimate load was defined as the highest load level obtained on the load-deflection curve. The two sets of results show good agreement for the smaller hole size, but the current work of the author gives substantially lower values for the larger holes. There may be two possible explanations for these discrepancies. Firstly, the use of the relatively stiff CST element, used by UENOYA and REDWOOD, to calculate the stress distribution may have resulted in an underestimation of the true stress field. If the stresses were under estimated the Ralyleigh-Ritz minimization would produce a higher plate capacity. Secondly, by restricting the displacement 53 1.2 1.1 T 1 1 1 1 1 1 1 1 1— 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Hole Size Ratio D/b Fig. 4.9: Variation of Ultimate Elastic-Plastic Buckling Capacity of Simply Supported Perforated Plate with Concentric Hole Size. 54 Fig. 4.10: E-P Buckling with von-Mises Stress function to the first eight terms of the Fourier series, UENOYA assumes that the resulting displacement field may be accurately represented by a combination of these modes. This assumption was sufficient for the full plate buckling mode to which it was compared. However, if higher terms become more dominate as the hole size increases, the displacement function may not be able to accurately represent the lowest buckling mode. The result would be to over-predict the buckling capacity. 4.2 ' Variation of Hole Location The second parameter studied was the location of a hole on the plate capacity. The center of a standard hole was placed in various locations on the plate surface. For each location the ultimate in-plane, elastic buckling, and elastic-plastic buckling capacities were determined. 4.2.1 Plate Geometry •b/2 b/2 b/2 b/2 Fig. 4.11: Plate Geometry and Loading used in the Analysis of the Variation of Hole Location Parameter 56 A standard hole with a diameter of 0.2b was centered at various locations in the plate. The geometry and loading are shown in figure [4.11]. The parameters used to define the hole location, ANG. and R/b, were varied as follows; ANG=45 to 135 deg., 5/6=0.15 and 0.3. The hole only needed to be considered in one quarter of the plate area due to the symmetry in loading and geometry. 4.2.2 Finite Element Model Two basic element meshes were required to do the analysis for the eccentric hole locations. When the hole center was located along a plate diagonal there is an axis of symmetry about this diagonal. By applying the appropriate displacement boundary conditions along this axis, only a one half-plate model was required. However, if the hole was not located on a diagonal there was no axis of symmetry and a full plate model required. The half model was used whenever possible as it used significantly less CPU time for the analysis. The half and full plate models were generated by the same method. The hole center was selected as the origin of a polar coordinate system. Radial lines were then established to each corner of the plate, dividing the model into quarter sections. Each section was then modeled by using nine isoparametric elements in a 3x3 grid. This resulted in a half and full model made up of eighteen and thirty-six elements respectively. The nodal spacing in each quarter of the model was determined in a similar manner as was used for the center hole models. With the polar coordinate system each quarter was divided into nine sections by ten equally spaced radial lines. Ten nodes were then 57 CC Fig. 4.12: Finite Element Model of Half the Plate Fig. 4.13: Finite Element Model of the Total Plate placed along each radial line with proportional spacing as given by equation [4.1]. Typical examples of the resulting finite element mesh for the half and full models are illustrated in figure [4.12] and [4.13]. These element meshes were used through-out the analyses, with some distortion to accommodate the plate geometry. 4.2.3 Results 4.2.3.1 In-plane Yielding Results from the in-plane yielding analysis are illustrated in figure [4.14]. As the 0.2b diameter hole was moved over the plate surface there was little change in the ultimate in-plane yield capacity. In fact only one location had more than a 2% change. This location had the hole near the plate boundary and had a 4% change in capacity. The model defined by the parameters, 72/6 =0.3, ANG.=90, D/b =0.2, produced a significantly lower capacity than the other models. Yielding started around the hole in this model as with other models. Once sufficient yielding had occurred between the inner hole boundary and the outer plate boundary, the plate reached its ultimate capacity. The failure mechanism was similar to the other configurations but yielding was restricted to the small area between the hole boundary and the plate outer boundary. Therefore this configuration was considered a local material failure rather than an ultimate plate capacity limit. In theory this local material failure mode may be of interest, but in reality shear webs would not likely experience this type of failure. Most shear webs have a flange or stiffener located on all sides. This would provide a mechanism for 60 '6 ID •> v • 1 *' 0.962 * .•••«•.. N .•0.980 / •-.0.980 S-.^"6.989 1.03 '•0.989 / 0.995'-' ".'0595 R/b=0.3 10 ^R/b.= O.I5 \, b/2 \ b/2 . Fig. 4.14: Ultimate In-plane Yield Capacity Results for Various Hole Locations Normalize to the Concentric Hole Ultimate In-plane Capacity redistribution of forces. The areas of high strain would transfer forces into the boundary stiffeners. These forces would travel through the stiffeners and then be transferred back into the plate at regions of lower strain. 4.2.3.2 3-Dimensional Elastic Buckling Unlike the in-plane capacity, significant changes in the ultimate elastic buckling capacity occur when the hole location moves away from the center of the plate. If the hole was moved from the or tension diagonal to the compression diagonal the elastic buckling capacity could be increased by as much as 50%. A number of results for the elastic buckling capacities are illustrated in figure [4.17]. All the capacities are given as a factor of the elastic buckling capacity of a plate with a concentric hole so that direct comparisons can be made between locations. If the hole was located in the tension diagonal the 61 62 Fig. 4.16: Full Plate, von-Mises Stress factor ranged from 0.992 to 1.017. This is essentially constant. However, along the compression diagonal the factor increased to 1.52. This means that there was a 52% increase in the elastic buckling capacity by moving the hole from the tension diagonal to thecompresion diagonal. R/b=0.3 R/b = 0.!5 b/2 •b/2 Fig. 4.17: Elastic Buckling Capacity Factors for Various Hole Hole Locations Normalize to the Concentric Hole Capacity. The elastic buckling mode of the plate shown in figure [4.18] and [4.19] illustrates, the different profiles of the tension diagonal and compression diagonal respective. The displacement profile along the tension diagonal was very smooth, similar in shape to a simple sine wave. However, across the compression diagonal, the profile appeared to be made up of more complicated shapes. The displacements were more pronounced at the hole boundary then at the edges of the plate. At the corners of the plate the displacements were very small and examination of the eigenvector showed that some were actually negative in these areas. 64 4.2.3.3 3-Dimensional Elastic-Plastic Buckling The results of the calculations for the ultimate elastic-plastic buckling capacity are detailed in figure [4.20]. The capacities are expressed as a factor of the concentrically perforated plate ultimate elastic plastic capacity. The ultimate elastic-plastic buckling plate capacity for each hole location is related to the magnitude of change in the elastic buckling capacity. The results showed that if there was little or no change in the elastic buckling capacity, the plate would undergo a typical elastic-plastic buckling (zone II) failure. The combined effects of nonlinear material and geometry would produced an ultimate elastic-plastic capacity lower than if either of the two factors were consider alone. However, if the elastic buckling capacity was increased significantly, as when the hole is moved to compression diagonal, the buckling capacity was so much higher than the in-plane yield capacity that the plate would fail without buckling. Thus the plate would experience an in-plane yield failure (zone I). Examples of models that underwhen zone I and a zone II failure modes, are given in table [4.1]. Table 4.1: Failure Mode Classification Model Character i st ics Failure D/b R/b ANG. MPa. MPa. MPa. Class A 0.2 0.3 45 136.9 133.5 128.3 zone II 2C 0.2 0.0 - 138.4 131.3 126.2 zone II E 0.2 0.3 135 136.9 207.5 138.8 zone I 65 Fig. 4.18: Profile of the Tension Diagonal : Profile of the Compression Diagonal *6 • b/2 -j. b/2 Fig 4.20: Ultimate Elastic-Plastic Capacity for Various Hole Locations Normalized to the Concentric Hole Ultimate Elastic-Plastic Buckling Capacity Since the elastic buckling capacity of the plate with eccentric hole always higher than that of a plate with a concentric hole and, no examples of elastic buckling (zone III) failure were found in this analysis. The elastic buckling capacity was never significantly lower than the in-plane capacity. However, this does not mean that a zone III failure could not occur. If the slenderness ratio of the plate were decreased, the elastic buckling capacity would be reduced, while the in-plane capacity would remain unchanged, forcing the plate into a zone III failure mode. The concentrically perforated plate provides a lower bound value for all other hole locations. With reference to figure [4.20], all but one hole location produced an ultimate 68 capacity factor greater than that of the concentric hole. The hole location that resulted in a lower ultimate capacity was the same one that had a low in-plane capacity. As discussed in section [4.2.3.1] this low value was due to a local material failure. It is unlikely that this mode of failure will occur in reality, since the flanges and stiffeners around the web will provide a mechanism for the redistribution of forces. If local material failures are prevented, then the ultimate in-plane capacity, for any hole location, will be approximately equal to the in-plane yield capacity of a plate with concentric hole. The elastic buckling analysis indicated that the concentric hole had the lowest capacity. The ultimate elastic-plastic buckling capacity was governed by a combination of the in-plane yield capacity and elastic buckling capacities. Moving the hole away from the center of the plate increased the elastic buckling capacity and had little effect on the in-plane capacity. Thereforethe minimum ultimate elastic-plastic buckling capacity for any location of a given hole size is given by the concentric hole configuration. 4.3 Optimum Doubler Plate The final parameter investigated in the study was the effectiveness of the addition of a doubler plate around the perforation. The objective of the doubler plate was to lower the stresses around the hole, limiting the material yielding and restoring some or all of the plates in-plane yield capacity. The goal of this study was to evaluate the effectiveness of various doubler plate shapes. 70 4.3.1 Plate Geometry The standard square plate with a concentric 0.26 diameter circular perforation was used throughout this part of the study. A doubler plate in the circular shape was attached to the perforated plate around the hole. The diameter and thickness of the doubler plate was varied. The diameter varied from 1 .3D to 2.425Z>and the thickness from 0.25* to 1.50t. 4.3.2 Finite Element Model Like the other models for the concentrically perforated plate only the one-quarter model was required. The axes of symmetry were again the plate diagonals. The same displacement boundary conditions were applied along these axes of symmetry and all displacements were restricted to in-plane movements only. The one-quarter plate model shown in figure [4.23] consisted of an 8x5 mesh of plane stress elements discussed in section [3.1.4.2], The model was generated by dividing the quarter plate into eight sections by nine equally spaced radial lines. Six nodes were then placed on each radial line. The first node placed on each line was set at the hole boundary. The next four nodes on each radial line were set a constant radii of 1.3£J 1.75D, 2.0D and 2.425D. The last node was then located along the outer plate boundary. The doubler plate was .modeled by specifying a thicker plate for the internal element rings. For the smallest diameter doubler plates only the first element ring was thickened. For larger doubler plate sizes two or three element rings were thickened. This allowed the same model to be used throughout the analysis with minimal changes between models. 71 T JL X\ \Dd \\200 \ E = 200 000 MPa (Ty - 300 MPa v = 0.3 t = 10 500 500 Fig. 4.22: Geometry and Loading of Perforated Plate with Doubler Plate Fig. 4.23: Finite Element Mesh of Perforated Plate with Doubler Plate 72 4.3.3 Results 4.3.3.1 In-plane Yielding The doubler plate analysis showed that, given the same cross sectional areas, a wide, thin doubler plate was more effective than a narrow, thick plate. In figure [4.24] a plot of the effective capacity restoration factor vs. a nondimensional doubler plate area is shown for various doubler plate diameters. As the doubler plate diameter ratio, $ was increased there was a significant increase in the capacity restoration factor. The capacity of a plate with doubler plate is given by equation [4.4]. If the capacity restoration factor is 1.0, the plate was restored to its original unperforated capacity. In figure [4.24] it is shown that a doubler plate of Dd = 2.0D and Ad/A=1 »0 has a capacity restoration factor of almost 1.0. Solving for the doubler plate thickness * td> gives the parameters for an optimum doubler plate size as, Dd=2.0D, td=t. Uy = fy + 1> (rv ~ fy) (4.4) Without the doubler plate, yielding started at the inner boundary of the perforation arid propagated around the hole and up into the body of the plate. Once the yielding had extended from the inner perforation boundary to the outer plate edges, the plate became unstable and the ultimate load had been attaned. With the addition of a doubler plate there was a significant modification to the yield pattern. Again, yielding started at the perforation boundary, but it did not extend into the plate. Instead, a second yield zone developed at the outer doubler plate boundary. This yielding spread rapidly from the edge of the 73 doubler plate to the plate boundary as the load level increased. Once the yielding had extended inward to the inner perforation boundary there was no further increase in plate capacity. A comparison of the two ultimate yield patterns, for the standard and the reinforced perforated plates as described above -is given in figure [4.25], 0.5 1.0 Doubler Plate Area AA/A Fig. 4.24: Effective Capacity Restoration Factor vs Doubler Plate Area for Various Doubler Plate Diameters 74 doubler plate reinforcement Dd = 400 j td = 20 without reinforcement Fig. 4.25: Spread of Yield Zones for Standard and Reinforced Perforated Plates 4.1 Convergence with Mesh Refinement An approximate solution method is enhanced if it can be rigorously proven that, as the step size or element size is reduced, the method will render the exact solution. The approximate method should provide some sort of bound on the exact solution and show an asymptotic convergence to the solution. The finite element formulation satisfies all of these requirements. From finite element theory it can be shown that the strain energy of the finite element solution provides a lower bound value for the system strain energy. Furthermore it has been shown that an eigenvalue analysis of the finite element stiffness matrix will provide an upper bound value for the system's elastic buckling load. The theory also states that the finite element strain energy will asymptotically converge to the exact solution like (l/n)p . where n is the number of elements in any one direction 75 and p depends on the governing differential equation of the continuum problem and the interpolation functions used to formulate the element. Proofs of these properties can be made under certain conditions. The consistent load must be applied over the system domain as well as between elements and exact integration of the element area is assumed. The consistent load requirement between elements and on the element boundaries has been satisfied throughout the analysis. The loading between elements was automatically satisfied by the finite element formulation. Along the element boundaries the consistent shear load vector for each element subjected to uniform shear was calculated and applied. Exact integration requires a lot of CPU time if numerical integration is used. The accuracy of the integration must be sufficient to exactly evaluate the element stiffness integral, including all term in the determinant of the Jacobian coordinate transformation matrix. Terms in the Jacobian matrix may be to the second and third power and the determinant will have terms to the fourth, fifth and sixth power. Combining this with other terms in the element stiffness integral would require integration be suficient to fully evaluate a polynomial of order nine or ten. Instead of exactly integrating the stiffness integral for all cases, BATHE [10] has shown it is sufficient to use a reduced integration, however, the integration must be sufficient to exactly evaluate the stiffness integral if the determinant of the Jacobian matrix is a constant. Since exact integration of the element stiffness integral is 76 not realistic a lower order integration is used. Using a lower order integration the finite element solution may still asymptotically converge to the exact value. However, the rate of convergence will be governed by the integration error and not the order of the element. By assuming that there was an asymptotic convergence rate with mesh refinement, the accuracy of the elastic buckling capacities calculated by NISA83 were estimated. The work was performed using the standard 1/4 plate model of a simply supported perforated plate with a hole diameter of 0.26. The element mesh for the plate was varied from 1x1 to 6x6. A bifurcation analysis was done on each of these permutations to approximate the elastic buckling capacity of the perforated plate. A value was then assume for the exact elastic buckling capacity of the plate and the relative error calculated for the buckling capacity give by each mesh. A Log-Log plot of relative error vs number of elements in one direction was set up and straight line passed through the points using a least squares fit algorithm. A new exact solution was then assumed and this process was repeated. By keeping recording of the cumulative error of each least squares fit associated with an assumed exact solution, the best value for the exact elastic buckling capacity was established. The results of this work are shown in figure [4.26]. The slope of the line in figure [4.26], representing the convergents rate, is approximately -3. Therefore, the convergence rate for this element and model is n to the power -3 or (1/n)3 . 77 Also given by the figure is the relative error associated with each mesh. It shows that the 3x3 element mesh has an error of 5%. A more refined grid would have provided a more accurate solution for the ultimate elastic and elastic-plastic buckling capacities, but would require more storage and CPU time which was not available on the VAX 11/730. The full plate model with the present mesh already used the capacity of the VAX. By using a consistent element mesh, with varying degrees of distortion depending on hole size and location, the calculated ultimate capacities that express the same relative error. Therefore, results from all the analyses were directly compared. Any changes in ultimate capacities was attributed to the parameters studied and not due to changes in the modeling technique. 0.001 Number Of Elements In One Direction n Fig. 4.26: Convergencs of the Elastic Buckling Load with Mesh Refinement for a Concentrically Holed Plate with Z>/6 = 0.2, 1/6 = 0.01 78 5 CONCLUSIONS The behavior of a square shear plate with a circular perforation at its ultimate load can be described by one of three different failure modes. The parameter which determines the failure mode is the plate slenderness ratio, t/b ( thickness / width). For stocky plates the ultimate capacity is limited by the in-plane material yield capacity. For slender plates the elastic buckling capacity of the plate is much lower than the material in-plane yield capacity. The ultimate capacity is therefore controlled by the elastic buckling capacity. Finally, for intermediate slender plates, both the in-plane yield capacity and the elastic buckling capacity are of the same magnitude. The failure mode is a function of both material yielding and buckling so the ultimate capacity is governed by the elastic-plastic buckling capacity. No analytical solution exists for determining the ultimate capacity of perforated plates with this type of failure. Therefore, numerical methods are required to estimate these capacities. The program NISA83 was used to carry out a parameter study on perforated plates. For each parameter configuration the program calculated the ultimate in-plane yield capacity, elastic buckling capacity and the ultimate elastic-plastic buckling capacity. The restart option of the program and ability to change the time step control and iteration methods contributed to the efficiency of these calculations. The "constant arclength" time step control proved to be well suited to follow the load-deflection path of the plate into the post buckling region. 79 Some minor changes were required to get NISA83 running on the Civil Engineering VAX 11/730. The computer speed and storage capacity was sufficient to handle the nonlinear problem modeling a 1/4 plate. However, when the full plate model was used the CPU time requirement became large and the memory storage became critical. If a much larger nonlinear problem were attempted a larger computer would be required The plotting program, NISPLOT, was developed for the output from NISA83. Although, the program ran under the EUNICE operating system it was found to be completely compatible with the VMS output files from NISA83. The information that is presented in a graphical form for both data checks and post-processing displacements or stresses make a graphics program a necessity for any finite element program. It is reccommended that NISPLOT (or a similar program) be extended to include all the elements in the NISA83 library. The first parameter investigated in the study was the variation of a concentric hole size. The results for the ultimate in-plane yield capacity showed a straight line correlation between increasing hole size and decreasing plate capacity. The results were also correlated to the ASCE Suggested Design Guidelines. It was found that these design rules tended to overestimate the in-plane capacity of the plate. The elastic buckling capacities calculated for each parameter variation were in agreement with other published results. For holes larger than 0.4 of the plate width, the elastic buckling capacity was reduced significantly below a 80 straight line correlation. This reduction should be taken into account in the design of any web where buckling could occur. Finally, the ultimate elastic-plastic buckling capacities were calculated for each hole size. The results were compared to other published work. The current results showed that the capacity of the plate have been slightly overestimated by others. The ultimate elastic-plastic buckling capacity was lower than either the in-plane yield or the elastic buckling values. Yielding of the material around the hole reduces the plate buckling load thus lowering the ultimate plate capacity. The second parameter in the analysis was the hole location. The ultimate in-plane yield capacity was calculated for each location. The results indicate that there is little variation in the plate capacity with hole location. The capacity of a plate with a concentric hole seemed to provide a good approximation of a perforated plate even with a hole in any location except close to the plate edge. If the hole was located too close to the boundary of the plate there was the possibility that local material yielding between the hole and the near boundary could reduce the plate in-plane capacity. The calculation of the elastic buckling load for each variation in hole location yielded some unexpected results. The elastic buckling load was found to increase by up to 50% if the hole was moved from the plate tension diagonal to the compression diagonal. The concentric hole produced the lowest elastic buckling capacity of any of the hole locations. 81 The ultimate elastic-plastic buckling capacities appeared to be a combination of the first two capacities. By moving the hole away from the center of the plate the ultimate plate capacity would be governed by the in-plane yield capacity, if the elastic buckling load became significantly higher than the in-plane yield capacity. Since the concentric hole provided the lowest elastic buckling capacity it was not a surprise to find that this location also had the lowest ultimate elastic-plastic buckling capacity. Thus the plate with a concentric hole provides a lower bound value for the ultimate elastic-plastic capacity of a plates with the same hole size. The effectiveness of a doubler plate as reinforcement was also investigated. The diameter and thickness of circular doubler plates were varied. The doubler plate dimensions, *«//t = 1.0, and •0/6=2.0, gave a perforated plate the same in-plane yield capacity as it would have without the hole. An attempt was made to establish the accuracy and convergence rate of the element mesh used in all the analyses. The elastic buckling load was calculated for the one quarter plate model using a mesh ranging from 1x1 to 6x6 elements. By passing the best straight line through the Log-Log plot of relative error vs number of elements in one direction, (n), a convergence rate of (1/n)3 was determined. The error in the 3x3 element grid, used throughout this work, was estimated from this plot to be 5%. 82 REFERENCE 1 Wang, Chu-Kia. "Theoretical Analysis of Perforated Shear Webs", Presented at a meeting of the ASME Cincinnati  Section, Cincinnati, Ohio Oct. 2-3, 1945. 2 Rockey, K. C, Anderson, R. G., and Cheung, Y. K. "The Behavior of Square shear Webs Having A Circular Hole", Symp.  on Thin Walled Steel Structures, University Colledge of Swansea, Crosby Lockwood and Sons Ltd. 1969, pp.148-169. 3 Uenoya, M. and Redwood, R. G. "Elaso-PLastic Shear Buckling of Square Plates with Circular Holes", Computers and  Structures, Vol.8, pp. 291-300, Pergamon Press Ltd., 1978. 4 Uenoya, M. and Redwood, R. G. "Buckling of Webs With Openings", Computers and Structures, Vol.9, pp. 191-199, Pergamon Press Ltd., 1978. 5 Redwood, R. G., and Uenoya, M. "Critical Loads for Webs with Holes", Journal of the Structural Division, ASCE, Vol.105, No. 105, pp. 2053-2067, Oct. 1979. 6 Janssen, T. L. "A Simple Efficient Hidden Line Algorithm", Computers and Structures, Vol.17, pp. 563-571, Pergamon Press Ltd., 1983. 7 H'afner, L., Ramm, E., Sattele, J. M., and Stegmuller, H. "NISA80 Proqrammdokumentation Programmsystem", Bericht des Instituts fur Baustik, Universitat Stuttgart, 1981. 8 Precision Visuals Inc. "DI3000 Users Guide", Precision Visuals, 6260 Lookout Road, Boulder, Colorado, 80301 USA, March 1984. 9 Marco, Renzo "Buckling of Plates with Circular Holes", B.Ap.Sc. Thesis, University of Maitoba, 1984. 10 Bath^ Klaus-Jurgen "Finite Element Procedures in  Engineering Analysis", Prentice Hall Inc., Englewood Cliffs, N. J. 1982 11 Subcommitty on Beams with Web Openings "Suggested Design Guides for Beams with Web Holes", ASCE, Journal of the  Structures Div., Vol.97, pp.2707-2728, Nov. 1971. 12 Bathe^ Klaus-Jurgen and Bolourchi, Said "A Geometric and Material Nonlinear Plate Shell Element", Computers and  Structures, Vol.11, pp. 23-48, Pergamon Press Ltd., 1980. 83 APPENDIX A Derivation of the Consistent Shear Load Vector for the Bicubic Isoparametric Element. Using a higher order element such as the bicubic isoparametric element gives excellent results to a problem modeled with a small number of elements. However, using fewer elements means that each element is more distorted and hence, more dependent on the applied loads. Work done by BATHE [12] has shown that the bicubic isoparametric element is not affected as much by distortion as many of the lower order elements. Nevertheless, he does reccommend using the element consistent load vector in order to minimize any errors caused by the element distortion. The following is the derivation of the consistent load vector for the bicubic isoparametric element, with a uniform shear applied along one boundary. The load vector is defined in terms of the element coordinates. The resulting equations can be used to determine the consistent shear vector for any distorted element. Fig. A.1: Bicubic Isoparametric Plate Shell Element with Uniform Shear Loading along One Edge 84 The definition of the ith term of the element consistent load vector P is give by equation [A.1]. Pi = JJ q(x,y)Ni(r,s)dxdy (A.l) Area where <l{xiy) = applied traction on the element surface Pi = i'Aterm of consistent load vector ^(^3)= shape function for node x For the isoparametric element this can be converted into the r,s local coordinate system and the volume integration reduced to the limits of -1 to 1. +1+1 P{ = J J q (z, y) 7Y, (r, s) detJ dr ds (A.2) -1-1 where detJ = determinant of the Jacobian transformation matrix If <j(r,s) is uniform along along the edge s=1 and zero everywhere else, equation [A.2] can be reduced to. +1 Pi = q j Ni detJ dr i = 2,6,5,1 (4.3) -1 and ^ dN-. dr 3 j=2,6,5,l or converting the indices to 1 through 4, gives the line integral, +1 Pk = q f Nkiy^x^dr k = 1,2,3,4 x2 x3 85 Table A.I: Four Cubic Shape Functions along the Bicubic Element Boundary s=1, -1<r<1 TV* * 1 6 fact. (r+1) (3r+1) (3r-1 ) (1-r) Ni 1 - 1 1 1 N2 -9 1 — 1 1 N3 9 1 1 — 1 1 1 1 1 -Multiplying these out and taking the derivatives gives: ^=(^)(-3r3 + r2 + 3r-l) ^4=(^)(9r3 + 9r2-r-l) ffi.(±)(-^ + ttr+1) dN2 /-9\ , 9 . -«r = (») (-9r2+2r+3) Expanding equation [A.3] . -l -l +X3fNt9-£dr + XlJNt 9Jl±dr -1 -1 (4.4) (4.5) (4.6) 86 Each of the four terms in equation [A.6] is of the form where k dr Nk = Ak (akr* + bkr2 + ckr+d) ^ = B,(eir2 + /ir + ri Thus, the integration for these terms need only be done once, +1 +1 Nk-g~r dr = AkBi j (akr3 + bkr2 + ckr + dk) for2 + fir + gt) dr -l = 2AkBt \{h€l ~ a"fl) + (<**e< + ckfi + bk9l) + dkei ^5 3 {A.7) Substituting the corresponding values for ak,bk,cktdk,ei,fi gt from equations [A.4] and [A.5] into equation [A.6] results in four simple linear expressions for the consistent load vector for an element with uniform shear stress along one side. Px = -0.500ZJ + 0.7125x2 - 0.300i3 + 0.0875x4 ' P2 = -0.7125Z! + 0.0x2 + 1.0125z3 - 0.300x4 P3 = 0.300Z! - 1.0125x2 + O.O13 + 0.7125x4 P4 = -0.0875X! + 0.300x2 - 0.7125x3 + 0.500x4 , U.8) 87 Appendix B ASCE Suggested Design Guides for Beams with Web Holes The ultimate in-plane shear stress of a perforated web as given by equation [18] in Reference [11] is given by: f» - ('" f) vrE r» —a* (§)'(*m) For a "circular hole H = 0.9D A = 0A5D Substituting the above into equation [B.1] gives: If the substitution S = 0.90(D/6) is made in equations [B.3] and [B.4] and the two equations are combined and simplified, a single expression, equation [B.7], is obtained for the ultimate in-plane shear stress for a square plate perforated by a circular hole. a = 3.o(i-l) (B.6) fy _^(l-S)(l/S-l) T» \Jl + 3(l/S- l)2 88 Hole Size Roitio D/b Fig. B.1: Ultimate In-plane Capacity of Shear Web with a Circular Perforation as Proposed in Reference [11] 89 APPENDIX C Modification of NISA80 at U.B.C. NISA80 The following is a summary of the changes that were made to the October 83 verion of NISA80, called NISA83 on the Civil Engineering VAX 11/730, at U.B.C. PROGRAM NISA83C Change file names to suit the disk names at U. from SLF( 1 ) = SLF(2) = SLF(3) = SLF(4) = SLF(5) = SLF(6); SLF(7): ' DRA2 ' DRA2 ' DRA2 ' DRA2 ' DRA2 ' DRA2 ' DRA2 [SCRATCH [SCRATCH tSCRATCH [SCRATCH [SCRATCH [SCRATCH [SCRATCH ]NISA20, ]NISA21, ]NISA22, ]NISA2 3, ]NISA24 ]NISA25 ]NISA26 SCR* SCR' SCR' SCR' SCR' SCR' SCR' F1='DRA2:[SCRATCH]NS01.SCR' F2='DRA2:[SCRATCH]NS 0 2.SCR' F3='DRA2:[SCRATCH]NS03.SCR' F8='DRA2:[SCRATCH]NS08.SCR' F9='DRA2:[SCRATCH]NS09.SCR' to SLF(1)='NISA20.SCR* SLF(2)='NISA21.SCR' SLF(3)='NISA22.SCR* SLF(4)='NISA23.SCR' SLF(5)='NISA24.SCR' SLF(6)='NISA25.SCR' SLF(7)='NISA26.SCR' C F1='NS01.SCR' F2='NS02.SCR' F3='NS03.SCR* F8='NS08.SCR' F9='NS09.SCR' SUBROUTINE INPUT Changes have been made to the cylindrical nodal input routines so that the user can select the normal axis. NAXIS=0 y-z plane is specified in polar coordinates x is the nomal axis NAXIS=1 same as NAXIS=0 NAXIS=2 x-y plane is specified in polar coordinates z is the nomal axis from 20 READ(INP,1000) - - - - ,Z(N),KN,IT WRITE(IOUT,2002) - - - - ,Z(N),KN,IT 1000 FORMAT( ----,15,12) 2001 FORMAT( - - - - ,5X,2HIT/) 2002 FORMAT( - - - - ,I5,2X,I5) to 20 READ(INP,1000) - - - - ,Z(N),KN,IT,NAXIS WRITE(IOUT,2002) - - - - ,Z(N),KN,IT,NAXIS 90 1000 FORMAT( - - - - ,15,12,12) 2001 FORMAT( ,5X,2HIT,3X,5HNAXIS/) 2002 FORMAT( ,I 5,2X,I 5,2X15) from C C CYLINDRICAL COORDINATES C 50 DUM = Z(N) * RAD Z(N) = Y(N) * SIN(DUM) Y(N) = Y(N) * COS(DUM) to C C CYLINDRICAL COORDINATES C 50 CONTINUE IF (NAXIS.EQ.2) THEN DUM = Z(N) * RAD Z(N) = X(N) X(N) = Y(N) * COS(DUM) Y(N) = Y(N) * SIN(DUM) ELSE DUM = Z(N) * RAD Z(N) = Y(N) * SIN(DUM) Y(N) = Y(N) * COS(DUM) END IF from DX = (X(N)-X(NOLD)) / XNUM to IF (NAXIS.EQ.2) THEN DZ = (Z(N)-Z(NOLD)) / XNUM ELSE DX = (X(N)-X(NOLD)) / XNUM END IF from C C CYLINDRICAL COORDINATES C 60 ROLD = Y(NOLD) / COS(DUMOLD) RNEW = Y(N) / COS(DUM) to C C CYLINDRICAL COORDINATES C 60 IF (NAXIS.EQ.2) THENN ROLD = X(NOLD) / COS(DUMOLD) RNEW = X(N) / COS(DUM) ELSE ROLD = Y(NOLD) / COS(DUMOLD) RNEW = Y(N) / COS (DUM) END IF SUBROUTINE MAIN Change the open statements for files NPLOT and NPLOT1 from unformated to formated. The reason for this change is the EUNIC FORTRAN 77 compiler gave an error when it read from an unformated file. This problem seems to be specific to the U. B. C. VAX and is related to running VMS with a EUNIC emulator rather than native UNIX. from IF(IPG.EQ.I) OPEN(UNIT=NPLOT,- - - ,FORM=UFM) IF(MODEX.NE.O) OPEN(UNIT=NPLOT1,- - - ,FORM=UFM) to IF(IPG.EQ.I) OPEN(UNIT=NPLOT,- - - ,FORM=FM) IF(MODEX.NE.O) OPEN(UNIT=NPLOT1,- - - ,FORM=FM) 91 SUBROUTINE D3INP Change the write statements from unformated to formated. from WRITE (NPLOT)NUMEP,MN to WRITE (NPLOT,2010)NUMEP,MN 2010 FORMAT(2I5) SUBROUTINE PLN Change the write statements from unformated to formated. from WRITE (NPLOT) HED WRITE (NPLOT) NUMNP,NUMEG WRITE (NPLOT) X1,Y1,Z1 to WRITE (NPLOT,2000) HED WRITE (NPLOT,2010) NUMNP,NUMEG WRITE (NPLOT,2020) X1,Y1,Z1 2000 FORMAT (A72) 2010 FORMAT (215) 2020 FORMAT (1P,3E15.6) SUBROUTINE D3PLOT Changes where made to all write statements from unformated to formated output. f rom to 2000 2010 2015 2020 2030 WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE WRITE FORMAT FORMAT FORMAT FORMAT FORMAT NPLOT) N,13,(INODE(I),I=1,13) NPLOT) N,9,(INODE(I),1=1,9) PLOT) N,5,(INODE(I),1=1,5) NPLOT) N,4,(INODE(I NPLOT) N,4,(INODE(I NPLOT) N,4,(INODE(I NPLOT) N,4,(INODE(I NPLOT) N,3,(INODE(I NPLOT) N,3,(INODE(I = 1,4) = 1,4) = 1,4) = 1,4) = 1,3) = 1 ,3) NPLOT NPLOT NPLOT NPLOT NPLOT NPLOT NPLOT NPLOT NPLOT (215, (215, (215, (215, (215, ,2000) N,13 ,2010) N,9, ,2015) N,5, ,2020) N,4, ,2020) N,4, ,2020) N,4, ,2020) N,4, ,2030) N,3, ,2030) N,3, 2X,13(15)) 2X,9(I5)) 2X,5(I5)) 2X,4(I5)) 2X,3(I5)) ,(INODE(I),1=1,13) (INODE(I),1=1,9) (INODE(I),1=1,5) (INODE(I),1=1,4) (INODE(I),1=1,4) (INODE(I),1=1,4) (INODE(I),1=1,4) (INODE(I),1=1,3) (INODE(I),1=1,3) Note; in line 3 WIRTE(PLOT)N,- - PLOT is not a typing error, This is how it appeared in the original vertion. SUBROUTINE PLOTGEO Change the write statement from unformated to formated. f rom WRITE (NPLOT) N,IEP,(INODE(I),I=1,IEP) 92 to WRITE (NPLOT,) N,IEP, (I NODE(I),I = 1 ,1EP) 2000 FORMAT (215,2X,<IEP>,15) SUBROUTINE WRITE Change the write statements from unformated to formated. from WRITE (NPLOT)PHED WRITE (NPLOT)DSI,DS2,DS3 to WRITE (NPLOT,2060)PHED WRITE (NPLOT,2070)DS1,DS2,DS3 2060 FORMAT (A70) 2070 FORMAT (1P,3E15.6) ************************************************** NISA80.2 The following is a summary of the changes that were made to the June 84 verion of NISA80, called NISA84 on the Civil Engineering VAX 11/730, at U.B.C. SUBROUTINE FNAMES from SFL(1)='DRA2:"SCRATCH:NISA.SCR' SFL(2)='DRA2:"SCRATCH:NISA.SCR' SFL(3)='DRA2:"SCRATCH:NISA.SCR' SFL( 4 ) = 'DRA2:"SCRATCH:NISA.SCR' SFL(5)='DRA2:"SCRATCH:NISA.SCR' to SFL(1)='SCRATCH:NISA.SCR' SFL(2)='SCRATCH:NISA.SCR' SFL(3)='SCRATCH:NISA.SCR' SFL(4)='SCRATCH:NISA.SCR' SFL(5)='SCRATCH:NISA.SCR' from 2000 FORMAT( German text ) 2010 FORMAT( German text ) to 2000 FORMAT( English text ) 2010 FORMAT( English text ) SUBROUTINE OPENRF from to SLF( 1 SLF(2 SLF(3 SLF(4 SLF(5 SLF(6 SLF(7 SLF( 1 SLF(2 SLF(3 SLF(4 SLF(5 SLF(6 SLF(7 ='DRA2: ='DRA2: ='DRA2; ='DRA2: ='DRA2: ='DRA2: ='DRA2: 'SCRATCH' "SCRATCH" "SCRATCH" "SCRATCH" "SCRATCH" "SCRATCH" "SCRATCH" NISA.RN1' NISA.RN2' NISA.RN3' NISA.RN4' NISA.RN5' NISA.RN6' NISA.RN7' 'SCRATCH: 'SCRATCH: 'SCRATCH: :'SCRATCH: ••' SCRATCH: •'SCRATCH: •'SCRATCH: NISA.RN1' NISA.RN2' NISA.RN3' NISA.RN4' NISA.RN5' NISA.RN6' NISA.RN7' 93 C's have been placed in the first column of the second version of OPENRF so that it is not compiled by the U. B. C. Vax 11/730 This second version is for the Cray computer. SUBROUTINE FOPEN from CHARACTER NAME*40 to CHARACTER NAME*(*) SUBROUTINE HEDIN Modified the output heading to acknowledge that the work is being done at the U. B. C. site on the Civil Engineering Vax 11/730. SUBROUTINE INPUT Changes have been made to the cylindrical nodal input routines so that the user can select the normal axis. NAXIS=0 y-z plane is specified in polar coordinates x is the nomal axis NAXIS=1 same as NAXIS=0 NAXIS=2 x-y plane is specified in polar coordinates z is the nomal axis from 20 READ(INP,1000) - - - - ,Z(N),KN,IT WRITE(IOUT,2002) - - - - ,Z(N),KN,IT 1000 FORMAT( - - - - ,15,12) 2001 FORMAT( - - - - ,5X,2HIT/) 2002 FORMAT( - - - - ,I5,2X,I5) to 20 READ(INP,1000) - - - - ,Z(N),KN,IT,NAXIS WRITE(IOUT,2002) - - - - ,Z(N),KN,IT,NAXIS 1000 FORMAT( - - - - ,15,12,12) 2001 FORMAT( - - - - ,5X,2HIT,3X,5HNAXIS/) 2002 FORMAT( - - - - ,I 5,2X,I 5,2X15) from C C to C C c CYLINDRICAL COORDINATES 50 DUM = Z(N) * RAD Z(N) = Y(N) * SIN(DUM) Y(N) = Y(N) * COS(DUM) CYLINDRICAL COORDINATES 50 CONTINUE IF (NAXIS.EQ.2) THEN DUM = Z(N) * RAD Z(N) = X(N) X(N) = Y(N) * COS(DUM) Y(N) = Y(N) * SIN(DUM) ELSE DUM = Z(N) * RAD Z(N) = Y(N) * SIN(DUM) Y(N) = Y(N) * COS(DUM) 94 END IF from DX = (X(N)-X(NOLD)) / XNUM to IF (NAXIS.EQ.2) THEN DZ = (Z(N)-Z(NOLD)) / XNUM ELSE DX = (X(N)-X(NOLD)) / XNUM END IF from C C CYLINDRICAL COORDINATES C 60 ROLD = Y(NOLD) / COS(DUMOLD) RNEW = Y(N) / COS(DUM) to C C CYLINDRICAL COORDINATES C 60 IF (NAXIS.EQ.2) THENN ROLD = X(NOLD) / COS(DUMOLD) RNEW = X(N) / COS(DUM) ELSE ROLD = Y(NOLD) / COS(DUMOLD) RNEW = Y(N) / COS(DUM) END IF SUBROUTINE DKTM There was a compile time error bacause of the multiple declaration of the variable ICON. from COMMON /PRECIS/ NDP,ICON to COMMON /PRECIS/ NDP,NPY ************************************** NISA80.2 UPDATE The following is a summary of the changes that were made to the update of the 3D-PLATE SHELL ELEMENT installed Sept. 17 85. on Civil Engineering VAX 11/730, at U.B.C. TRANSFER FILE D3DMAIN.FOR When the Transfer file D3DMAIN.FOR was read from the IBM diskette using the program KERMIT a non ASCII character was found on line No. 1117. The character was edited from the file using the PC editor, EDLIN, and replaced by ?. Line 1112 to 1119 SUBROUTINE D3LSS (A,G,GI,IT) C Q ****************************************************** c * * C * TRANSFORM STRESS AND STRAIN LOCAL-GLOBAL * C * LOCAL 3-DIRECTION IS ZERO ??? * C *C * A ... VECTOR TO BE TRANSFORMED95 TRANSFER FILE D3DINP.FOR The same problem of a non ASCII character occured in line No. 393 of the transfer program D3DINP.for. Again the chacter was edited from the file before transfer of the file was completed to the VAX. Line 390 to 393 2020 FORMAT (1H1,15X,'E LEMENT INFORMATIO N'// 1 5X,'IEL = NUMBER OF NODES FOR THIS ELEMENT'/ 1 5X,'IPS = STRESS OUTPUT CONTROL NUMBER'/ 2 5X,'KG = NODE INCREMENT FOR GENERATION ( SECOND CARD ?)' SUBROUTINE D3STIF When NPAR(5) was selected as 1 (commplete thickness integration) the program stopped because of an error in the varriable array dimesion. Which array was never determined, however, the varriable NBO is never appears in the parrameter string in the subroutine D3DISD. from C 120 CALL D3DISD (DISD,DDISD,B,ALFN,EDIS,DC,DCA(1,1,N),NC(1,N), 1 HTET(1,1,N),IEL,MN,NBO,ND,IFORM,HHI) C to C 120 CALL D3DISD (DISD,DDISD,B,ALFN,EDIS,DC,DCA(1,1,N),NC(1,N), 1 HTET(1,1,N),IEL,MN,ND,I FORM,HHI) C 96 APPENDIX D Program Listings APPENDIX D.l NISPLOT Q ********************************************************* C C UNIX VERSION C C N I S A 8 3 C C PLOT c C THIS PROGRAM USES "DI-3000" TO PLOT THE FINITE ELEMENT GRID OUTPUT C AND THE DEFLECTED SHAPE FROM "NISA83" IN A 3-D FORM. C c p CREATE ME T A F I L E S c INPUT AND OUTPUT FILES. c GEO. INPUT FILE (FORMATED) 'INF ILE' = 1 (FROM USER) c DISP. INPUT FILE (FORMATED) 'INF ILE ' = 4 (FROM USER) c STRESS INPUT FILE (FORMATED) 'INFILE ' = 3 (FROM USER) c OUTPUT PLOTED TO 'NOUT' = 6 c OUTPUT SCRATCH FILE 'I SCR' = 7 c FULL PLATE GEO. ' ' = 8 c DEVICE TYPE 'MDEV' = 0 (METAFILES) c 'NDEV' = 1 (GRAPHICS) Q ***************************************************************************** IMPLICIT REAL*4(A-H.O-Z) COMMON / PLT / INODE(100,13,5), IEL(100,5), N(100,5). NMAX(5) COMMON / MAX / RMIN(3), RAVE(3), RMAX(3), RATIO COMMON / STR / NF(6, 16),RS(2,25),NPOINT(4, 16),FACT(16),ICOL(7 ) COMMON / SVIEW / D(3), U(3) COMMON / HEAD / PHEAD INTEGER NUMP, NUMEG LOGICAL D3STR, D3PLT, METST, PHEAD CHARACTER*45 HED. PHED CHARACTER EOF DIMENSION STRESS(16,30,5), 1X(500), Y(500), Z(500), DX(500), DY(500). DZ(500). 2RX(500), RY(500). RZ(500). STRMAX(2) NPGEO = 1 NPDIS = 4 NPSTR = 3 NOUT = 6 ISCR = 7 NPLATE= 8 MDEV = 0 NDEV = 1 PHEAD = CALL SETUP (NPGEO,NPDIS,NPSTR,NPLATE,ISCR,D3STR.D3PLT,METST) c  C SET UP THE SCREEN FOR PLOTTING c  CALL JBEGIN CALL JDINIT (NDEV) CALL JDEVON (NDEV) IF (METST)THEN CALL JFSOPN (3,0,0,'NISAPLOT.MFL') CALL JDINIT (MDEV) CALL UDEVON (MDEV) END IF CALL JASPEK ( 1 .RATIO) IF (RATIO.LT. 1) THEN CALL JVSPAC (-1.0, 1.0, -RAT 10,' RAT 10 ) ELSE IF (RATIO.GT. 1 ) THEN CALL JVSPAC (-1.0/RATIO, 1.0/RATIO, -1.0, 1.0) ELSE RATI0=1.0 END IF CALL JSETDB (0) 97 NISPLOT Listing READ IN AND PLOT FULL PLATE AND ELEMENT MODEL. IF (NPLATE .NE.O) THEN NVIEW=0 NPLOT=NPLATE CALL READNO CALL MAXMIN CALL VIEW CALL READEL CALL JOPEN CALL PLTELE CALL PLTHED CALL JCLOSE CALL JPAUSE CALL JFRAME END IF (X,Y,Z,HED,NUMP,NUMEG. 1,NPLOT.I SCR (X.Y.Z.NUMP) (NVIEW) (NUMEG.NPLOT,ISCR) (X,Y,Z.NUMEG,NVIEW) (HED) (NDEV) ISTOP) READ IN FULL MESH NODE POINTS AND ELEMENTS. NVIEW=1 NPLOT=NPGEO CALL READNO (X,Y,Z,HED,NUMP,NUMEG,1.NPLOT.ISCR.ISTOP) CALL MAXMIN (X,Y,Z,NUMP) CALL VIEW (NVIEW) CALL A PLOT SUBROUTINE TO PLOT THE NODE POINTS. (SUBROUTINE PLOTNO ) IF (METST) THEN CALL JOPEN CALL PLTHED (HED) CALL JCLOSE CALL JPAUSE (NDEV) CALL JFRAME PHEAD=.FALSE. END IF CALL JOPEN CALL PLTNOD (X,Y.Z,NUMP,NOUT) CALL PLTHED (HED) CALL JCLOSE CALL JPAUSE (NDEV) READ ELEMENT DATA AND PLOT ELEMENT (SUBROUTINE PLOTEL) CALL READEL (NUMEG,NPLOT,I SCR) CALL JFRAME CALL JOPEN CALL.PLTELE (X,Y,Z,NUMEG.NVIEW) CALL PLTHED (HED) CALL JCLOSE CALL JPAUSE (NDEV) IF D3STR .TRUE. READ IN THE STRESS FILE AND PLOT THE UNDEFLECTED SHAPE WITH A STRESS COLOR FILL IF (D3STR) THEN NPLOT=NPSTR CALL VIEW (NVIEW) CALL READST (STRESS,NMAX,STRMAX,NUMEG,NPLOT,I SCR) CALL JFRAME CALL JOPEN CALL ELESTR (X , Y.Z.STRESS,STRMAX,NUMEG.2) CALL PLTELE (X,Y,Z,NUMEG.2) CALL PLTHED (HED) CALL JCLOSE CALL LEGEND (STRMAX) CALL JPAUSE (NDEV) END IF IF D3PLT .TRUE. THAN PLOT THE DEFLECTED AND THE ORIGINAL SHAPE IN 3-D. 98 NISPLOT Listing IF (D3PLT) THEN NPLOT=NPDIS CALL READNO (DX,DY,DZ,PHED.NUMP,NUMEG,2,NPLOT.I SCR,I STOP) IF (ISTOP.LT.O) GOTO 50 DO 39 ITER=1.10 NVIEW=2 CALL VIEW (NVIEW) CALL ADDDIS (X,Y,Z. DX.DY.DZ. RX,RY,RZ,NUMP,- 1) CALL JFRAME CALL JOPEN CALL PLTELE (X.Y,Z.NUMEG.3 ) CALL PLTHED (HED) CALL JCLOSE CALL JPAUSE (NDEV) CALL ADDDIS (X,Y,Z. DX.DY.DZ, RX,RY,RZ.NUMP,1) CALL JOPEN CALL FILLEL (RX.RY,RZ,NUMEG, 1 ) CALL FILLEL (RX,RY,RZ,NUMEG,2) CALL PLTELE (RX,RY.RZ,NUMEG,2) CALL JCLOSE WRITE (NOUT,1100) READ (5, ' (A1 ) ' )EOF IF (EOF.EC'S' .OR. EOF.EC's') GOTO 50 c  C IF D3STR AND D3PLT ARE TRUE PLOT THE DEFLECTED SHAPE WITH C A STRESS COLOR FILL 39 50 IF (D3STR) THEN IF (EOF.EC'C .OR. EOF.EC'C) THEN CALL JFRAME CALL JOPEN CALL PLTELE (X,Y,Z.NUMEG,3) CALL ELESTR (RX,RY,RZ,STRESS,STRMAX,NUMEG , 1 ) CALL ELESTR (RX,RY,RZ,STRESS,STRMAX,NUMEG.2) CALL PLTELE (RX,RY.RZ,NUMEG,2) CALL PLTHED (HED) CALL JCLOSE CALL LEGEND (STRMAX) WRITE (NOUT,1100) READ (5, ' (A1 ) •• )EOF IF (EOF.EO.'S' .OR. END IF END IF CONTINUE CONTINUE END IF EOF.EO.'S') GOTO 50 CLOSE PLOT ROUTINE IF (METST) THEN CALL JDEVOF (MDEV) CALL JDEND (MDEV) END IF CALL JDEVOF (NDEV) CALL JDEND (NDEV) CALL JEND CLOSE (UNIT=NPGEO) CLOSE (UNIT=ISCR) 1100 FORMAT(/,' <RETURN> TO CONTINUE, S<RETURN> TO STOP.') STOP END Q ******************************* ************ ********************************* C SUBROUTINE VIEW c **************************************************************************+* SUBROUTINE VIEW (NVIEW) IMPLICIT REAL*4 (A-H.O-Z) COMMON / MAX / RMIN(3), RAVE(3), RMAX(3), RATIO COMMON / SVIEW / D(3).U(3) INTEGER NVIEW CALL JRIGHT (.TRUE.) CALL JVUPNT (RAVE(1),RAVE(2).RAVE(3)) 99 NISPLOT Listing IF (NVIEW.EO.1) THEN D( 1 ) = -3.0 D(2)=-3.0 D(3) = -1 .0 U( 1 ) = - 1 .0 U(2)=-1.0 U(3) = 4.0 UMIN = RMIN( 1)-RAVE( 1 ) UMAX = RMAX( 1)-RAVE( 1 ) VMIN=RMIN(2)-RAVE(2) VMAX=RMAX(2)-RAVE(2) CALL JNORML(0.0.0.0.-1.0) CALL JUPVEC(0.0,1.0.0.0) CALL JWINDO (UMIN,UMAX,VMIN.VMAX) CALL JPERSP (-10.0) ELSE IF (NVIEW.EO.O) THEN UMIN = RMIN( 1 )-RAVE(1)*0.7 UMAX=RMAX( 1 )-RAVE( 1 )*0.7 VMIN=RMIN(2)-RAVE(2)*0.7 VMAX=RMAX(2)-RAVE(2)*0.7 CALL JNORML(0.0,0.0.-1.0) CALL JUPVEC(1.0,1.0,0.0) CALL JWINDO (UMIN,UMAX,VMIN,VMAX) CALL JPERSP (-10.0) ELSE IF (NVIEW.EO.2) THEN DUM=RMAX(1)-RMIN(1) UMIN=-0.65*DUM UMAX= 0.65*DUM VMIN=-0.65*DUM VMAX= O.G5*DUM DIST=(RMAX(1)-RMIN(1))*0.90 c  100 CONTINUE WRITE (G.1010) (D(I),I=1,3) READ (5,1000,ERR=100) BX.BY.BZ IF (BX.EO.0.0 .AND. BY.EO.0.0 .AND. BZ.EO.0.0) THEN ELSE D(1)= BX D(2)= BY D(3)= BZ 110 CONTINUE WRITE (6,1020,ERR=110) (U(I),I=1,3) READ (5,1000,ERR=110) BX.BY.BZ IF (BX.EO.0.0 .AND. BY.EO.0.0 .AND. BZ.EO.0.0) THEN ELSE U(1)=BX U(2)=BY U(3)=BZ END IF END IF 1000 FORMAT (3G12.6) 1010 FORMAT (/,' NORMAL VECTOR X.Y.Z ?',3F10.3) 1020 FORMAT (' UP VECTOR X.Y.Z ?',3F10.3) c  CALL JNORML (D(1 ) ,D(2 ) ,D(3 ) ) CALL JUPVEC (U(1 ) ,U(2) ,U(3) ) CALL JWINDO (UMIN,UMAX,VMIN,VMAX) CALL JVUPLN (DIST) CALL JPERSP (DIST*-3.0) ENDIF CALL JWCLIP ( .TRUE . ) RETURN END c **************************** C SUBROUTINE SETUP Q ***************************************************************************** SUBROUTINE SETUP (NPGEO.NPDIS,NPSTR,NPLATE,I SCR, 1 D3STR.D3PLT.METST) IMPLICIT REAL*4 (A-H.O-Z) CHARACTER CHAR LOGICAL D3STR,D3PLT,METST,STAT CHARACTER * 20 INFILE 100 NISPLOT Listing c c c 100 1 10 120 GEOMETRY AND SCRATCH FILE WRITE(6,2000) CALL IFILE (INFILE,STAT) IF ( .NOT. STAT ) GOTO 100 OPEN (UNIT=ISCR,FILE='for007.dat',STATUS='scratch' ) OPEN (UNIT=NPGEO,FILE=INFILE,STATUS='old' ) IF (INFILE .EO. '2f.geo') THEN OPEN (UNIT = NPLATE,FILE='2f.plt',STATUS='old' ) REWIND NPLATE ELSE NPLATE=0 END IF REWIND ISCR REWIND NPGEO DISPLACEMENT FILE CONTINUE WRITE(6,2010) READ (5, '(A1 ) ' ,ERR=110) CHAR IF (CHAR.EO.'y' D3PLT=.TRUE . WRITE(6,2020) CALL IFILE (INFILE,STAT) IF (.NOT. STAT) GOTO 120 OPEN (UNIT=NPDIS.FILE=INFILE.STATUS='o1d') REWIND NPDIS ELSE D3PLT=.FALSE. END IF .OR. CHAR.EO.'Y') THEN C C C 130 140 STRESS FILE CONTINUE WRITE(6,2030) READ (5, '(A 1 ) ' ,ERR=130) CHAR IF (CHAR.EQ.'y' .OR. CHAR.EO.'Y') THEN D3STR=.TRUE. WRITE(6,2040) CALL IFILE (INFILE,STAT) IF (.NOT. STAT) GOTO 140 OPEN (UNIT=NPSTR ,FILE = INFILE, STATUS= 'old') REWIND NPSTR ELSE D3STR=.FALSE . END IF C C METAFILE C WRITE (6,2050) READ (5, ' (A1 ) ' ) CHAR IF (CHAR.EO.'Y' .OR. CHAR.EQ.'y') THEN METST=.TRUE. ELSE METST=.FALSE. END IF 2000 FORMAT ( 'GEOMETRIC INPUT FILE NAME?'.$) 2010 FORMAT (//.'DO YOU HAVE A DISPLACEMENT FILE Y/N ? '.$) 2020 FORMAT ( 'DISPLACEMENT INPUT FILE NAME? '.$) 2030 FORMAT (//,'DO YOU HAVE A STRESS FILE Y/N ? '.$) 2040 FORMAT ( 'STRESS INPUT FILE NAME? ',$) 2050 FORMAT (//,'DO YOU WANT TO CREATE A METAFILE Y/N? '.$) RETURN END Q *************** ******************************************* C SUBROUTINE IFILE Q *****************************************************************S SUBROUTINE I FILE(INFILE.STAT) CHARACTER*20 INFILE LOGICAL STAT *********** ******* 101 NISPLOT Listing READ(5,'(A20)') INFILE INQUIRE (FILE=INFILE,EXIST=STAT) IF (STAT.EO. .FALSE. ) THEN WRITE(6,*) ' **** ERROR **** WRITE(6,*) ' FILE DOES NOT EXIST' WRITE(6,*) ' TRY AGAIN ' END IF RETURN END c *****************************,****,**«*,»***********************«* C SUBROUTINE READNO Q ******************************************************************* SUBROUTINE READNO (X,Y,Z,HED,NUMP,NUMEG,ICORD,NPLOT,I SCR,I STOP ) IMPLICIT REALM (A-H.O-Z) INTEGER NUMP,NUMEG,ICORD,NPLOT.ISCR CHARACTER*45 HED DIMENSION X( 1 ) , Y( 1 ) ,Z(1) READ(NPLOT,2000,ERR=240,IOSTAT=ISTOP) HED WRITE(ISCR,2000) HED IF (ICORD. EQ..1) THEN READ(NPLOT.2010) NUMP,NUMEG WRITE(ISCR,2010) NUMP,NUMEG END IF DO 220 1=1, NUMP READ(NPLOT,2020,ERR = 240,IOSTAT = I STOP) X(I),Y(I),Z(I) WRITE(ISCR,2020) X( I ) ,Y(I),Z(I) 220 CONTINUE 240 CONTINUE 2000 FORMAT (A45) 2010 FORMAT (215) 2020 FORMAT (1P.3E15.6) RETURN END c ************************************************************* C SUBROUTINE READEL c ****************************************************************** SUBROUTINE READEL (NUMEG.NPLOT,I SCR ) IMPLICIT REAL*4(A-H,0-Z) COMMON / PLT / IN0DE(100,13.5), IEL(100,5), N(100,5), NMAX(5) INTEGER NUMEG,NPLOT,ISCR c  C READ NODES OF ELEMENT AND STORE IN INODE(300,13,5) C FOR ELEMNT GROUP NUM. c  DO 420 NUM=1,NUMEG READ(NPLOT,2000) NMAX(NUM),NUMEL WRITE(ISCR,2000) NMAX(NUM),NUMEL DO 410 LOOP = 1 ,NMAX(NUM) READ(NPLOT,2010) I EL(LOOP,NUM),N(LOOP,NUM), 1 ( INODE(LOOP,J,NUM),J=1,N(LOOP,NUM)) WRITE(ISCR,2010) I EL(LOOP,NUM),N(LOOP,NUM), 1 (INODE(LOOP,J,NUM),J=1,N(LOOP,NUM ) ) 410 CONTINUE 420 CONTINUE WRITE (ISCR,2020) 2000 FORMAT (215) 2010 FORMAT (21 5,2X, 13( 115, : ) ) 2020 FORMAT (' ** COMPLETED READING IN NODE AND ELEM. DATA. **'./) RETURN END c ****************************************************************** C SUBROUTINE READST C ****************************************************************** SUBROUTINE READST (STRESS,NMAX,STRMAX,NUMEG,NPLOT.I SCR) CHARACTER*45 SHED DIMENSION STRESS( 16,30,5),NMAX(5 ) ,STRMAX(2) STRMAX(1)=10000.0 STRMAX(2)=-10.0 READ (NPLOT.1020)SHED DO 720 NUM=1,NUMEG DO 710 IEL=1,NMAX(NUM)/5 102 NISPLOT Listing READ (NPLOT, 1000) (STRESS(I,I EL,NUM),I = 1 , 16) WRITE (ISCR. 1000) (STRESS(I,I EL.NUM),I = 1, 1G) DO 700 1 = 1 , 1G IF (STRESS(I. I EL,NUM) .LT.STRMAX( 1 ) ) 1 STRMAXf1)=STRESS(I,IEL.NUM) IF ( STRESS(I.I EL,NUM).GT.STRMAX(2) ) 1 STRMAX(2)=STRESS(I.IEL.NUM) 700 CONTINUE 710 CONTINUE 720 CONTINUE IF ((STRMAX(2)-STRMAX(1)) .LT. 0.001 ) STRMAX(2)=STRMAX(2)+1.0 1000 FORMAT (4(4(2X, 1PE 1 2 . 5 ) /)) 1020 FORMAT (2X.A45) RETURN END C SUBROUTINE MAXMIN SUBROUTINE MAXMIN (X,Y,Z,NUMP) IMPLICIT REAL*4 (A-H.O-Z) COMMON / MAX / RMIN(3), RAVE(3), RMAX(3), RATIO DIMENSION X(500),Y(500),Z(500),R(3) DATA R / 1.0,1.0,1.0 / RMIN(1)=X(1) RMIN(2)=Y(1 ) RMIN(3)=Z( 1 ) RMAX(1 )=X( 1 ) RMAX(2)=Y( 1 ) RMAX(3) = Z( 1 ) DO 510 I=2,NUMP IF (X(I) .LT .RMIN(1)) RMIN(1) = X(I IF (X(I).GT .RMAX( 1 ) ) RMAX(1) = X(I IF (Y(I).LT .RMIN(2)) RMIN(2) = Y(I IF (Y(I).GT .RMAX(2) ) RMAX(2) = Y(I IF (Z(I ) .LT .RMIN(3)) RMIN(3) = Z(I IF (Z(I).GT .RMAX(3) ) RMAX(3) = Z(I CONTINUE DELTX = RMAX( 1 ) -RMIN(1 ) DELTY =RMAX(2) -RMIN(2) DELTZ =RMAX(3) -RMIN(3) IF ((DELTY/RATIO) .GT. DELTX) THEN DELT=DELTY R(1)=1.0/RATIO ELSE DELT=DELTX R(2)=RATI0 ENDIF DO 520 1=1,3 RAVE(I) = (RMAX(I) + RMIN(I ) )/2 .0 RMAX(I)=RAVE(I) + R(I)*DE LT*0.59 RMIN(I)=RAVE(I ) - R(I)*DELT*0.59 520 CONTINUE RETURN END C A**************.*.********************.********.** C SUBROUTINE ADDDIS SUBROUTINE ADDDIS (X,Y,Z,DX.DY,DZ,RX.RY , RZ , NUMP.NCASE) IMPLICIT REAL*4 (A-H.O-Z) COMMON / MAX / RMIN(3), RAVE(3), RMAX(3), RATIO DIMENSION X(500),Y(500),Z(500) , NCASE = -2 NO SCALING IS DONE. NCASE = -1 ADD ONLY CONST. TO Z DISPLACEMENTS. NCASE = 0 SCALE THE Z-DISPLACEMENT (DZ) AND ADD TO ORIGINAL COODINATES. NCASE = 1 DO BOTH THE ABOVE. 1DX(500),DY(500),DZ(500), RX(500),RY(500),RZ(500) CONST = NCASE * 0.065 * (RMAX( 1)-RMIN( 1 ) ) IF (NCASE .EO. -1) THEN C0NST=16.0 * CONST DMULT=0.0 ELSE IF (NCASE .EQ. -2) THEN 103 c c c c c c NISPLOT Listing 600 CONST=0.0 DMULT=0.0 ELSE SEARCH THROUGH THE Z DISPLACEMENTS AND THE SCALE THEM SO THAT THEY HAVE APPROX. 20% SCREEN WINDOW. RMAX(3)=DZ( 1 ) RMIN(3)=DZ( 1 ) DO 600 I=2,NUMP IF (DZ(I).GT.RMAX(3)) RMAX(3 ) = DZ( I ) IF (DZ(I).LT,RMIN(3 ) ) RMIN(3)=DZ(I ) CONTINUE IF (RMAX(3).LE,RMIN(3)) THEN DMULT=1 .0 ELSE DMULT= 0.25 * (RMAX(1)-RMIN(1 ) ) /(RMAX(3)-RMIN(3)) END IF END IF DO 610 1 = 1 ,NUMP RX(I ) =X(I)+DX(I) RY(I)=Y(I)+DY(I) RZ(I)=Z(I) + DZ(I)*DMULT + CONST 610 CONTINUE RETURN END Q *********************** C SUBROUTINE PLTHED C ***************************************************************************** SUBROUTINE PLTHED (HED) IMPLICIT REAL*4(A-H.O-Z) ' COMMON / MAX / RMIN(3), RAVE(3), RMAX(3), RATIO COMMON / HEAD / PHEAD CHARACTER*45 HED. NEWHED CHARACTER*1 FLAG LOGICAL PHEAD DIMENSION WX(4),WY(4),WZ(4) c  C LOCATE THE VIEWPLANE AND THE SET THE TEXT ATTRIBUTES C SO THE HEADING WILL APPEAR AT THE TOP OF THE PAGE. IF (PHEAD .EO. IF (RATIO .LE. CALL JCONVW JCONVW JCONVW JCONVW JCONVW JCONVW JCONVW JCONVW .FALSE . ) RETURN 1.0) THEN (-1 .0,RAT 10,WX(1),WY( 1 ) ,WZ( 1 ) ) (O.O.RATI0,WX(2),WY(2),WZ(2) ) (0.0,0.0,WX(3),WY(3),WZ(3) ) (1.0,-RATI0,WX(4),WY(4),WZ(4)) (- 1 ,0/RATIO, 1 .0,WX(1),WY(1),WZ(1)) (0.0, 1 .0,WX(2),WY(2),WZ(2) ) (0.0,0.0,WX(3),WY(3),WZ(3) ) (1,0/RATIO,-1.0,WX(4),WY(4).WZ(4)) CALL CALL CALL ELSE CALL CALL CALL CALL END IF CXBASE=WX(2)-WX(1) CYBASE=WY(2)-WY(1) CZBASE=WZ(2)-WZ(1) CXPLAN=WX(1)-WX(3) CYPLAN=WY(1)-WY(3) CZPLAN = WZ( 1 )-WZ(3) CXSIZE=0.065*RATI0*SQRT(CXBASE*CXBASE+CYBASE*CYBASE+CZBASE*CZBASE) CYSIZE=0.055*SQRT(CXPLAN*CXPLAN+CYPLAN*CYPLAN+CZPLAN*CZPLAN) CALL JUPDAT WRITE (6,2000) 0 READ (5,1000) FLAG IF (FLAG.EO.'Y' .OR WRITE (6,2001) . READ (5,1001) NEWHED END IF CALL JBASE (CXBASE.CYBASE, CALL JPLANE(CXPLAN,CYPLAN, CALL JSIZE (CXSIZE.CYSIZE) CALL JCOLOR (0) CALL JJUST (2,3) FLAG.EO.'y') THEN CZBASE) CZPLAN) 1 04 NISPLOT Listing CALL J3MOVE (WX(2),WY(2).WZ(2 ) ) CALL JFONT (18) CALL JFATTR(1,1.0.1.3,16383) IF (FLAG.EO.'Y' .OR. FLAG.EO.'y') THEN CALL JFSTRG (NEWHED) ELSE CALL JFSTRG (HED) END IF C C WRITE FOOTNOTE AT BOTTOM OF PAGE C CXSIZE=CXSIZE/2 .0 CYSIZE=CYSIZE/2 .0 CALL JLWIDE(SOOO) CALL ddUST (3,1) CALL JCOLOR (2) CALL JSIZE (CXSIZE.CYSIZE) CALL J3M0VE (WX(4),WY(4),WZ(4)) CALL JFONT (1) CALL J3STRG ('U. B. C. CIVIL ENGINEERING') 1000 FORMAT (A1 ) 1001 FORMAT (A45) 2000 FORMAT (/,' Do you want a new title? y/n') 2001 FORMAT (/,' Enter new Title',/, 1 ' B 1 2 3 4 E' ) RETURN END c *****»****************************^ C SUBROUTINE PLTNOD C »*****«».»*•.******•«*..««**««*«****»**.»*«.*************•**.*.,*•* SUBROUTINE PLTNOD (X.Y.Z.NUMP) IMPLICIT REAL*4(A-H.O-Z) CHARACTER*1 YES COMMON / MAX / RMIN(3), RAVE(3). RMAX(3) , RATIO INTEGER NUMP DIMENSION X(500),Y(500),Z(500) CALL JCMARK (2) CALL JCOLOR (0) DO 720 1=1,NUMP CALL J3MARK (X(I ) .Y(I ) .Z(I ) ) 720 CONTINUE c  C NUMBER THE NODE POINTS. c  CALL JUPDAT WRITE(6,2000) READ (5,1000) YES IF (YES.EO.'Y' .OR. YES.EO. 'y ' )CALL NODNUM (X.Y.Z.NUMP) 1000 FORMAT (A1) 2000 FORMAT (/,' NODE NUMBERING? y/n ',$) RETURN END c **************************************************************** C SUBROUTINE PLTELE c **************************************************************** SUBROUTINE PLTELE (X,Y,Z.NUMEG,NWRITE) IMPLICIT REAL*4(A-H.O-Z) COMMON / PLT / INODE(100,13,5). IEL(100.5). N(100.5), NMAX(5) COMMON / MAX / RMIN(3),RAVE(3) ,RMAX(3 ) , RATIO INTEGER NUMEG DIMENSION X(500).Y(500),Z(500) c  C CALL ELEMENT PLOT SUBROUTINE TO C PLOT OUTSIDE OF THE ELEMENT IN SOLID LINE, C AND PLOT ANY INTERNAL LINES WITH DASHED LINES. C LOOP OVER ELEMENT GROUPS. c  ICOLOR = 1 IF (NWRITE .EO. 3) IC0L0R=2 DO 810 NUM=1,NUMEG IF (NWRITE .EO. 1) IC0L0R=NUM IF (ICOLOR GE. 3) I COLOR = ICOLOR+1 CALL JCOLOR(ICOLOR) 105 NISPLOT Listing IELO=0 OO 800 K =1,NMAX(NUM) IF (IELO.EO.IEL(K.NUM)) THEN ISTYL=3 ELSE ISTYL=0 IELO=IEL(K,NUM) END IF IF (NWRITE.E0.1 .OR. NWRITE.E0.2 .OR. ISTYL.EQ.O) THEN CALL JLSTYL (ISTYL) CALL ORAY (X.Y,Z,K,NUM) END IF 800 CONTINUE 810 CONTINUE c  C WRITE ELEMENT NO. IN THE MIDLE OF THE ELEMENT. c  CALL ELENUM (X , Y,Z,NUMEG,NWRITE) RETURN END Q ****************************************^ C SUBROUTINE DRAY Q **************************** SUBROUTINE DRAY (X,Y,Z,K,NUM) IMPLICIT REAL*4(A-H.O-Z) COMMON / PLT / INODE(100. 13,5) , IEL(100,5). N(100,5), NMAX(5) INTEGER NRAY.K.NUM DIMENSION X(500),Y(500),Z(500), 1 XARRAY(12),YARRAY(12),ZARRAY(12) NRAY=N(K,NUM)-1 DO 300 J=2,N(K,NUM) XARRAY(J- 1 ) = X(INODE(K,J.NUM)) YARRAY(J- 1 ) = Y(IN0DE(K.J.NUM)) ZARRAY(J-I) = Z(INODE(K,J,NUM)) 300 CONTINUE CALL J3M0VE(X(IN0DE(K, 1.NUM)),Y(INODE(K, 1.NUM) ), 1 Z(INOD E(K, 1 .NUM))) CALL J3P0LY(XARRAY,YARRAY,ZARRAY,NRAY) RETURN END Q ***************************************************************************** C SUBROUTINE NODNUM Q ***************************************************************************** SUBROUTINE NODNUM (X.Y.Z,NUMP) IMPLICIT REAL*4(A-H,0-Z) COMMON / MAX / RMIN(3),RAVE(3),RMAX(3), RATIO CHARACTER*3 CHAR INTEGER NUMP DIMENSION X(500),Y(500),Z(500),WX( 2),WY ( 2) , WZ( 2 ) IF (NUMP.GE.200) RETURN C CALL JCONVW (0.0.0.0,WX(1),WY(1),WZ(1)) CALL JCONVW (0.5,0.0.WX(2),WY(2) ,WZ(2)) CXBASE=WX(2)-WX(1) CYBASE=WY(2)-WY(1) CZBASE=WZ(2)-WZ(1) CALL JBASE (CXBASE.CYBASE.CZBASE) CALL JPATH (1) XSIZE=0.013*(RMAX(1)-RMIN(1)) YSIZE=0.013*(RMAX(2)-RMIN(2)) CALL JSIZE(XSIZE,YSIZE) CALL JCOLOR (1) CALL JJUST (3,3) DO 700 1=1,NUMP WRITE(CHAR, ' ( 13) ' ) I CALL J3M0VE(X(I),Y(I),Z(I)) CALL J3STRG(CHAR) 700 CONTINUE RETURN END 106 NISPLOT Listing C SUBROUTINE ELENUM c «««**«**.«.*****.**.*.****^^ SUBROUTINE ELENUM (X,Y,Z.NUMEG,NWRITE) IMPLICIT REAL*4(A-H.O-Z) COMMON / PLT / INODE(100, 13,5) , IEL(100,5), N( 100.5), NMAX(5) COMMON / MAX / RMIN(3),RAVE(3),RMAX(3), RATIO CHARACTER*3 CHAR DIMENSION X(500).Y(500),Z(500).WX(3),WY(3),WZ( 3) IF (NWRITE.EO.2 .OR. NWRITE.EO.O) THEN RETURN END I F C CALL JCONVW (0.0,0.0,WX(1),WY( 1 ) ,WZ( 1)) CALL JCONVW (0.5,0.0,WX(2),WY(2),WZ(2)) CALL JCONVW (0.0,0.5,WX(3),WY(3),WZ(3)) CXBASE=WX(2)-WX(1) CYBASE=WY(2)-WY(1) CZBASE=WZ(2)-WZ(1) CXPLAN=WX(3)-WX(1) CYPLAN=WY(3)-WY(1) IF (CXPLAN.EO.O .AND. CYPLAN.EO.O) RETURN C CZPLAN=0.0 CALL JBASE (CXBASE.CYBASE.CZBASE ) CALL JPLANE(CXPLAN,CYPLAN,CZPLAN) CALL JPATH ( 1 ) CALL JJUST (2,2) XSIZE=0.022*(RMAX(1)-RMIN(1)) YSIZE=0.022*(RMAX(2)-RMIN(2)) CALL JSIZE(XSIZE,YSIZE) IC0L0R=2 DO 840 NUM=1,NUMEG IF (NWRITE ,E0. 1) ICOLOR=NUM IF (ICOLOR .GE. 3) ICOLOR=ICOLOR+1 CALL JCOLOR (ICOLOR) IELO=0 DO 830 K=1,NMAX(NUM) IF (IELO.NE.IEL(K,NUM)) THEN IELO=IEL(K,NUM) RELX=0.0 RELY=0.0 RELZ=0.0 NN=N(K,NUM)-1 DO 820 1 = 1 ,NN RELX=RELX+X(INODE(K,I.NUM)) RELY=RELY+Y(INODE(K,I,NUM)) RELZ=RELZ+Z(INODE(K,I,NUM)) 820 CONTINUE RELX=RELX/NN RELY=RELY/NN RELZ=RELZ/NN WRITE(CHAR.'(13)') IEL(K.NUM) CALL J3M0VE (RELX,RELY,RELZ) CALL J3STRG (CHAR) END IF 830 CONTINUE 840 CONTINUE RETURN END c »*«***»*******•****************.***»*^^ C SUBROUTINE FILLELE c *********************».***^^ SUBROUTINE FILLEL (X,Y,Z,NUMEG,NSURF) IMPLICIT REAL*4(A-H.O-Z) LOGICAL VISBLE. OK COMMON / PLT / INODE( 100, 13. 5) , IELC100.5), N( 100.5). NMAX(5) DIMENSION NOD(9,4).X( 1) , Y ( 1) . Z ( 1 ).I NUMB( 16 ) .DX(4 ) .DY(4).DZ(4) DATA NOD / 1, 2, 3, 5, 16. 13. 11, 14, 15, 1 2, 3, 4. 6.15.14.14.15. 6, 2 13,16, 5.15,14.11, 9. 8, 7 , 3 12. 13, 16, 16, 13. 12. 10. 9. 8/ IF (NSURF.EO.1) THEN 107 NISPLOT Listing IC0L0R=4 ELSE IC0L0R=6 ENDIF INTEN =16384 CALL JPINTR ( 1 ) CALL JCOLOR (ICOLOR) CALL JPIDEX (ICOLOR,INTEN) DO 950 NUM=1,NUMEG C C COLLECT ALL THE NODE NUMBERS FOR ONE ELEMENT IN TO ONE STRING. C DO 940 LEL=1.NMAX(NUM) ,5 DO 900 1=1,12 INUMB(I)=INODE(LEL,I,NUM) 900 CONTINUE DO 910 1=2,3 INUMB(1 + 11) = INODE(LEL+1 ,I,NUM) INUMB(I+13)=IN0DE(LEL+2,I,NUM) 910 CONTINUE C C START FILLING ELEMENT C DO 930 1=1,9 DO 920 J=1 ,4 DX(J)=X(INUMB(NOD(I, J) ) ) DY(J)=Y(INUMB(NOD(I,J) ) ) DZ(J)=Z(INUMB(NOD(I,J) ) ) 920 CONTINUE OK = VISBLE (DX,DY,DZ,NSURF) IF (OK) CALL J3PLGN (DX.DY.DZ,4) 930 CONTINUE 940 CONTINUE 950 CONTINUE CALL JPINTR (0) RETURN END Q ************************** C SUBROUTINE ELESTR C ****************************************************************** *'* ********* SUBROUTINE ELESTR (X,Y,Z,STRESS,STRMAX,NUMEG,NSURF) IMPLICIT REAL*4(A-H,0-Z) LOGICAL VISBLE, OK COMMON / PLT / INODE(100. 13,5), IEL(100,5), N( 100,5), NMAX(5) COMMON / STR / NF(6, 16),RS(2,25),NPOINT(4, 16 ) ,FACT( 16), ICOL(7 ) DIMENSION X(1).Y(1),Z(1),ELX(4),ELY(4).ELZ(4) DIMENSION XNEW(25),YNEW(25),ZNEW(25) DIMENSION STRESS(16,30,5), STRMAX(2) CALL DATIN IC0LD=O INTEN =16384 IF (NSURF.EO.1) THEN IVALUE=8 ELSE IVALUE= 1 END IF CALL JPINTR (IVALUE) C C LOOP OVER ALL ELEM. GROUPS. AND EACH ELEMENT IN THE GROUP. C DO 890 NUM=1.NUMEG DO 880 NE LE =1,NMAX(NUM ) / 5 LEL= (NELE-1)*5 + 1 C C CALL SHAPE TO GENERATE A 25 NODED ELEMENT WITH THE C SAME OUTER BOUNDARIES C AS THE 16 NODE ISPARAMETRIC ELEM. C CALL SHAPE(X,Y,Z, XNEW,YNEW,ZNEW,LEL,NUM) C C FILL THE 16 SECTIONS OF THE ELEMENT WITH THE COLOR ASSOCIATED WITH C THE STRESS LEVEL IN EACH SECTION. C 108 NISPLOT Listing DO 8GO IAREA=1,16 ICNEW=6* (STRESS(IAREA,NELE,NUM)-STRMAX(1) +0.05) / 1 (STRMAX(2)-STRMAX(1)) + 1 IF (ICNEW.NE.ICOLD) THEN IF (ICNEW.GT.7) ICNEW=7 IF (ICNEW.LT.1) ICNEW=1 ICOLD=ICNEW CALL JCOLOR (ICOL(ICNEW)) CALL JPIDEX (ICOL(ICNEW).INTEN) END IF C C START FILLING ELEMENT C DO 850 J= 1 .4 ELX(J)=XNEW(NPOINT(J, IAREA) ) ELY(J)=YNEW(NPOINT(J,IAREA)) ELZ(J)=ZNEW(NPOINT(J,IAREA)) 850 CONTINUE OK = VISBLE(ELX,ELY.ELZ,NSURF) IF (OK) CALL J3PLGN (ELX,ELY,ELZ,4) 860 CONTINUE C C END LOOP OVER ELEMENTS IN GROUP, AND END LOOP OVER ALL GROUPS. C 880 CONTINUE 890 CONTINUE CALL JPINTR (0) RETURN END Q *************************** C SUBROUTINE SHAPE Q ***************************************************************************** SUBROUTINE SHAPE( X,Y,Z, XNEW,YNEW,ZNEW, LEL.NUM) COMMON / STR / NF(6, 16 ) ,RS(2,25 ) .NPOINT(4,16),FACT(16),ICOL(7) COMMON / PLT / INODE( 100, 13,5) . IEL(100.5), N( 100,5). NMAX(5) DIMENSION X( 1 ) ,Y( 1 ) ,Z( 1) DIMENSION XNEW(25), YNEW(25). ZNEW(25), INUMB(16) DIMENSION A(8) ,SF( 16) DO 901 1=1,25 XNEW(I)=0.0 YNEW(I)=0.O ZNEW(l)=0.0 • 901 CONTINUE C C COLLECT ALL THE NODE NUMBERS FOR ONE ELEMENT IN TO ONE STRING. C DO 900 1=1,12 INUMB(I)=INODE(LEL,I,NUM) '• 900 CONTINUE DO 9 10 1=2,3 INUMB(1+11)=INODE(LEL+1,I,NUM) INUMB(I+13) = INODE(LE L + 2,I,NUM) 910 CONTINUE C C CALCULATE NEW COORDINATES OF A 25 NODED ELEMENT. C DO 940 ICOORD=1,25 C C FIND THE SHAPE FUNCTIONS GIVE THE LOCAL COORDINATES R AND S. C R=RS(1.ICOORD) S=RS(2,ICOORD) A( 1 ) = ( 1 + R) A(2) = (3*R+1 ) A(3) = (3*R-1 ) A(4)=(1-R) A(5)=(1+S) A(6) = (3*S+1 ) A(7) = (3*S-1 ) A(8)=(1-S) DO 920 J=1, 16 SF(J)=A(NF(1,J)) * A(NF(2,J)) * A(NF(3.J)) * 1 A(NF(4,J)) * A(NF(5,J))*A(NF(6,J))*FACT(J)/256.0 109 NISPLOT Listing 920 CONTINUE DO 930 1=1,16 XNEW(ICOORD)=XNEW(ICOORD) + SF(I ) *X(INUMB(I ) ) YNEW(ICOORD)=YNEW(ICOORD) + SF(I)*Y(INUMB(I ) ) ZNEW(ICOORD ) = ZNEW(ICOORD) + S F(I)* Z(INUMB(I)) 930 CONTINUE 940 CONTINUE RETURN END Q **************************************************** C SUBROUTINE LEGEND Q ***************************************************************************** SUBROUTINE LEGEND (STRMAX) COMMON / STR / NF(6, 16),RS(2,25),NPOINT(4, 16 ) ,FACT( 16 ) ,ICOL(7 ) DIMENSION STRMAX(2),RELX(4),RELY(4),RELZ(4) CHARACTER*17 CHAR C DATA RELX / -0.01, -0.10, 0.0, 0.10 / DATA RELY / -0.025, 0.0. 0.05, 0.0 / DATA RELZ / 0.0, 0.0, 0.0. 0.0 / CALL JRESET CALL JRIGHT ( .TRUE . ) CALL JVPORT ( .45,1.0,-1.0,1.0) CALL JVUPNT (0.0,0.0,0.0) CALL JNORML (0.0,0.0,-1.0) CALL dUPVEC (0.0,1.0,0.0) CALL dWINDO (-.275, .275. -1.0, 1.0) CALL dPERSP (-1.0) CALL JWCLIP ( .FALSE . ) CALL dOPEN CALL JPINTR (1) CALL JSIZE (0.04,0.04) CALL JdUST (1,2) CALL JCOLOR (0) CALL J3M0VE (-0.275,0.87,0.0) CALL JHSTRG ( ' [BUND]STRESS*[BLC]T[ELC] M[BLC]PA*MM' ) CALL dSIZE (0.025.0.025) PNTX = -0.160 DO 100 1=1, 7 PNTY = 0.87-1*0.075 51 = STRMAX(1) + (1-1)*(STRMAX(2)-STRMAX(1)) /6.0 52 = STRMAX(1) + (I)*(STRMAX(2)-STRMAX( 1 ) ) /6.0 - 0.1 IF (I.E0.7) THEN WRITE(CHAR,1010) S1 ELSE WRITE(CHAR,1000) S1,S2 END IF CALL dPIDEX ( ICOL(I), 1500) CALL d3M0VE (PNTX,PNTY,O.0) CALL dR3PGN (RELX,RELY,RELZ,4) CALL d3STRG (CHAR) 100 CONTINUE CALL dPINTR (0) CALL dCLOSE C 1000 FORMAT (F7.1,' T0',F7.1) 1010 FORMAT (F7.1) RETURN END c ******************************** C LOGICAL FUNCTION VISBLE Q **************************************************************** LOGICAL FUNCTION VISBLE (ELX.ELY,ELZ,NSURF) DIMENSION ELX(4) ,ELY(4),ELZ(4 ) ,VX(4),VY(4),VALUE(2) C C SEE IF THE PLANE DEFINED BY THE FOUR PASSED POINTS IS C VISIBLE UNDER THE CURRENT VIEWING TRANSFORMATION C C NSURF =1 FILL UNDERSIDE (-Z) C =2 FILL TOPSIDE (+Z) C 110 NISPLOT Listing DO 100 1=1,4 CALL UCONWV (ELX(I ) ,ELY(I),ELZ(I ) ,VX(I ) .VY(I ) ) 100 CONTINUE C DO 110 1=1,2 DDX1 = VX(I+1) - VX(I) DDY1 = VY(1 + 1 ) - VY(I) DDX2 = VX(I+2) - VX(I+1) DDY2 = VY(I+2) - VY(I+1) VALUE(I) = -DDX1*DDY2 + DDX2*DDY 1 IF (NSURF.EO.2) VALUE(I) = VALUE(I) * -1.0 110 CONTINUE C VISBLE = VALUE(1) .GT. 0.00 .OR. VALUE(2) .GT. 0.00 C RETURN END Q ***************************************************** C SUBROUTINE DATIN Q **************************************************************************** SUBROUTINE DATIN COMMON / STR / NF(6.16),RS(2,25),NPOINT(4,16),FACT(16),ICOL(7) DATA NF / 1,2,3,5,6,7, 1,2,4.5.6,7, 1.3,4.5.6.7, 1 2,3,4,5,6,7, 2,3,4,5,6,8, 2,3.4,5.7,8, 2 2,3,4,6,7,8, 1.3,4,6,7,8, 1,2.4.6.7.8, 3 1,2,3,6,7,8, 1,2,3,5,7,8, 1,2,3,5.6,8, 4 1.2,4,5.6,8. 1.2.4.5.7,8. 1,3.4.5.7.8. 5 1,3,4.5,6,8^/ DATA RS / 1.0,1.0. 6.5,1.0, 0.0.1.0, -0.5,1.0, -1.0,1.0, 1 -1.0,0.5, -1.0,0.0, -1.0,-0.5, -1.0,-1.0, -0.5,-1.0, 2 0.0,-1.0, 0.5,-1.0, 1.0,-1.0, 1.0,-0.5, 1.0.0.0, 3 1.0,0.5, 0.5,0.5, 0.0,0.5, -0.5,0.5, -0.5,0.0, 4 0.0,0.0, 0.5,0.0, -0.5,-0.5, 0.0,-0.5, 0.5,-0.5 / DATA FACT / 1.0, 9.0, -9.0, 1.0, 9.0, -9.0, 1.0, -9.0, 1 9.0, 1.0, -9.0, 9.0, 81.0,-81.0, 81.0,-81.0 / DATA NPOINT / 8,9,10,23, 7,8,23,20, 6,7,20,19, 5,6,19,4, 1 23,10,11,24, 20,23,24,21, 19,20,21,18, 4,19,18,3, 2 24,11,12,25, 21,24,25,22, 18.21.22,17, 3,18,17,2. 3 25.12.13.14. 22.25,14,15, 17,22,15,16, 2,17,16,1 / DATA ICOL / 4,6,2.3,5,1,7 / RETURN END 1 1 1 APPENDIX D.2 MESHGEN c ********************************************************* c C PROGRAM TO GENNARATE NODE GRID C C This version generates a grid of elements for a holed plate C A 1/4, 1/2, or fullplate model can be generated if the no. C of sides specified (NSS) is 1,2,or 4 respectfully. C Fixed boudaries can be specified if NSS is negative. C C NSS = NO. OF SIDES C NER = NO. OF ELEMENTS RADIALLY C NEA(4) = NO. OF ELEMENTS PER ARC. C NN(4) = CORNER NODE NO. C c ********************************************************************** DIMENSION NL(36),RL(36),AL(36),ANG(5),XY(5) LOGICAL FLAG COMMON NEA(4),NN(4),NSS,NER,NNA,NNR,FLAG FLAG=.TRUE. INPUT=1 IOUT=7 PI=3. 141592654 DEG=PI/180.0 NUMEL=0 OPEN (UNIT=INPUT,FILE='NODE.IN' ,STATUS= ' OLD' ) OPEN (UNIT= IOUT,FILE ='NODE.OUT' ,STATUS='NEW ) READ (INPUT,1000) NSS,NER,R1 IF (NSS.LT.O) THEN NSS=-1*NSS FLAG=.FALSE. END IF DO 90 ISIDE=1,NSS READ (INPUT, 1010) NEA(I SIDE) READ (INPUT, 1020) ANG(ISIDE),XY(ISIDE ) NUMEL = NUMEL + NEA(ISIDE) 90 CONTINUE IF (NSS.E0.4) THEN NUMNO = NUMEL * 3 * (NER*3+1) ELSE NUMNO = (NUMEL*3+1) * (NER*3+1) END IF WRITE (I0UT.2OOO) NUMNO,NER ANG(NSS+1) = ANG( 1 ) + 90.0*NSS IF (NSS.E0.2) THEN XY(NSS+1) = XY(1)-1000.0 ELSE XY(NSS+1)=XY(1). END IF NN(1)= NEA(1)*3 IF (NSS.NE.1) THEN DO 95 1=2,NSS NN(I)= NN(I-1) + NEA(I)*3 95 CONTINUE END IF DO 110 IS=1,NSS NNR=NER*3 NNA=NEA(IS)*3 IIS=IS+1 DO 100 IANG=1,NNA INODE=NN(IS) - NNA + IANG ANGLE = (IANG-1 ) * (ANG(11S)-ANG(IS)) / NNA + ANG(IS) IF (IS.EO.1)THEN R2=XY(IS) / SIN(ANGLE*DEG) ELSE IF (IS.E0.2)THEN R2 = XY(IS) / COS(ANGLE*DEG) ELSE IF (IS.EQ.3)THEN R2 = XY(IS) / SIN(ANGLE*DEG) ELSE IF (IS.E0.4)THEN R2 = XY(IS) / COS(ANGLE*DEG) END IF IF (NSS.NE.4) ANGLE = ANGLE-45.00 CALL PSPACE (NL.RL.AL,INODE,ANGLE,R1,R2,IOUT) 100 CONTINUE 1 12 MESHGEN Listing 110 CONTINUE IF (NSS.E0.1 .OR. NSS.E0.2) THEN INODE = NN(NSS) + 1 IF (NSS.EO. 1 ) THEN R2=XY(1) / SIN(ANG(2)*DEG) ELSE IF (NSS.EQ.2)THEN R2=XY(2) / C0S(ANG(3)*DEG) END IF ANGLE = ANG(NSS+1) - 45.00 CALL PSPACE (NL.RL.AL,INODE,ANGLE,R1,R2,IOUT) END IF CALL LOAD (NL.RL,AL,DEG.IOUT) 1000 FORMAT (2I5.F12.5) 1010 FORMAT (15) 1020 FORMAT (2F15.9) 2000 FORMAT ( ' Title',/, 2 14, ' , 5. 0, ' 13 ' . 3, ' ,/. 3 ' 1, 2, O, 0, O,Restart',/, 4 ' 0. 0.15. 0,',/, 5 ' 0.0, 0.0, 0.00001 , ,300000. , ',/, 6 ' , , 3, 1, 0,',//) STOP END SUBROUTINE PSPACE it********:) SUBROUTINE PSPACE (NL,RL,AL,INODE,ANGLE,RO,R2.IOUT) DIMENSION NL( 1 ) ,RL( 1 ) ,AL( 1 ) LOGICAL FLAG COMMON NEA(4),NN(4 ) ,NSS,NER,NNA,NNR,FLAG POWER=1.O/NNR R1=R0 J=INODE C0NST=(R2/R1)**POWER IF (NSS.EQ.1 .AND. J.E0.1) THEN WRITE (IOUT,2540) INODE,R1,ANGLE .AND. J.EO . (NN(NSS) + 1 ) ) THEN INODE,R1.ANGLE .AND. J.EO. 1) THEN INODE,R1,ANGLE .AND. J.EO.(NN(NSS)+1)) THEN INODE,R1,ANGLE INODE,R1,ANGLE THEN + NN(NSS) ELSE IF (NSS.EO.1 WRITE (IOUT,2560) ELSE IF (NSS.EO.2 WRITE (IOUT,2540) ELSE IF (NSS.EO. 1 WRITE (IOUT,2540) ELSE WRITE (IOUT,2500) END IF DO 200 1=1,NNR IF (NSS.EO.4) INODE=INODE ELSE INODE=INODE + NN(NSS) + 1 END IF R1=R1*C0NST IF (I.EO.NNR) THEN NL(J)=INODE RL(J)=R 1 AL(J)=ANGLE END IF C C SIMPLY SUPPORTED BOUNDARY C IF (FLAG) THEN C C ONE QUARTER PLATE C IF (NSS.EQ. 1 .AND. J.EQ. 1 . WRITE (IOUT,2540) INODE ELSE IF (NSS.EO.1 .AND. J WRITE (IOUT,2550) INODE ELSE IF (NSS.EO.1 .AND. J WRITE (IOUT,2560) INODE ELSE IF (NSS.EQ.1 .AND. d WRITE (IOUT,2570) INODE .AND. I .NE.NNR ) THEN R1.ANGLE E0.1 .AND. I.EO.NNR) THEN R1.ANGLE EO.(NN(NSS)+1) .AND R1.ANGLE EO.(NN(NSS ) + 1 ) .AND R1.ANGLE I.NE.NNR) THEN I.EO.NNR) THEN 113 MESHGEN Listing ELSE IF (NSS.EQ.1 .AND. I.EQ.NNR) THEN WRITE (I0UT.2530) INODE,R1,ANGLE ONE HALF PLATE ELSE IF WRITE ELSE IF WRITE ELSE IF WRITE ELSE IF WRITE ELSE IF WRITE FULL PLATE I .NE.NNR) THEN (NSS.E0.2 .AND. J.EQ.1 .AND. (IOUT.2540) INODE,R1.ANGLE (NSS.EO.2 .AND. J.EQ.1 .AND. I.EQ.NNR) THEN (I0UT.255O) INODE.R1.ANGLE (NSS.EQ.2 .AND. J.EQ.(NN(NSS)+1) .AND. I.NE.NNR) THEN (IOUT.2540) INODE,R1.ANGLE (NSS.EQ.2 .AND. J.EQ.(NN(NSS )+1 ) .AND. I.EQ.NNR) THEN (I0UT.258O) INODE.R1.ANGLE (NSS.EO.2 .AND. I.EQ.NNR) THEN (IOUT.2530) INODE.R1.ANGLE ELSE IF (NSS.EQ.4 .AND. IF (J.EQ.(NN(2)+1)) WRITE (I0UT.251O) ELSE IF( J.GT.(NN(2)+1) WRITE (I0UT.252O) INODE ELSE WRITE (I0UT.253O) INODE END IF I.EQ.NNR) THEN THEN INODE,R1.ANGLE .AND. J.LE. R1,ANGLE R1 ,ANGLE (NN(3)+1) ) THEN INTERNAL ELSE WRITE (IOUT.2500) INODE,R1,ANGLE END IF FIXED BOUNDARY ELSE ONE QUARTER PLATE IF (NSS.EQ.1 .AND. J.EQ.1 .AND. I.NE.NNR) THEN WRITE (I0UT,264O) INODE,R1,ANGLE ELSE IF (NSS.EQ.1 .AND. J.EQ WRITE (I0UT.2650) INODE,R1 ELSE IF (NSS.EQ.1 .AND. J.EQ WRITE (I0UT.26G0) INODE,R1 ELSE IF (NSS.EQ.1 .AND. J.EQ WRITE (I0UT.267O) INODE,R1 ELSE IF (NSS.EQ.1 .AND. I.EQ.NNR) THEN WRITE (I0UT.2630) INODE,R1,ANGLE 1 .AND. I.EQ.NNR) THEN ANGLE (NN(NSS)+1) .AND. I.NE.NNR) THEN ANGLE (NN(NSS)+1) .AND. I.EQ.NNR) THEN ANGLE ONE HALF PLATE ELSE IF WRITE ELSE IF WRITE ELSE IF WRITE ELSE IF WRITE ELSE IF WRITE FULL PLATE (NSS.EQ.2 .AND. J. (I0UT.264O) INODE. (NSS.EQ.2 .AND. J. (I0UT.265O) INODE, (NSS.EQ.2 .AND. J. (I0UT.2640) INODE. (NSS.EQ.2 .AND. J. (I0UT.2G8O) INODE, (NSS.EQ.2 .AND. I. (I0UT.2G3O) INODE, EQ.1 .AND. I.NE.NNR) THEN R1.ANGLE EQ.1 .AND. I.EQ.NNR) THEN R1.ANGLE EQ.(NN(NSS)+1) .AND. I.NE.NNR) THEN R1.ANGLE EQ. (NNf NSS) + 1) .AND. I.EQ.NNR) THEN R1.ANGLE EQ.NNR) THEN R1 .ANGLE ELSE IF (NSS.EQ.4 .AND. I.EQ.NNR) THEN IF (J.EQ.(NN(2)+1)) THEN WRITE (IOUT.2610) INODE,R1,ANGLE ELSE IF( J.GT.(NN(2)+1) .AND. J.LE. WRITE (IOUT.2620) INODE,R1,ANGLE ELSE WRITE (I0UT.263O) INODE,R1,ANGLE END IF (NN(3)+1) ) THEN 114 MESHGEN Listing c C INTERNAL C ELSE WRITE (IOUT.2500) INODE,R1,ANGLE END IF END IF 200 CONTINUE C c SIMPLY SUPPRTED BOUNDARY L 2500 FORMAT (14. ',0,0,0,0,0, 1, 0 00, , 2( F8 3 , . ' ). ' 0,1 ,2 2510 FORMAT (14, ',1.1,1,0.0, 1, 0 00, ,2(F8 3. . ' ), ' 0.1 , 2 2520 FORMAT (14. ',0.1,1.0.0, 1, 0 00, . 2( F8 3, .').' 0,1 , 2 2530 FORMAT (14, ',0,0.1,0.0. 1. 0 00. , 2( F8 3 , .'),' 0,1 , 2 2540 FORMAT (14. ',0,1.0,1,0, 1. 0 00, ,2(F8 3. . ' ). ' 0.1 , 2 2550 FORMAT (14, ',0,1,1,1,0, 1, 0 oo-. ,2(F8 3. . ' ). ' 0.1 . 2 2560 FORMAT (14, ',1,0,0,0,1. 1. 0 00, , 2( F8 3, , ' ), ' 0.1 , 2 2570 FORMAT (14, '.1,0,1,0,1. 1. 0 00, ,2(F8 3. .'),' 0,1 . 2 2580 /•* FORMAT (14, ',1,1.1.1,0, 1, 0 00, ,2(F8 3. . ' ),' 0.1 , 2 c CLAMPED BOUNDARY 2610 FORMAT (14, ',1,1,1,1.1. i, 0 00. , 2( F8 3, ,'),' 0,1 , 2 2620 FORMAT (14. ',0,1,1.1,1, 1, 0 00, , 2(F8 3 . ,').' 0,1 , 2 2630 FORMAT (14, ',0,0.1,1,1. 1 , 0 00, ,2(F8 3, ,'),' 0.1 , 2 2640 FORMAT (14. ',0.1.0.1,0, 1 , 0 00, ,2(F8 3. .'),' 0.1 , 2 2650 FORMAT (14, '.0,1,1.1.1, 1, 0 00. ,2(F8 3, .').' 0.1 , 2 2660 FORMAT (14, ',1,0.0,0,1, 1 , 0 00. ,2(F8 3 , .').' 0,1 , 2 2670 FORMAT (14, ',1,0,1,1,1, 1 , 0 00, ,2(F8 3, . ' ),' 0,1 .2 2680 FORMAT (14, ',1,1,1,1.1. 1 • 0 00.. . 2(F8 3. , ' ) , ' 0,1 , 2 RETURN c END Q ************ S U B R 0 U T I N E L 0 A D ********* * * C SUBROUTINE LOAD (NL,RL,AL.DEG,IOUT) DIMENSION NL(1),RL(1),AL(1),DIS(36),P(30) 1 ,SUM(4),PELE(4),IFIRST(16),I LAST( 16) LOGICAL FLAG COMMON NEA(4),NN(4),NSS.NER,NNA,NNR,FLAG C C PRINT ELEMENT NODE NUMBERING C DO 290 IR=1,NER WRITE (IOUT,2000) NN(NSS)/3 CALL ELENO (IFIRST.I LAST,IR,NN,NSS) WRITE (IOUT.2010) ( I F I RST(I ) ,I = 1 . 16) IF (NEA(NSS).NE.1) 1 WRITE (IOUT,2020) NN(NSS)/3,(I LAST(I),I = 1, 16) 290 CONTINUE C C PRINT NUMBER OF LOAD POINTS C IF (NSS.E0.4) THEN NUMLP = 2 * NN(NSS) + NSS + 1 ELSE NUMLP = 2 * (NN(NSS) + 1) END IF WRITE (IOUT,2030) NUMLP THICK=10.0 C C CALCULATE THE DISTANCE BETWEEN POINTS C DO 300 I=1,NN(NSS) NN5=NN(NSS)+1 11=1+1 IF (II.E0.NN5 .AND. NSS.EO.4) 11=1 X1=RL(I) * COS(AL(I ) *DEG) 115 MESHGEN Listing Y1=RL(I) * SIN(AL(I)*DEG) X2 = RL(II) * COS(AL(II )*DEG) Y2 = RL(II) * SIN(AL(11)*DEG) DIS(I) = SORT( (X2-X1 )*(X2-X1 ) + (Y2-Y 1 ) *(Y2-Y 1 ) ) 300 CONTINUE C C FOR EACH SIDE CALCULATE THE CONSISTENT LOAD VECTOR C DO 330 ISIDE=1 .NSS IDIR = 1 IDIR2 = 2 DMULT =1.0 SUM(ISIDE) = 0.0 DO 305 1 = 1 , 30 P(I)=0.0 305 CONTINUE DO 320 IEL=1,NEA(ISIDE) ID = NN(ISIDE) - NEA(ISI0E)*3 + (IEL-1)*3 + 1 A = DIS(ID) B = DIS(ID+1)+A C = DIS(ID+2)+B CALL CONSTL (PELE.A.B.C) DO 31.0 d=1 .4 NP = (I EL-1)*3 + J P(NP) = P(NP) + PELE(J)*THICK 310 CONTINUE 320 CONTINUE C C SET THE CORRECT SIGN AND OIRECTION FOR EACH SIDE, C THEN CHECK THE SUM OF THE LOAD VECTOR. C IF (NSS.E0.4) THEN IF (ISIDE.EO.2 .OR. ISIDE.EO.4) IDIR=2 IF (ISIDE.EQ.2 .OR. ISIDE.EO.3) DMULT = -1 ELSE DMULT=1.0/SORT(2.0) END IF C C PRINT THE NODAL LOADS C DO 325 1=1,(NEA(ISIDE)*3+1) P(I) = P(I)*DMULT SUM(ISIDE) = SUM(ISIDE)+P(I) K = NN(ISIDE)-NEA(ISIDE)*3 + I IF (K.E0.NN5 .AND. NSS.E0.4) K=1 WRITE (I0UT.2040) NL(K).IDIR.P(I) IF (NSS.EO.1) THEN P(I) = -1.0 * P(I) WRITE (I0UT.2040) NL(K),IDIR2,P(I) END IF 325 CONTINUE 330 CONTINUE C C PRINT LATERAL LOADS C C DO 340 1=1,NN(NSS),3 C IF (NSS.NE.4) THEN C F (I.EO.1) THEN C WRITE (IOUT.2050) I,1+1,1+2 C ELSE C WRITE (I0UT.2O6O) 1,1+1.1+2 C END IF C ELSE C WRITE (I0UT.2O6O) 1,1+1,1+2 C END IF C 340 CONTINUE C IF (NSS.NE.4) WRITE (IOUT.2070) NN(NSS)+1 WRITE (IOUT.2080) WRITE (IOUT.2090) (I.SUM(I).I=1,4) 1 16 MESHGEN Listing 2000 FORMAT ( 1' 7, '.12, ' ,3,0,0,0, 16, ,4,4,5,0,1,1,1,1,',/, 2' 1,7.70E-05.0.0.',/, 3' 200000.0.0.3.1.2.300.0,0.0,') 2010 FORMAT ( 1' 1 , 16,334, 1,0,0,0,10.0,',/, 1614) 2020 FORMAT ( 113, ',16,334,1,3,0,0,10.0,',/,1614) 2030 FORMAT ( 114, ',1,3, ',/, 2' 1,3',/, 3' 0.0, 0.0',/, 4' 1.0. 1.0.'./. 5' 2.0, 2.0,') 2040 FORMAT (14, ',' ,12, ', 1,',F10.4) 2050 FORMAT (14, ', 3, 1, 0.04',/, + 14, ', 3. 1,0.12',/, + 14, ', 3, 1,0.12') 2060 FORMAT (14, ', 3, 1, 0.08',/, +14,'.3,1,0.12',/, + 14, ', 3, 1,0.12') 2070 FORMAT (14, ',3, 1, 0.04') 2080 FORMAT ('1, 1, 1,0,4, ,,1.0,',/, 1 ' 1 ' ) 2090 FORMAT 1 (4(/,' SUM OF THE FORCES FOR SIDE',12,' IS',F15.6.' sq mm' RETURN END C c *************** SUBROUTINE CONSTL ********** C SUBROUTINE CONSTL (PELE,A,B,C) DIMENSION PELE(4) C C CALCULATE THE CONSISTENT LOAD FOR THE CUBIC SHAPE FUNCTION C C . . • • C 0 A B C C PELE( 1 ) = 0.7125*A - 0.3000*B + 0.0875*C PELE(2) = 1.0125*B - 0.3000*C PELE(3) =-1.0125*A + 0.7125*C PELE(4) = 0.3000*A - 0.7125*B + 0.5000*C RETURN END C Q ************** 5 (j g R • u j i N E E L E N 0 ************** C SUBROUTINE ELENO ( IFIRST,ILAST,IR,NN,NSS ) DIMENSION IFIRST(16), ILAST(16), IMULT(16), JMULT(16), 1 IADD(16). JADD(16),NN(4) DATA IMULT / 3,0,0,3,2,1,0,0.1.2,3,3,2,1,1,2 / DATA IADD / 4,4,1,1,4,4,3,2,1,1,2,3,3,3,2,2 / DATA JMULT / 2,-1,0,3,1,0,0,0,1,2,3.3.2,1,1.2 / DATA vJADD / 0, 0,3,3,0,0,1,2,3,3,2,1,1,-1,2,2 / NO=NN(NSS) IF (NSS.NE.4) N0=N0+1 ISTART=NO*(IR-1)*3 JSTART=NO*(IR-1)*3+N0+1 DO 500 1=1,16 IFIRST(I)= IMULT(I)*NO + ISTART + IADD(I) IF (NSS.EO.4) THEN ILAST(I) = JMULT(I)*NO + JSTART - JADD(I) ELSE ILAST(I) = IFIRST(I)+NN(NSS)-3 END IF 500 CONTINUE RETURN END 117 APPENDIX E Computer Communications APPENDIX E.l WORDSTAR Output on the MTS Zerox 9700 The following commands will transfer a WORDSTAR file on the IBM PC to MTS and then print the file on the Xerox 9700 lazer printer. Require: WORDSTAR DISKETTE WINDOW DISKETTE obtainable from the UBC book store I or G account on the UBC MTS system In WORDSTAR print the file to a disk file. P "filename" y "fileprint" RETURN RETURN Y RETURN X // exit WORDSTAR when printing is finished // // change to WINDOW diskette // A>WSCLRBIT "fileprint" "fileclear" // get the attention of the smart switch with kermit // A>KERMIT SET BAUD 4800 CONN *6 // MTS on port 6, VAX VMS on port 5 // ctrl]C EXIT // using the same diskette run WINDOW // A>WINDOW G // or I depending on the location MTS account // SIG "ccid" "password" CREATE "filenew" %T PC "fileclear" MTS "-filetemp" ASCII RUN PC:WPPRINT SCARDS="-filetemp" SPRINT="filenew" PAR=WORDSTAR SET PROUTE=CNTR SET DELIVERY=CIVL // or CNTR //) CON *PRINT* PORTRAIT ONESIDE COPY "filenew" *PRINT* SIG %EXIT Other commonly used MTS comands are: DIS *PRINT* // displays print statistics // REL *PRINT* // release print to printer // SYS QUE USER // show que or time of printing // "filetype" file names provided by the user // // enclose comments 118 APPENDIX E.2 Transfer of a VAX-VMS File to the UBC/MTS System The following commands will transfer a file from the Civil Engineering VAX 11/730 to the UBC/MTS system. RETURN "VAXid" // sign on to the VAX // "password" SD "defaultDIR" ALLOC TXAO // dial up UBC (228-1401) on modem // KERMIT SET SPEED 1200 CONN G (I) SIG "ccid" // sign on to the UBC MTS system // "password" RUN *KERMIT RECEIVE "filename" ctrlP // ctrl]C on the IBM PC // SEND "filename" CONN EXIT SIG // sign off the MTS system // EXIT LO / sign off the VAX VMS system // " " enclose user id's, filenames, directories and passwords // // enclose comments 119 

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 10 0
Japan 1 0
Nigeria 1 0
India 1 0
City Views Downloads
Unknown 4 0
Cherry Hill 2 0
Glassboro 2 0
Ashburn 2 0
Mountain View 1 0
Sunnyvale 1 0
Tokyo 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

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-0062574/manifest

Comment

Related Items