.V.. ylt.l._:. II )I IVY-l ‘ . .. Ar. C'QVQCIDO’ "II. I: is ’ 53rd"? 0" I 1“ " ' {0’0 gev-Ffipuafllofiay .3. "(LV (flsf' VII-11! I Sft' 0’§:,I.P;.l’ - I) it}! t!!! I Cl! . eh. fol! 5821157..) 'I'I‘ 2; ‘VJI (o 3.! hi £41 lift. 9 30.. r! g: tn... $4.38} I. x ' Q 5 . vino; z.s¢r.3dril.§u..$a :. fin! II 3!! 1. a .‘v I an I} I ’9 {.1 938.; .rlrx . , ‘ 0|- 1 .!. I O.- u- 5.... . U I. :533 . 2 . . urchfi...n.flo.nnmnmu at?! p x 15....“ ’ - (1”; u.- Iliyfl". 5")" (It? 2 {1. Sun». #5.; {Husvfi vi tirfff’llnflfiul!‘ ”WHY. $1135?“ tilt. it.“ 00.; 1.41- . «TI-‘9!- ta;.vld’vlb up; 7.01 19?). v .5?» I 1:) 14......in 71f. . . pa SI! ( : ofi‘lhfaihnvflfg‘v (’1 Infill?!“ Alt. :1... .. '00; to . 4a valueyxolll‘flvv-l In . 1.2.6.510 .34.}:trlpovla! 0901 10.6.93 51...: .z: A 5.31.101.- .ollltt‘)- {rag-r 11' 1: r— '11! I I‘irlllo . V 3.1.... . ... y .A u: .u- .\4.u, A ‘ H 96% MICHIGAN 1» WWW/11mm. This is to certify that the dissertation entitled Correlation of Stress, Strain, Nbrphological Shifts and Permeation-Rate Changes in Polymer Films presented by Scott Allan Nbrris has been accepted towards fulfillment of the requirements for Ph.D. degree in Agricultural Engineering j rofess Date W MS U is an Affirmative Action/Equal Opportunity Institution - 0- 12771 k ""W LIBRARY Michigan State 1 University ~ J *— PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. E DATE DUE DATE DUE DATE DUE mm 2 email {'1 e 'UDU; T j "T57 MSU I: An Affirmative ActioNEqueI Opportunity Institution CWMH.‘ Correlation of Stress, Strain, Morphological Shifts and Permeation-Rate Changes in Polymer Films. BY Scott Allan Morris A DISSERTATION Submitted to Michigan State University in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy Department of Agricultural Engineering 1992 ABSTRACT CORRELATION OF STRESS, STRAIN, MORPHOLOGICAL SHIFTS AND PERMEATION-RATE CHANGES IN POLYMER FILMS. BY Scott Allan Morris This study constructs a quantitative link between stress, strain, morphology changes and changes in the rate of permeation in a polymer film sample. The purpose of this is to provide a simple method of correlating mechanical input and performance changes in polymeric material using a series of simple tests rather than the repetitive construction and evaluation of prototypes. A step-loading test fixture was devised to apply varying levels of stress to a cruciform film sample. An optimization scheme based on the COMPLEX algorithm was devised to determine the material constants of the polymer used in the material and a finite element model for viscoelastic materials was constructed to resolve the state of strain, thinning and change in free volume occurring in the sample material. Small Angle Light Scattering (SALS), and Scanning Electron Microscopy (SEM) were used to inspect the material for gross changes in morphology as a result of mechanical input. Finally, co2 permeation was tested via the quasi- isostatic method to check the film for changes in permeation rate. The study revealed an increase of up to 20% in the rate of permeation in response to a small (less than 6%) increase in free-volume and an even smaller (less than 1%) amount of thinning in the sample. The correlation of these changes with a specific state of strain in the region tested represents a significant step forward since the stress-strain model may be used to resolve similar occurances in other structures made of the same material. More generally, the study represents the inititation of a modelling methodology which may be used to predict overall permeation changes in a particular structure before beginning prototype construction. Copyright by SCOTT ALLAN MORRIS 1992 Dedication This work is dedicated to Dr. Allan J. Morris Frances M. Stearns-Morris Pamela F. Morris for their pertinacious support. ACKNOWLEDGEMENTS I would like to acknowledge the persistent and crucial support of Dr. Bruce Harte in the construction, funding and realization of this project. Thanks also to my major professor, Dr. Larry Segerlind, and the other committee members; Dr. John Gerrish and Dr. Gary Cloud, all of whom have made valuable and knowledgeable contributions to this study. Other major contributors to this study are the A. H. Case Center for Computer-Aided Engineering and Manufacturing, and particularly Greg Fell and Fred Hall for their persistent good humor in the face of endless simple questions. At the Center for Electron Optics, Michigan State University, Stan Flegel and Peg Hogan are to be thanked for their help with the electron microscopy portion of the study. Dr. Dashin Liu, Department of Metallurgy, Mechanics, and Material Science, deserves thanks for his assistance with the vagaries of viscoelastic solid mechanics. Dr. Jack Giacin of the MSU School of Packaging donated crucial advice and equipment for the permeation analysis. Beyond that, the entire staff, faculty and graduate student bodies of both the School of Packaging and the Department of Agricultural Engineering are to be thanked for providing a great place to "push the envelope" and for providing a tough act to follow. Also I would like to thank Edna for fielding many penetrating insights. vi TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES INTRODUCTION OBJECTIVES LITERATURE REVIEW 3.1 Mechanics and Viscoelasticity 3.1.1 Effects of Loading 3.1.2 Material response to a Single Step Function 3.1.3 Numerical Solutions of the Stress- Strain Equations 3.2 Optimization and Parameter Estimation 3.2.1 The Complex Method of Function Minimization Gas Permeation in Deformed Polymer Films Material Analysis 3.4.1 Scanning Electron Microscopy 3.4.2 Small Angle Light Scattering, Spherulite Size and Orientation (AU Db) DERIVATION OF THE FINITE-ELEMENT FORMULATION 4.1 Derivation of the Material Constitutive Equations 4.2 Implementation of the Finite-Element viscoelastic Solid Modelling Program MATERIALS AND METHODS Derivation of Material Properties Constants Permeation Testing Small Angle Light Scattering Analysis Scanning Electron Microscopic Analysis U'IUTU'IUI uhbJNH vii Page xi 10 ll 14 15 17 20 20 21 25 25 37 39 39 45 50 50 6.0 RESULTS AND DISCUSSION 6.1 6 2 Parameter Estimation Method Evaluation Material Parameter Estimation 6.2.1 Error in Material Models and Estimated Parameters Due to Loading History Deviations From the Assumed Heaviside Function Strain Within the Specimens Permeation Changes in the Material Small Angle Light Scattering Analysis Scanning Electron Microscope Evaluation 7.0 SUMMARY AND CONCLUSIONS 7.1 7.2 Summary Conclusions 8.0 RECOMMENDATIONS APPENDIX A. PROGRAM LISTINGS APPENDIX B. MATERIAL DEFLECTION DATA APPENDIX C. PERMEATION DATA BIBLIOGRAPHY viii 51 51 57 6O 66 75 78 81 84 84 85 86 88 174 179 183 Table Table Table Table Table Table Table Table Table Table Table Table 10. 11. 12. LIST OF TABLES Parameter Estimation Evaluation Results. Material Parameter Estimation Results. VISC01.FOR VISC02.FOR P21.FOR P21a.FOR DUMERGE.FOR Documentation for Programs. Change in Position of Indeces at 100% Loading. Change in Position of Indeces at 82% Loading. Change in Position of Indeces at 47% Loading. Summary of Permeation Data at Various Loading Levels. ix 54 58 90 117 148 167 168 170 176 177 178 179 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 1. 2. 3. 8. 10. 11. 12. 13. 14. 15. 16. 20. 21. 22. LIST OF FIGURES Chart of Constitutive Equations (Flugge, 1968). Response of a Three Parameter Solid to an Applied Step Loading. Schematic of Small-Angle Light Scattering Apparatus. The Experimental Arrangement for Photographic Light Scattering From Films (Stein and Rhodes, 1960). Tensile Loading Fixture Diagram Photograph of Biaxial Tensile Loading Fixture. Diagram of Film Sample and Finite Element Grid Used for Material Parameter Estimation. Permeation Cell Schematic. Permeation Cell Used in the Study. Permeation Cell With Film Sample Under Tension. Finite Element Grid Used to Evaluate the Parameter Estimation Method. Strain Gauge Load Cell. Uniaxial Test Fixture With Load Cell Installed. Load vs. Time in Tensile Testing Fixture. General Sequence of Programs Used in Parameter Estimation and Material Response Calculations. Finite Element Grid Used in Calculating Thinning and Changes in Free Volume. Thinning Calculations for Film at 100% Loading. Thinning Calculations for Film at 82% Loading. Thinning Calculations for Film at 47% Loading. Calculated Change in Free Volume at 100% Loading. Calculated Change in Free Volume at 82% Loading. Calculated Change in Free Volume at 47% Loading. 12 22 23 4o 41 42 46 47 48 52 61 62 64 67 68 69 7O 71 72 73 74 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. Graph of Permeation Rate Versus Applied Load. SALS Hv Pattern Before Applied Load. SALS My Pattern After Applied Load. SALS H,Pettern Before Applied Load. SALS Hrpettern After Applied Load. Photomicrograph of Film Before Deformation. Photomicrograph of Film After Deformation. Finite Element Grid Used in Parameter Estimation Evaluation Model. Change of CO2 Concentration in Permeation Cell Versus Time at 100% Loading. Change of CO2 Concentration in Permeation Cell Versus Time at 100% Loading. Change of CO2 Concentration in Permeation Cell Versus Time at 100% Loading. xi 76 79 79 80 80 82 83 175 180 181 182 1 . 0 INTRODUCTION The use of barrier polymers in packages designed for long shelf-life products is increasing. With this increase comes the need for a quantitative assessment of the effects of stress and strain imposed during package manufacture, filling, and distribution on the gas-transport properties of a material. Mechanically induced changes in barrier properties can result in the reduced shelf-life of a packaged product, or the reduced efficacy of a pharmaceutical product even if there is no overt sign of deterioration. An improved understanding of the relationship between mechanical, morphological, and gas-barrier properties of materials will allow producers and users to optimize their use of these materials by obtaining desired barrier properties with lower volumes of more carefully engineered materials. The relationships between loading, deformation, structural orientation, and changes in gas transport properties (diffusion and permeation) in polymers are poorly correlated. Some studies have characterized microstructural changes and changes in gas transport properties of highly crystalline polymers (Polycarbonate, PVC, and polyimide) under simple uniaxial tensile strain. These results often 2 conflict with existing theories of free-volume reduction, however, and do not apply to cases of multiaxial strain, nor to the semicrystalline polymers prevalent in the packaging industry such as PE and PP [O'Brian et al.,1987; Smith and Adam, 1981; El-Hibri and Paul 1985]. Studies have been conducted to characterize changes in the crystalline/amorphous structure of deformed semi- crystalline polymers in terms of an accurately measured state of strain. These studies did not consider the gas-transfer properties of the polymers. [Segula and Rietsch, 1985; Kendra et al., 1985; Burke and Weatherly, 1987; Pan et al., 1987; Schultz,1984]. Other investigators have attempted to correlate strain and gas transport properties, but they have only considered the strain field imposed on the samples in terms of simple uniaxial elasticity/plasticity relationships. They did not consider viscoelastic (time-dependent) relaxation effects, nor did they make an effort to characterize the resultant changes in morphology [Yasuda and Peterlin, 1974; Paulos and Thomas, 1980; El-Hibri and Paul, 1986; Morris and Lee, 1987]. The relationship between stress, strain, changes in morphology, and changes in gas transport properties is in need of a comprehensive, multifaceted study in order to approach mechanically caused barrier property change problems with the proper analytical tools. The significance of the quantitative understanding of this relationship reaches 3 beyond food and pharmaceutical packaging into the fields of: -Separation and purification operations (optimization of membrane structure during the fabrication of separator/filter cartridge structures. -Biomedical processes (such as dialysis and oxygenation via the strain modification of the requisite semipermeable membrane structures). -Extended storage of whole blood/ blood platelets (modification of the polymeric container for optimum gas transport) (NASA MSG-21157). 2 . 0 OBJECTIVES This study was designed to determine the correlation between mechanical stresses, viscoelastic response, shifts in polymer morphology, and changes in gas-transport properties in the semicrystalline polymers used in packaging. A quantitative link between mechanical stress and shifts in the gas-barrier properties in these materials was to have been defined. Specifically, the objectives of the study are: 1. To experimentally determine the viscoelastiijroperties of a material sample and apply the coefficients to a numerical model of the biaxially strained sample in order to predict the stress, strain and time relationships of the material. 2. To characterize changes, if any, in the morphology of the polymer. 3 To measure the gas-transport rate of the material before and after deformation. 4. To link the information gathered from the previous objectives into a comprehensive picture of the functional changes occurring in the polymer film as a result of two-dimensional mechanical stress. 3 . 0 LITERATURE REVIEW 3.1 Mechanics and Viscoelasticity The primary quantities which are considered in the mechanics of materials are stresses, strains, and material properties. Stresses result from forces acting within the body being considered. Strains are the result of differential movements within the body, and are usually described in terms of percentage of a referential measure (eg. final length divided by original length in the linear case). Displacements are the movement of referential point(s) within the body relative to some external coordinate system and are related to strain. The field of the analytical mechanics of deformable bodies is primarily concerned with the relationship of these three conceptual quantities, and constructs three broad classifications of families of equations to relate the quantities to one another: equilibrium conditions (or equations of motion for dynamic problems), kinematic relations, and constitutive equations. Equilibrium conditions are material independent, and define the stresses (or conditions of motion) acting on the body in question. Kinematic relations, which are also material-independent 5 6 relate the strains within a body to its displacements. The constitutive equations are material dependent and relate the stress and strain within the body. Since the materials may vary widely in their response, there are many types of constitutive equations falling into four broad categories; elastic, viscous, plastic, and viscoelastic. Elastic materials may be defined as those materials which store energy without significant loss (as with a stretched metal spring), and stress is linearly proportional to strain. Viscous materials, by contrast, dissipate energy without significant energy-storage capacity (demonstrated by pouring water from one glass to another) and exhibit stress in proportion to strain rate. Plastic behavior combines the preceding two concepts in the following manner: A plastic material will store energy in an elastic fashion until the "yield point" is reached, after which the material begins to irreversibly deform without further energy storage. Plastic behavior usually does not account for the rate of strain on the material (Flugge, 1967). Viscoelasticity is concerned with the study of the materials which have a stress-strain relationship with characteristics that do not exactly fit the concepts of plasticity, viscosity or elasticity. Viscoelastic materials can both store and dissipate energy in a manner where the rate of strain is dependent on the time-history applied stress. Viscoelastic constitutive equations are often 7 composed of both viscous and elastic equation elements (hence the term "visco-elasticity") , and the constitutive equations of these types of materials may be linear, with constant material coefficients, or non-linear (where the material coefficients contain terms which are a function of stress, strain, or their'derivatives) (Ferry, 1980). In this study, for the sake of simplicity, the response of the material was assumed to be linearly viscoelastic. For most constitutive equations arising from the study of deformable bodies, the two basic elements of viscoelastic models are the elastic spring element where Stress = Constant ° Strain (1) 0= I?'8 (The constant in this equation is usually the elastic [Young’s] Modulus, E). The viscous dashpot element used is Stress Constant ° Strain Rate (2) a constant ° t In fluid mechanics, the constant is often taken to be u, the viscosity coefficient relating shear stress and shear rate. The more complex material responses are often modelled by combinations of simple elements, some of which are reproduced in Figure 1. The three-parameter solid is the constitutive model that therential equation Model Name Inequalities o———vvw——O elastic solid 0 - 4.4 O—E—__—o viscous fluid 0 "' 9d o—ww—(E—O Maxwell fluid a *7)" - 9.. 0 Kelvin solid 0 a 9" + q“ o +p,a‘r - q.¢ + gt 3-paramcter solid 91 > [’39. Figure 1. Chart of Constitutive Equations. (Flugge, 1968) 9 will be used in this study since its compliance curve most nearly matches that of the types of polymers under consideration (Ferry, 1980), and it encompasses all of the preceding models in the table. 3.1.1 Effects of Loading The type and degree of loading is critical in predicting the response of viscoelastic materials. Since the analytical solution of the constitutive differential equations is often accomplished using Laplace transforms loading regimes considered are usually those which can be constructed of the more common functions considered within the context of Laplace transforms: Dirac delta (spike) functions, Heaviside unit (stepped) functions and ”ramp" functions. Superposition and the time-shift principles (the so-called t- and s-shifts) may be used to combine these functions to mimic many "real-world" situations (Speigel, 1965). In this study, the loading will be limited to a single step function. 10 3.1.2 Material Response to a Single Step Function. The response of a three-parameter solid may be shown by solving the model's governing equation, o+plo = qoc +q1¢ (3) using the Laplace Transform, this operation gives 1 0 Qo’qia or — l 0:2(i11- where .q_° = A (5) 0'; sum) g1 where o. is the initial step loading value. When rearranged and operated on by the inverse transform, the equation becomes e=%:[%(1-e'“) + pledt] (6) 11 Note that =T=Lfl e(t 0) ca (7) g ( t=eo) 3.5% which is illustrated in Figure(2). 3.1.3 Numerical Solutions of the Stress-Strain Equations In all but the most simple geometrical configurations, the exact analysis of the state of stress and strain at a chosen point within the body to be considered verges on the impossible. Sections with irregular geometries, sharp corners, inclusions such as holes or slots, or boundary conditions that are less than ideal all render the analytical formulations of the states of stress and strain unsolvable in any practical sense (Cook a Young, 1985). Numerical methods are viewed as a means to closely approximate the solutions of the partial differential equations of stress and strain of the body to be analyzed. The finite difference method replaces the partial differential equation (PDE) and boundary condition terms with a system of equivalent difference equations. The solution of this system of simultaneous equations then yields the solution to the PDE at each node. This method is useful in W P a L G i GfiGz 2 o——J\N\/\———— + -- ELI- G 43 o qo 1+6; 1 g G n “1 c c. to o u: For a stepped load 00 applied at time t=0 Figure 2. Response of a Three-Parameter Solid to an Applied Step Loading. 13 that the formulation is relatively simple. Unfortunately, the rectangular grid geometry for the points used in the analysis is simple as well--and inflexible (Ugural, 1981). Although it is possible to produce a mesh of triangular elements, or re-map them into circular or spherical coordinates, extremely irregular boundaries or unusual shapes will cause the model to match the spatial coordinates of the actual object very poorly (Paulsen, 1992). The finite-element method, although burdened with a more difficult formulation, has the advantage of not being limited to a strictly rectangular or triangular approximation of the body's geometry. Rectangular and triangular elements may be mixed, curved elements may be produced (or closely approximated), and the dimensions of the elements may be allowed to vary in order to more closely approximate the geometry of the subject. In the system used to model the viscoelastic response of the materials in this study, the finite element formulation also will be shown to conveniently accommodate the time-dependency of the material as well as the irregular boundary conditions of the test samples (Segerlind, 1984, Boresi and Sidebottom, 1984). 14 3.2 Optimization and Parameter Estimation Determining the material coefficients for the defining differential equation of the three-parameter model is a troublesome aspect of the mechanics of viscoelastic materials for all but the simplest types of tests and sample geometries. There are several standard ASTM tests for these types of materials, but the results are quite simplified and are of limited use for the modelling of complex systems. A great deal of effort has been devoted to producing a good set of material coefficients for many types of constitutive models, both linear and nonlinear (Augl and Land, 1985) (Perry, 1980) . Most methods require simplified test protocols or sample geometries and may not return values that are useful in large-scale modelling of the two-dimensional films used in this study. Investigators in the system optimization field have occasionally used finite element modelling to describe the objective function of the optimization scheme at hand (Jehle & Mlejveck, 1990) . The usual scheme in such optimization models is to computationally fit the material coefficients of the model so that the difference between the output of the model and a set of observed data, or desired outcomes, is minimized. This scheme, although computationally intensive, is useful in the optimization of a system that is too 15 complex to describe analytically, or does not possess the simple geometries so often used in standard measurement techniques. This method lends itself well to the problem at hand since the already developed finite element model can be used as the objective function. The data taken from optical measurements in existing test fixtures may be used as the data to be matched. The method developed for this study was that of minimizing the difference between the finite element model's calculated values for strain and the experimentally observed strain data over a series of equally spaced time-steps. The advantage of this method, particularly in a biaxially stressed system, was that the material parameters can be retrofitted to the observed. data. without the need for additional testing equipment. Additionally, this method shows great promise for the analysis of material properties in other situations where the material is in some unusual configuration, is in use and may not be removed for analysis, or' may only' be studied 'via some non-contact methodology. 3.3.1 The Complex Method of Function Minimization The Complex method, an extension of the Simplex method, is an efficient and versatile tool for finding global optima within the domain of possible solutions (Box, 1969). It 16 avoids some of the convergence problems which occur with simplex methods and may be used with almost any number of variables or type of objective function. It also has the notable advantage of being able to accommodate restrictions on both the limits of the estimated parameters and the regions to be considered feasible. It has the further advantage of not requiring the calculation of derivatives. In operation, the Complex method randomly seeds'a‘number of vertex points about the allowable domain of the function to be minimized (creating a complex polygon) and calculates the value of the objective functions at each of these points. The difference between these points is considered, and if the difference exceeds a preset parameter value then the vertex returning the lowest value is replaced by another point located a distance along a line located through the rejected vertex and the centroid of the complex polygon. 'This process repeats until the difference between the vertex points is less than the minimum (B). At this point, the centroid of the isoplethic polygon is considered to be the minimum (ringed, in a sense, by the vertices of the complex polygon), and the computation stops, returning the values of the parameters at the centroid, and the centroid objective function value. 17 3.3 Gas Permeation in Deformed Polymer Films The rate of gas permeation through materials can be correlated to several phenomena which occur during the deformation process. The simplest is the result of changes in section thickness which occur as a material is strained (half the thickness yields roughly twice the permeation). In most polymers, however the processes which occur are much more complex, and the changes are not nearly as simple as with simple thinning phenomenon. Under large strains, polymer films deform visco- elastically, changing their morphology from a mass of nearly random macromolecular chains to a well oriented series of parallel fibrils. Experimental research has shown that the fibrillar component of these oriented materials has a permeation rate that is orders of magnitude lower than that of the amorphous material surrounding it (Williams and Peterlin, 1971). From this, one may correctly assume that extensively orienting a film and creating a high fibril content. will cause the ‘permeation rate of 'the film to decrease drastically. This phenomenon is commonly used in industry to tailor the permeation rates of polymer films to a specific application, and to reduce the amount of raw material needed to perform the barrier function of a particular package. 18 Currently, the most visible example of this method is in the polyethylene terpthalate (PET) bottles, used by the soft- drink industry, which are oriented both axially and radially in order to retain carbonation for the required period of time. Careful manipulation of both the molding and orientation process during production has resulted in a substantial material reduction since the introduction of the bottle in the 1980’s. Under lesser strains, the material behaves much differently than in the large-strain orientation described above. Under' small strains, polymer free—volume' and diffusion increase (provided that Poisson’s ratio is less than 0.5, which is approximately the case for most organic polymers other than rubber) until a plateau is reached at relatively low levels of strain (0.05 < 8 < 0.20 for polyethylene). The increase in free volume is accounted for in the small-strain equation for thinning 82:4: (ox+oy) (8) =-p(ex+ey) (Gere and Timoshenko, 1984). 'The "Free Volume" of a polymer may be thought of as the volume of the void-free solid versus the space occupied by the actual molecular chains comprising it” .A long chain polymer may be folded and bundled such that there is sufficient space between the folds to allow the slow passage of sufficiently small gas molecules. 19 For a square area of material x= y p=o.2s ‘9’ 020.1 y: e,--0.25(o.1+o.1)=-o.os ’- New Area (x+e,,) -(y+ey) (x+0.1x) '(y+0.1y) (1°) 1 . 21 (Original Area) New Thickness (241,) =(z-0.Sz) =(o.9Sz) = 0.95 (Original Thickness) New Volume=(1.21 (Original Area) - 0.95 (Original Thickness) =1 . 14 Original Volume (11) Note that for large values of the elastic modulus, E, the effect is more pronounced to a limiting value equal to the ratio of new surface area to original surface area. Additionally E may be replaced with E(t) by superposition to give a small-strain linear viscoelastic formulation (Flugge, 1968) with similar results. After the ”plateau" of maximum free-volume expansion has been reached, the aforementioned effects of orientation begin to dominate the system and permeation will drop off to well 20 below its unstrained value. Sorption steadily increases as well and appears to approach an asymptote as the levels of strain increase (Yasuda et al, 1964; Yasuda and Peterlin, 1974). The studies mentioned above consider the viscoelastic polymer in terms of an elastic-plastic model and do not account for the time-dependance of the material response in terms of any overall material constitutive model or system. 3.4 Material Analysis 3.4.1 Scanning Electron Microscopy Scanning electron microscopy (SEM) is an often-used technique for the characterization of the surface structure of polymeric materials (Roulin-Maloney, 1990). The usual procedure utilizes either an as-produced surface or a surface created for the analysis (through cryomicrotomy or similar method) as the specimen over which a replicant surface is cast using a high resolution casting material. Often the surface to be examined is treated with some type of etchant such as p-xylene or chromic acid to increase the resolution of the differing regions within the distorted material (Jang et al, 1985, Hashimoto et al, 1976). Once the surface (or its replica) is prepared, the material is desiccated using critical point extraction to remove the 21 solvents and residual moisture which may contaminate the interior of the microscope. The surface to be examined is electroplated (to dissipate any charge that may be built up from the scanning electron beam), and an image may then be recorded. Recent developments in scanning force electron microscopy (SFEM) and scanning tunnelling electron microscopy (STEM)‘have led to very high.resolution of surface characteristics and have made quantifiable roughness measurements possible (Reiss et al, 1991; Howells et al, 1991). The equipment for these SFEM and STEM, however, was not available for use in this study. Although the surface profile information ws not be as good as with SFEM and STEM methods, the geometry and degree of asymmetry of the spherulites should have been accurately measurable. 3.4.2 Small Angle Light Scattering, Spherulite Size and Orientation. Small angle light scattering (SALS) is commonly used to determine the degree and direction of orientation of many types of materials. In its simplest form, the apparatus is simply a highly collimated (preferably coherent) source of monochromatic light which is passed through a polarizer then through the material to be analyzed and through a second, adjustable polarizer (analyzer) (Figures 3 and 4). In the 22 . monoummmc Emanumom £qu magnfigm no 00828 .m 38E cont—20m .893 .. .........J..... Ezm m .oei/«aéiiiiSsai US$92 580m 23 ‘ Polorizer )“w Photographic Film Figure 4 . 'Ihe Ebrperimental Arrangement for Photographic Light Scattering From Films (Stein and Rhodes, 1960). 24 case of high molecular weight polymers, the spherulite size may be estimated (using a HeNe laser with A=632.8 nm) as a, = 4.03 / 41: sin(0./2) (12) where R0 = average radius of the spherulite and 0. - maximum scattering angle (Stein & Rhodes, 1969). Orientation may be obtained by observing that the drawing of the polymer will pull the polymer chains in both the crystalline and amorphous regions of the polymer into an oriented state which is most often compared to a series of oriented rods. The scattering that occurs from this fibrillar component will produce an bias in the scattering pattern normal to the direction of orientation of the polymer (Rhodes and Stein, 1961; Stein and Rhodes, 1960) which may be used (when corrected for sample thickness) to correlate the degree of orientation at the particular sampling area to the rotation of the biased maximum. The usual reason for attempting to measure these changes in orientation is to observe any fibril formation or re- formation which might be occurring in the polymer during deformation. Since the onset of a decrease of permeation rate is linked to this fibril formation, the determination of the point at which these changes begin to occur will be useful in the construction of an accurate predictive model for permeation changes. 4.0 DERIVATION OF THE FINITE ELEMENT FORMULATION Starting with the equilibrium equation aUJ+F1=o i,j=1,2,3 (13) and the viscoelastic constitutive equation 3 aij'f(G1jn(t‘t))(-bi%:t—))dt (1" 0 where 1010k01-102' 3" Using the general convolution notation, which states that for any functions A' and B', 1: [AN t-t)(6—B%fl)dt . a’w’ (16) 0 the viscoelastic constitiutive equation may be expressed as ainguunu (17) 26 where Gun"; [62“) -Gl(t)]6116k1+%61(t) [811611+61161k] Rearranging the terms of equation (17) gives it be (1:) 011.]‘G11u ( t-t)(—k$l?_)dt C =1f[G(t)-G(t)]0115 (27%)“ wl no be (t) *611!G(t) [511°11‘5115111(_—%r—)dt NI For equation (19a) 611%?ij for j =k 51:51:19: (“PM =3»... 8']. Using the identity equation (19a) becomes (18) (flkfl 09b) (20) (21) hie 27 C 0 6t (23) C = l _ 5(‘11"'322"'¢3:s)) -3-‘[[Gz(t) G1(t)]6u( 5: dt Using the identities 31:1“11: t: symme ry 61136.11 (23’ and producing the identity (51k511+51151k)511:51k5j1511+51151k511 =51k51153k+511515521 =6 +6 F’ 11:53:: 11 11 ‘2” “bquj+511Fij =25}: allows (19b) to be rewritten as 28 combining (22) and (25) gives which may be expressed in convolution notation as Defining and O p 0 ll 01.. 02.. 030: 0‘08 6 t C 1 be (t) 5‘ 1 3&1“) [51:511*°1151t] +dt .2031”) a: at t C a,,-%f[G3(t)-G1(t)]6,,a?" + fe.%1dx 0 0 aij'% [63(t) -Gl(t)]611*ekk+g1(t) .311 .3: 1;}: a”; 3:1 7 ‘3': 1:2! l(£aé+iaz 2 a. a, 011' °xx 0338 a”, “33 on 0128 ox}, (25) (2‘) (27) (23) (29) 29 (27) can be written as 0 =11 are k kl - (3°) k,fl - 1'.2°,3°,4° where '§(c,+2c1) %(Gz-G1) §(G,-Gl) 0 ‘ items) l(c-G) o 2 1 2 1 [AJ- 3 3 (31) §(c,+2cl) o 1 L 5Y1". 361‘ Using the functional I-f[%Gunteuteu-F1tu,] dV-[(Sw,)dA (32, v I (Christensen, 1971) the solution to (32) may be found, 30 provided that the following boundary conditions are satisfied: Ganja?! on B. ‘33) uisAi on Bu where B, is the boundary over which the tractions S, are specified, B. is the boundary where the displacements A. are specified, and n1 describes the boundary unit normal vector (positive outward). When the boundary conditions are satisfied, the first variation, 61, vanishes and the solution may be found by finding the stationary value of I. The functional I may be discretized over the relevant region, as a set of subregions: 0'1 I-E-gf (ekAmte;)dV-£ I (11.55?) dV-E I (u.‘¢S.‘)dA .‘1 0.1 v‘ v' a: (34) k'1°,2',3°,4' I'LZ where e represents the subregions. The displacements may be found via the matrix formulations (Zienkiewicz, 1971) {2'} = [B‘] {U} (35a) {11'} . [N‘JIUi (35b) f: 31 where [N] is the shape function matrix, {u} are the elemental displacements, {U} are the nodal displacements, [B'] is the matrix relating strains to nodal displacements, and {e'} is the strain vector. Substituting the equations (35a) and (35b) into (34), integrating, then performing the convolution gives I=3‘-{Ul"#[10*{U}-{U}’*R (as) where the stiffness matrix [K] and force vector {R} are, respectively, [It] =Elk'1=§1f [a :1 ’[A] [3'] av v' (37) {R} am: 0} =£(f [N‘] TF°dv+f [N‘] 1use} (111) v. a. and [ ]T represents the transpose of a matrix The first variation of (36) gives a result similar in appearance to the elasticity finite element formulation: OI=5{U}7* [K] *{Ui-GIUV" {R} =0 m *{U} -{R} -0 (3” [K] *{U}={R} From the preceding derivation the elasticity finite-element formulation 32 [K] {Displacement} = {Force} [K] {U} ={F} (39) then can be shown to have a corresponding time-dependent formulation, [K] *{m-{R} (40) where (40) represents the integral t I[K(t-t)]d{U(t)}d‘t'-{R(t)} m) t-O Discretization of (41) over n equally spaced time-steps 2 [K(tn-t.)] {AU(t.)}-{R(t,)} (42) Pl 'where (AU(t_)} represents individual displacements from time tn to th' This formulation gives each displacement as a function of all previous displacements, plus a possible initial step displacement although this is often taken to be zero. The other individual components of equation (42) are [K] The stiffness matrix which contains spatial data and the viscoelastic modulii 33 [K] =2: IK°1 .. (m =2: f [la-1'14] [3'] av o-i it Where [B‘] is the shape function for element e and [A] is the matrix containing the time dependent material properties derived from the material constitutive equations. The coefficients in [A] are . .Em an 3:: 1_p3 ‘aiz=azi=%tl§Ill (44) “F _a11(1-u) 33“_—j§——' and {R(t)} is the force vector at time t. 34 4.1 Derivation of Material Constitutive Equations The three-parameter (Maxwell) model for viscoelastic solids, has the defining differential equation o+p16 =q°e+q1é (45) operating on (45) with the Laplace-transformed Heaviside unit step function Sloan“) ] =ao(8) -[aou(t)e'“dt‘ (46) 1 I—a ' 0 :33 produces the transformed differential equation [%+p,]a—o=[qo+qis]'e' (47) From this, and the superposition.principle, one may define a viscoelastic shear modulus (Flugge, 1968) which may be used in place of the elastic shear modulus in the material properties matrix of the finite element formulations 35 1+ zc(s)- ' 10’ tee) €k+th The shear*modulus may be expressed (after rearrangement and taking the inverse transform) as a function of three groups of coefficients: -q.(c) (‘9) =_l. - -22! % G(t) ml (1 ‘1)9 substituting .; 2g. = -29. (50) K; (1 m ) =3! e gives G(t)=K1[l-Kze'x’m] (51) It is worth noting that the value for G(t) when t = w is simply G(oo) =fi- = K1 (5:) This equation for G(t) may be substituted into the general 36 equation for the elastic modulus E = 26(1+u) (53) where Poisson's ratio, u, is constant to give em =(19t1-k,e‘*s“’1 (1+u) (so which are used in the previously described material properties matrix [A] defined in (44). 37 4.2 Implementation of the Finite-Element Viscoelastic Solid Modelling Program To implement the previously derived finite element method, the program VISC02 was written in FORTRAN-88, and implemented on the Michigan State University Case Engineering Center VAX-8650 as well as the CRAY Y-MP4/464 at the National Center for Supercomputing Applications at the University of Illinois, Urbana-Champaign?. The implementation of this method is rather straight forward since many of the components of the algorithm are taken directly from elasticity finite element formulations. The substantial difference is that the residual term.must be calculated for each time-step, and is a function of all of the preceding time-steps. This particular version calculates many of the components of these terms at the beginning of the program, then holds them in memory arrays, so that the program is not constantly recalculating the same values for the residual term components. It should be noted that the large number of residual terms generated to account for the stress] strain ”history” of 2 The code for VISC02 is included as part of Appendix A as is an earlier version, VISCOl, which uses disk memory to store the residual matrices, and is usable (although very slow) on smaller computers. 38 the material demands a great deal of storage capacity. These requirements can be mitigated somewhat by the use of banded storage techniques and other data-compression algorithms, but the final result often remains memory intensive. 5.0 MATERIALS AND METHODS 5.1 Derivation of Material Properties Constants The primary material used in this study was 1 mil cast polyethylene film (PW-242 ”Flex-O-Film”, Flex-O-Glass, Inc. Chicago, IL 60651). Actual thickness of the film was checked using a TMI 549 M micrometer (Testing Machines Inc. Amityville, NY) and was found to be .85 mil (0.00085"; 0.000216cm) over all parts of the film. Material properties were derived using a cruciform sample subjected to a stepped loading in a biaxial tensile fixture (Figures 5 and 6). The samples were marked with a non-interactive acrylic ink (Hunt Mfg. Co., Statesville, NC 28677) in a pattern to conform to the grid design selected for use with the P21.FOR and P21a.FOR software (Figure 7). This pattern was chosen as a compromise between minimizing the number of data points to be collected and conformity to the sample during deformation. Since the practical limit for the finite element software is approximately 15% total strain, this limit was found by trial and error (using linear samples) to be approximately 2800 psi. (19.30 MPa) and was taken to be the limiting factor in the strain level applied to the film. The (39) 40 Test Specimen Supported Clamp r—fi ‘_ r—1 _J'_‘l ' 1 u u 1 If u (1 1 A. \ / Tension Cables / l j Y.uis—’D‘_M‘m Test Fixture Diagram Figure 5 . Tensile Loading Fixture Diagram. 41 Figure 6- Photograph of Biaxial Tensile loading Fixture. 42 ouwcwm coo mHmEcm Edam mo Ecumofio 8 ' m .cowuoEHumm Houofioucm Hangman: you some capo newsman .h ousmflm .. . W\\\\\\\\\\\ 1:, V JilJ "38m ..\ \\\\\\\\\\\\\ « 4 \\\\\\\\\\\\\ 'L\\" 43 samples were then loaded to 47%, 82%, and 100% of maximum strain (1316, 2296 and 2800 psi or 9.08, 15.83 and 19.30 MPa respectively). The choice of loading levels are a result of the standardized size of the dead-weights used to load the fixture's cross-beams. The deflection of the material was recorded against a background grid of 0.2" x 0.2” (0.51 x 0.51 cm) squares using a VHS format camcorder (Panasonic VHS Reporter) and then played back using the still-frame feature of an RCA-500 VCR. Data was recorded manually from the film ‘markings and measurement grid previously described. The VHS format records images at 1/30 sec. intervals which allows the deformation history to be recorded.over a span of 6 frames of tape (1/5 sec.) . Since the start of the test was not synchronized with the frame recording of the camera, the exact start and end-points of the deflection history were taken by extrapolation. Initial tests showed a significant amount of rebound-in the film due to aneunknown factor in the test fixture. By experimenting with the mass distribution on the loading beams of the tensile fixture it*was found that the rebound.could.be nearly eliminated by placing the dead-weights near the ends of the beams. This acted to reduce the amount of beam flexure and rebound during material deflection without affecting the level of load placed on the sample. Once the data points had been recorded, the deflections were calculated and averaged, exploiting the symmetry present 44 in the test sample to minimize the calculation necessary to construct models for the mapping of strain fields. 45 5.2 Permeation Testing Permeation testing was done using the quasi-isostatic method. A permeation cell was devised with extended clamping arms which can accommodate the sample held in the tensile fixture (Figure 8). Carbon dioxide flowed through the lower chamber of the cell at a regulated rate of approximately 50cc/min. The upper chamber of the permeation cell (figures 9 and 10) was ventilated with low pressure compressed air for a minimum of 60 minutes after the cell was clamped on the film specimen to allow the CO2 in the lower chamber to reach 100% concentration. At this point the upper cell was sealed and the headspace of the upper chamber was assayed at approximately seven minute intervals for CO, concentration by analyzing 1 cc. syringe-drawn aliquots of the headspace gas with a Carle 2153-B Gas Chromatograph and Spectra-Physics 2400 Computing Integrator. The concentration values were plotted using the integrator's least-squares3 curve-fitting software, and the rate of increase ascertained in order to calculate the rate 3It should be noted that the least-squares fit is used as a matter of convenience since it was part of the integrator software. The actual curve is more closely defined by an exponential time curve [ YaC,(1-e‘“) ] , but for the purpose of rate determination either will suffice. 46‘ .oflgmcom Sou coflmofiuwm .m. 0.53m .50 5:38.60. com vac—am 38333—8883....83. cementum » 5:”— 47 Figure 9. Permeation Cell Used In the Study. 48 Figure 10. Permeation Cell With Film Sample Under Tension. 49 of permeation. A standard gas sample (5.05 % C02, 21.3% 0,, balance N2: AGA Specialty Gasses, Maumee OH 43537) was used to calibrate the detector response, but due to the availability of only' the single. concentration all calculations of permeation rate are done at that concentration in order to minimize error due to detector non- linearity. The curve-fitting software gave the equation of the line plotting the rise in concentration versus time as Ax2 + Bx +6 8 y (55) where x represents time and y is the concentration at time x. The lesser root of the quadratic equation may be used to construct the point x’ at which the concentration passes through the standard value y’ x = -a-‘/a' 2: Mc-y’) (56) The time-rate of concentration change is taken at this point as Dc/Dt = Dy/Dx - 241 x’ + s (57) and used for permeation calculations. 50 5.3 Small Angle Light Scattering Analysis The Small Angle Light Scattering (SALS) analysis fixture was constructed such that material under tension in the tensile fixture could be analyzed. The Helium-Neon laser source produces a polarized beam of light which was used with an analyzer to record the H|I and H, patterns photographically. The 00 was measured against a grid on the projection screen and the u., or rotation of the flare pattern, was measured using a slit photometer. 5.4 Scanning Electron Microscopic Analysis Surface features of the film in both its unstressed and maximally stressed state were cast in the tensile fixture using Reprosil Type 1 light (L.D. Caulk Co., Milford DE 19963) and then a transfer casting was made using Spurr's epoxy resin (Klomparens et al, 1986). The replica was then sputter coated with gold, and examined in the JEOL-35 Scanning Electron Microscope (Center For Electron Optics, Michigan State University). The samples were checked for change of surface features at the center within the area of the sample covered by the permeation cell where the largest strains were predicted to occur. 6.0 RESULTS AND DISCUSSION 6.1 Parameter-Estimation Method Evaluation Two parameter estimation programs, P21.FOR and p21a.FOR were written to extract material property constants from deflection data using tests with known levels of applied loads. The programs utilize the COMPLEX function minimization methods incorporating the VISC02 finite element program as the objective function, attempting to minimize the difference between the deflection predicted by the finite element software and the "real world" deflection data. In order to evaluate the accuracy and utility of the parameter estimation software P21.FOR and P21A.FOR, the previously-described finite element software used dummy material property constants to generate nodal deflection values as a substitute for observed data over an arbitrary 10 time-step period (figure 11) . The parameter estimation software was then used to try and re-extract the original constant values. Records were kept of the number of function-evaluations required for the model to converge. The effect of varying the operating parameters on the final accuracy of the returned estimate was ascertained 51 52 Center of Syrmetry 1.0E4 2.0E4 1.034 A 7 Thickness 1.0 (l) (3) 2T 45 ji. (2) (4) /////////////////W L 1 J r 0 5.0 10.0 Figure 11. Finite Elanent Grid Used to Braluate the Parameter Estimation Method. 53 (table 1). Different numbers of complex polygon vertices were used in the optimization scheme, and the value of B (the intrapolygonal variation of the vertex values, below which the model is assumed to have converged) was varied as well. 54 hflfle»l. Parameter Estimation Evaluation Results Number of Elements: Number of Nodes: Estimation of: 4 9 K 1 Upper Limit of Estimation: 1.0x10‘ Lower Limit of Estimation: 1.0x102 Starting Value: 5.0x10‘ 3 Vertices Nodal Error 5 at_§entrgid 0.0 0.240% 5.0 0.002 2.5 0.002 1.0 0.002 0.5 0.002 Nodal Error 8 amasentrgid 0.0 0.082% 5.0 0.074 2.5 0.003 1.0 0.008 0.5 0.004 Nodal Error 6 W151 0.0 0.001% 5-0 0.003 2.5 0.003 1.0 0.002 0.5 0.003 Nodal Error 3 at Cegtggig 0.0 0.063% 5.0 0.029 2.5 0.008 1.0 0.003 0.5 0.000 Number of 30 38 38 38 38 Number of 34 38 56 64 75 Number of 56 56 62 89 93 Number of 81 89 99 111 126 55 Table 1 (cont'd). Estimation of: E, Upper Limit of Estimation: 1.0x106 Lower Limit of Estimation: 1.0x102 Starting Value: 5.0x10‘ 3 Vertices Nodal Error 8 10.0 2.214% 5.0 0.319 2.5 0.113 1.0 0.056 0.5 0.051 4 Vertices B Nodal Error 10.0 1.125% 5.0 0.024 2.5 0.042 1.0 0.006 0.5 0.003 5 Vertices Nodal Error 6 at Centroig 10.0 0.107% 5.0 0.189 2.5 0.189 1.0 0.036 0.5 0.005 6 Vertices Nodal Error 6 at Centzgig 10.0 0.060% 5.0 0.021 2.5 0.021 1.0 0.017 0.5 0.001 Number of 54 85 91 107 112 Number of Exaluatiens 99 102 102 154 187 Number of Exaluations 128 138 147 176 201 56 As is shown by the tabulated values, the model returned a very accurate estimate of the data with relatively few evaluations of the finite element objective function and a very coarse model grid. It must be noted that the ”error" value to be minimized is a sum of the errors of all points over all time-steps, and thus is extremely situation dependent. The simple study shown illustrates the relative efficiency of the method even when constraint values are used. This computation-intensive type of parameter estimation becomes sl w'as the number of nodes in the objective function model increases, so the software was implemented on the Cray Y/MP 4-464 supercomputer at the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign. The results of these trials indicated that with the parallel processor environment available the computation time shrunk by several orders of magnitude (from tens of minutes to less than three seconds). The software is extremely efficient at estimating simpler cases such as elastic (time-independent) deformation. Although many simpler methods exist for some elastic problems, this method may find use in evaluating unknown materials of unusual geometry. 57 6.2 Material Parameter Estimation The film used in this study was evaluated at three different load levels (the 100%, 82% and.47% levels.described in section 5.1) using the biaxial tensile fixture. Although the fixture is capable of axially independent loading of a sample, the small amount of material available for the study confined the tests to regimes where both axes' were loaded at the samelevel simultaneouly. Initial tests showed that the entire deflection of the material occurred within approximately one fifth of a second (6 frames of videotape) and.no further deformation was recorded.over'periods of up to a week. As will be subsequently shown, this is due to the load-time history of the material. Data taken from the video record (Appendix B) was normalized and converted into displacement files for use with the parameter estimation software P21.FOR and P21A.FOR. Results of the material parameter-estimation are shown in Table 2. 58 Table 2. mterial Parameter Estimation Results 100% ' 43335 136323 5.55 Error‘ 1 5.63% 81.10% 82% 53742 142661 12.51 Error J. 4.15% 79.18% ‘ . ~ ~ . .,_.l___- —- j- —-7—~— --‘—*———«——————~~r - W ~ ~ ‘ , , , 47% 1 40473 155974 ’ 96.02 Error 3 4.72% 40.41% I It should be noted that this parameter estimation method was designed for use where the time history of the material would be a significant factor in the results of the rest of the test (eg. that the permeation testing would occur before the material had reached an equilibrium level of deflection). In this particular case, all of the strain occurred over such a short time period and other material measurements occurred ‘Note that the errors referred to in table (2) are the average absolute cumulative deviation of all points over all timesteps relevant to that parameter estimation run. The K, parameter is derived from a single timestep as previously explained, and the K2 and K, parameters are derived from 6 timesteps. Further, the error in K, is carried through into the estimation run for R2 and K, as these values may seem high. 59 after such a relatively long time period, that the time- dependency of the material is almost irrelevant. A parameter-estimation method of this complexity would be a much more efficient analytical tool for a material which exhibits much slower deformation time history (on the order of days or weeks). 60 6.2.1 Error in Material Models and Estimated Parameters Due to Loading History Deviations From the Assumed Heaviside Function Since the tensile testing fixture used with a material with an extraordinarily short deformation history is essentially an impact loading device, the assumption that the load onset in the sample forms a perfect step-function becomes somewhat questionable since the load in a perfectly elastic (undamped) system forms an oscillator with a load amplitude twice that of the dead weight loading, and a system which.is overdamped to the point.of not exceeding the applied load will show a load onset more closely akin to a ramp function. Both of these conditions violate the assumptions uponwwhich the material model is based, and.call the validity of the parameter estimation method into question. In order to check the load vs. time curve of the film in the test fixture, a strain gauge load cell was fabricated (Figure 12) and put in series with one of the test fixture tension cables (Figure 13) . The output of the strain gauge was conditioned with a Daytronic 3170 Strain Gauge Conditioner and displayed on a Hewlett-Packard 54504A Digital Oscilloscope. The load vs. time history for the first 500 milliseconds of loading was measured at each of 61 /1~xy8"umm Shadntki e /— lg on g... 3 Figure 12. Strain Gauge Load Cell 62 coaaoumcH Haoo coon suds ououxflm umoa Howxcwcc .ma r\\\\....... 3...... -— 33%. ll. 0 Rm 1. i... m / 5:933 g 53am .3W3.\ ounce 63 the three loadings used in the study. The loading was shown to be a ramp of approximately 90ms. duration (nearly one-half of the information recording period used for parameter estimation), with a peak load value approximately 25% higher than the dead-weight loading of the beams, and with a significant relaxation (Figure 14)‘. The effects of these deviations from the assumed load vs. time curve on the estimated values for the material constants are worth noting more for their implications in the method than for this particular study. Since the deflection of the material used in this study occurs over such a short period of time and since the R, value thus dominates the material model, the error in the K, estimate may be corrected directly as a function of the peak loading value. In the case of an experiment with deformation occurring over a more substantial length of time, the K1 and IQ variables become more significant, and the error induced in these coefficients may become more significant as well. It may be that for an extended load-time relationship the "ramping" and relaxation would constitute a fairly small component of the overall time history of the material; the major cause for concern would be, once again, error in the K1 estimate caused by the peak loading being larger than the ‘It should be noted that the curves shown are smoothed significantly from the actual data. which showed significant electronic noise as well as resonance data from the tensile fixture cables. 64 cuouxwm mcflumoa oaflmcoa cfi mafia .m> coon :HEORRMHHBE mafia com ocv com com - p P omoa Ofluoum .mnH m¢.m » _ _ _ _ . _ _ _ . _ Goon Oflumum .mna mo.m coon Oflumum .mQH mm.h “ . . . . _ coxou mo: sumo ucoEoooHamflc cowcz uo>o cofluomlllilv_ . .ea magmas cod . 0H ma om (mmmm) andnno IIao 65 static level on the beam. A solution to the problems caused by the real-world loading history compared with the "perfect" Heaviside step-function would be to incorporate matching variable loading conditions in the finite-element models used as objective functions in the parameter- estimation method and as a modelling method for the strained film. It is unclear, however, at what point the deviation from the assumptions necessary for the development of the material model will begin to affect the accuracy of the analysis. A more practical solution might be to redevelop the material model using a ramp loading function, a combination of ramp and Dirac spike functions, or some other (more accommodating) model rather than a step loading function. 66 6.3 Strain Within the Specimen Finite element simulations of the test specimen were constructed and run to ascertain the state of strain in the film, particularly within the region tested for permeation and morphology changes (Figures 15 and 16) . The values returned by the material parameter estimation experiments were used as the material constants in previously described VISC02.FOR software. The model returned values which were then plotted to give an estimate of the state of strain and change of free volume occurring within the material at the time of measurements This model provides algood ”map” of the area of interest in terms of strain and change of volume, and the change of volume and thinning during the three leading regimes is plotted in Figures 17 to 22. As the figures clearly illustrate, there is an increase in the free volume of the material in the region of the specimen where the permeation testing took place. 67 GENERAL SEQUENCE OF PROGRAMS ( Material Deflection Data‘) [bunence.rg§] Sue-ad Dlsplaceaant Data Saspla Grid P21A.F K. Estisata (Elastic Constant) [22w Estlsata of Material Constants ‘K.. K.’ K.) [v l sc02. refit—- Material Deflection Calculations Change in Thickness AREA.FOR and Free Volume Figure 15. General Sequence of Programs Used in Parameter Estimation and Material Response Calculations 68 .QBHQw oofim ow women-5 can 8.555. 95338 a... Bo: Bo 0:208 038E . S 866E .28 .m.m .8... .9... .mA .8 — p h _ p p u m .N «a .1 low... lac as 3.8 a: 3d 2: 2 4. so as E Aomv . AGNV ow ha o~ .. so so a. :5 so 8: 5 2 m :3 a: :3 an «N mu m m 88 :5 .2 2 82 a: o E a: m~ m a .2 69 Tucoouoo 5.” who moan.» :fi .mcacooq wooa um ENE HON mcofluodooamu 8255. .5 0.59.5 h m n v n N p o P p _ [TO //N n.0,. f £an 302 05 «o 8308 S a 35.0891 3:. .9 v.0: 70. 39.3 132 05 «o COMUUOm v) a 3885 mg “9 Tucoouom 5... one menace, :3 .9688 mum on SE 08 9802538 8:255. .3 086E o n w n N F o — F b b b (HID 8.0] (\Ii’J 14¢ . R U0 In . s as . 3 o. o. O 6‘ 6.] (Au-N Oi] billt'g’l 71 39.8 352 05 no c0308 v) a 3:08.53. amt. ”9 7900qu 5 one menace, HHS .9803 who on 5.3m Home 9.0303300 $5. a v n N . 3 086E 72 Tucmouwa CH mum 83$! :5 65an 33 um 838 wwb 5 $5.6 ESBHS .8 $ng 73 Tucwouwm 5 mum 83mg 3.4V .mcémS «mm um 839 8E fl mgfi Bumasufio AN 33E h o n v n N p o Rx? /w/m ”NH/A/WWwfl. llll.\l..W \. A. w .firmu . - I n 74 Soimfissy.9nuo c0308 v) a 3:08.53. 35. "9 A.u:mouwm :a mum mwsam> Hdfiv .9588 E um 239 8Q 5 9986 8,3328 .mm 93$ 6! VI V 75 6.4 Permeation Changes in the Material The permeation rate of three different samples of the material were measured at three different strain levels each. The relatively small degree of strain occurring in the sample suggests that there would be a marked increase in the rate of permeation, and this was shown to be the case. Additionally the increase was not linear (Figure 23) but the increase was monotonic with increasing strain. These findings are of considerable interest, since the permeation rate of a material apparently can increase by at least 20% due to small loads in what might be considered the "elastic" range of the material. Although this increase is ostensibly limited by the onset of fibrillar orientation in the material, there is still the potential for this phenomenon to be of interest in the engineering of polymeric material structures, either in the construction of a safety factor where barrier properties are the foremost concern, or in the manipulation of a polymeric film to the desired permeation specifications. The curvlinear nature of the increase in permeation suggests that the nature of the mechanism by which the permeation increases in the polymer is either extremely sensitive to changes in free volume, or that there is some 76 .goo~“ P .mHgEmm ou nowagmd mmmHum.mcmum> mumm coflummsumm mo gamma A.flmmv mmmuum ooflaomc HKXUF .mm mafia Au 49 fl.096 r Immo r Inxuw e ruaw. i Teavp umxwp Imzwp .uop.p fur; ('“QVM/ zN/UTW/ 8‘10) 9393 nonmci 77 threshold level within the material above which the material has a much higher level of permeation. In the latter case, the permeation increase within the material tested would vary as the threshold-exceeded area increases, giving an approximately second-order permeation curve with.a biaxially loaded sample. 78 6.5 Small Angle Light Scattering Analysis The Small Angle Light Scattering (SALS) method used to test the film for changes in orientation showed that the material had a distribution of spherulites on the order of Sum which were oriented in the ”machine direction‘" of the film. This low level of orientation is to be expected due to the stresses applied to all films during production. Unfortunately, even the largest strain applied to the material failed to show any appreciable difference in orientation direction or magnitude (as shown by the direction and degree of extension of the "flare" in the pattern in the photographs) between the strained and unstrained specimens in either the H, (polarizers crossed) or H, (polarizers parallel) configurations (Figures 25 to 28). This was expected from the small-strain model which suggests that the changes which occur in the material are due solely to free- volume changes in the structure in the polymer rather than the formation of morphological artifacts which would significantly alter the optical activity of the material. 6Machine direction refers to the direction in which the material is rolled up on a spool or otherwise transported through processing machinery. Almost all materials show some degree of structural orientation along the machine direction unless the process has been designed to eliminate these. 79 > ES 829% ~32 gums m mam .mm 953m @mQH conga? muowwm cumuumm >m gm .vm whomflm 80 ES Edema H32 fiwfimm Mm aw . hm 95mg UmQH 8:QO whomwm gown mm gm .3 353m 81 6.6 Scanning Electron Microscope Evaluation The electron photomicrographs of the surface features of unstrained and maximally strained film showed no significant changes in the surface features of the film between the strained and unstrained state (Figures 29 and 30). Again, this is to be expected from a phenomenon which does not alter the morphological nature of the material. The apparent difference in surface texture is largely due to the difficulty in focusing the SEM on an almost featureless surface. An interesting observation to be made from these photomicrographs is the existance of what appear to be "pores" in the material, although these do not change with the strain in the material and.may'or’may not be artifacts of the replicating process. 82 .coflgnommp 80wa E3.“ mo cmwggonm . mN wash 83 .coflgomwp kumm gnaw mo fimgggm . mm wfififlm 7.0 SUMMARY AND CONCLUSIONS 7.1 Summary The preceeding study developed a numerical model of the mechanics of deformable solids to estimate the thinning and free-volume changes that occur as a polymer film sample is stretched under a "stepped" loading regime. A tensile testing fixture was developed to load the film to various degrees of tensile stress, and to measure the movement of index marks applied to the test samples for the elucidation of material property constant estimates for use in the numerical models. The tensile fixture was also used to hold the stressed film for measurement of permeation and changes in optical activity once the loading was applied to a sample. A permeation cell was modified for use with the abovementioned tensile tester to measure the C02 permeation rate at the center portion of the film sample. Further testing of the center region of the sample was accomplished using a Small Angle Light Scattering apparatus developed specifically for the study, and by casting film surface replicas for analysis by Scanning Electron Microscope. With these tests and methods, it has been shown that the linkage between applied force.and changes in permeation can 84 85 be quantified, with a much better understanding of the specific changes that are occurring in the sample. 7.2 Conclusions It is clear that a large change of permeation is possible with relatively small dimensional changes, as predicted in the literature and shown in this study. Further, the changes may be correlated to a close estimate of the state of strain, thinning and volume change of the material at the time of ‘permeation. although the exact microstructural mechanism by which these changes occur is not elucidated by the analytical methods used here. The utility of this type of analysis is immediately apparent; it is now possible to predict the changes in permeation occuring within a polymer film structure under mechanical loading without construction of a prototype and without testing beyond that which is necessary to ascertain a few simple material properties. Extension of the fundamental principles underlying this study--the construction of a quantitative linkage between mechanical input and barrier property changes--should allow more rapid, and therefore more cost-effective, design and development of materials, production processes and structures using these materials. 8 . O RECOMMENDATIONS The study described here has proven the feasibility of the linkage of several types of studies to produce a clearer picture of the circumstances surrounding the changes in permeation in strained material. To be a more practical engineering tool several improvements are in order: The collection methodology with which the material- property data was obtained must be improved. The method has the potential to be quite accurate, but it is hampered by the inaccuracy and error-producing tedium of visual measurements and the low resolution of the video system used. A high- speed imaging and digitization system would be very useful. A better understanding of the changes in strained polymer structure is necessary.‘ Since gross structural changes are not occurring within the polymer, a more sensitive measure of the changes occurring (such as wide- angle x-ray scattering) may provide additional information. A more robust tensile strain fixture is necessary if high-speed deformations are to be studied since the one used in this study exhibited a great deal of secondary motion when activated. A load-deformation time history needs to be recorded to give a more accurate picture of the load response of the 86 87 film, and perhaps a more complex material model is in order to accomodate the less than ideal load-time history produced by the equipment used. APPENDIX A Program Listings This appendix contains listings for the following main programs and (limited) supporting documentation: VISC01.FOR VISC02.FOR P21.FOR P21a.FOR DUMERGE.FOR File Input Format for VISC01.FOR and VISC02.FOR Additional File Input Format for P21.FOR and P21a.FOR 88 89 90 Table 3. VISC01.FOR PROGRAM VISCOI . FOR DECLARE COHON STATEMENTS DWELE PRECISIGI “8.8) DMLE PRECISICNI INVRS(200,200) DMLE PRECISION KINV(200,200) COM/MINNELS MINNBLOCK/NWNCDES MLE PRECISIN “200) MIRBLOIXI R DMLE PRECISIOI MFFM) COW/CfiFFBLOCK/COEFF MIIDINBLOCK/IDIN cm lNTS/NWBEROFTINESTEPS DWELE PRECISIGI TINSTART, TININCR DMLE PRECISION NWXQOO) , NO0Y(200) WINWEINWXJWY INTEGER EI(325),EJ(325),EK(325),EM(325) CONNON/ELNODES/EI,EJ,EK,EH DMLE PRECISIGI NASTERK(200,ZOO) COIN/BI GK/NASTERK DMLE PRECISIGI A(3,3) COINNUABLOCK/A INTEGER NWSPYK WINS/NWSTART , NUNSTOP INTEGER IBDHZOO) MIIBlNDX/IBDY DMLE PRECISIOI BVAL(200) MIBWNDVAL/BVAL INTEGER WY WINDY/MDT INTEGER IDPLUS‘I WI IDP/ IDPLUS‘I MLE PRECI S I GI (PAST ( 200 , 200) WIOASTBLOCK/KPAST 0000000 91 Table 3 (cont'd) DOUBLE PRECISION DELU(200) INTEGER IROUCZDD),JCOL(ZDO),JORD(ZOO) CONNON/IJJ/IROH,JCOL,JORD DOUBLE PRECISION KZERO(ZOO,ZOO) CONNON/KZ/KZERO DOUBLE PRECISION Y(ZOO) CONNON/NYE/Y DOUBLE PRECISION FINALK(ZOO,200),LASTR(ZDO) CONNON/ENDO/FINALK,LASTR SAVE 100 FORMATIA) If (NFE .ED. 1) THEN “‘15 (*,100)' eeeeeeeee pmv!sco 1 new" ENDIF m AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA OPEN FILES This block of cmnds opens, labels, and nunbers the appropriate files for use by VISCOI with the exception of the series Of files needed for [K(t)] storage, as those are created as needed. OPEN (3, FiLEa'GENERAL_DATA', srArus- 'OLD') RENIND 3 OPEN (9, FILEI'BOUNDARIES', STATUS= 'OLD') REHIND 9 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA READ IN INITAL DATA This reads in some of the necessary parameters to operate ease of the arrays used in this program. READ (3,‘)NUNELS,NumberOfNodes NUMNODESslumberOfNOdes*2 NODCOUNT=NUNNODES NUNN=NUNNODES READ (3,*)(COEFF(I).I=1,6) READ (3,*)NUNBEROFTINESTEPS READ (3,*)NUNSTART,NUNSTOP READ (3,*)TINSTART,TININCR 92 CALL MAKEARRAYS(IT) 130 FORMAT (I2,15x,012.6> IF ((NFE .EO. 1) .OR. (NPRNT .GE. 1)) THEN NRITE (*,100) ’ COEFFICIENT VALUE' DO J31,6 HRITE (*,13D)J,COEFF(J) ENDDO ENDIF C Read in the known boundary conditions C from the 'boundaries.dat file: IF (NFE .E0. 1) THEN URITE (',100)' ' URITE (',100) ' NUMBER OF KNOWN DISPLACEMENT VALUES:' ENDIF READ (9,.) NUMBDY IF (NFE .EO. 1) THEN UNITE (*.*) NUMBDY URITE (',100)' ' URITE (*,100) ' DIRECTION INDICATOR: X81 YIO' HRITE (',IDD) ' NODE DIRECTION VALUE' ENDIF oo I-1,Nunsov READ (9,*) IBNDX, lDlR, BVAL(I) IF (NFE .E0. 1) THEN URITE (*,*) IBNDX,IDIR,BVAL(I) NRITE (*,100)' ' ENDIF C IDIR: XII YSO FOR Z'D PROBLEMS IBDY(I)=(IBNDX‘Z)'IDIR ENDDO CLOSE (9) C " Once arrays are ready, construct the series of [K(t)] values. c for all of the timesteps in the problem. if (NFE .E0. 1) THEN HRITE (*,100)' STORING [Kl MATRICES FOR ' ENDIF DO s,NT=o,NunsERorTINESTEPs IF (NFE .20. 1) THEN NRITE <*.101)' TINESTEp-I,NI ENDIF 101 FORMAT (A,IZ) CALL MAKEA(NT,TIMINCR,TIMSTART) C I 0 CALL MAKESNAPES( MASTERK,NUMNODES) ’ O O 0 CALL STOREMASTERKINT,MASTERK,NUMNODES) 93 5 CWTINUE C * Solve the tine-dependant problem, once all of the m values are c ready and stored. NNPLUSISNUMNODES+1 IDIMINUMNODES-NUMBDY IDPLUSIIIDIM+1 CALL SOLUTION(NUMNODES,NNPLUSl,KZERO,lDlM,lDPLUS1,FINALK, c LASTR,KINV) "RITE (’,100) ’ SOLUTION COMPLETED!ll!!' HRITE (*,100) ' ' CLOSE (3) 94 c AAAAAMAA‘AAAAAAAAAAAMAAAAAAAAA AAAA AA AAAAAAAAAAAAAAAAAAAAAMAAAAAAAAAAAAA SUBRQIT I NE MAKEARRAYSI IT) C This subroutine sets up the indexed array of node and ale-ant values used C by the [K] generation loops. INTEGER NODNO,IELNO COMMON/NUM/NUMELS DOUBLE PRECISION NODX(200),NOOY(200) MleE/wa,NmY INTEGER EI(325),EJI325),EK(325>,EM(325) COMMON/ELNODES/EI,EJ,EK,EM OPEN (’9, FILEa'ELENENT_oATA', sums: 'OLD') REIIINO I. C * Create an array of indexed x & y values associated with each node 100 20 26 FORMAT (A) IF (NFE .E0. 1) THEN WITEPJDO) ' NODE X Y' ENDIF CONTINUE READ (4,*) NODNO IF (NODNO .EO. '1) GO TO 23 READ (4,*) NODX(NODNO), NOOY(NODNO) IF (NFE .E0. 1) THEN NRITE (*,') NODNO,NODX(NODND),NODY(NODNO) ENDIF SS TO 20 IF (NFE .E0. 1) THEN URITE (*,IOO) ' ' URITE (*,100) ' ELEMENT I J K "I ENDIF CONTINUE * Produce an index of nodal values associated with each element 25 READ (6,*) IELNO IF (IELNO .E0. -1) GO TO 25 READ (6,*) EI(IELNO), EJ(IELNO), EK(IELNO), EM(IELNO) IF (NFE .E0. 1) THEN HRITE (*.*) IELNO,EI(IELND),EJ(IELNO),EK(IELNO),EM(IELNO) ENDIF GO TO 26 CONTINUE RETURN END C - C 100 101 A‘AALAAAA---AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA‘AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAQ SUBROUTINE MAKEA(NT,TIMINCR,TIMSTART) DOUBLE PRECISION TIMEVAL,TIMINCR,TIMSTART DOUBLE PRECISION A(3,3) DOUBLE PRECISION COEFF(6) DOUBLE PRECISION GI,G2,B,C,D COMMON/COEFFBLOCK/COEFF COMMON/ABLOCK/A FORMAT (A) FORMAT (A,I2) C Calculate the value of the timestep in seconds: an 000 C C 40 flflflflfi O TIMEVAL=(NT*TIMINCR)+TIMSTART NRITE (*,') TIMEVAL URITE (*.100) ' ' GO TO 30 This Subroutine computes the values for E and MU at some tine-step NT. It utilizes a 3-parameter exponential decay model for the elastic modulus and Poisson's ratio. The COEFF array contains the necessary coefficients. EMSCOEFF(I)tCOEFF(2)*DEXP(-I.ODO*(COEFF(3)*TIMEVAL)) PR8COEFF(6)+COEFF(5)*DEXP(-1.ODO*(COEFF(6)'TIMEVAL)) This part takes the E and MD values for the current time value and returns the (A) matrix for that time value. BtEM/(1.0DO-(PR**2.0DO)) A(1,I)-B A(2,2)=B A(3,3)8B*(1.ODD-PR)/2.ODO A(1,2)=PR*B A(2,1)8A(2,1) A(1,3)=O.DO A(Z,3)'O.DO A(3,I)=O.DO A(3,Z)=0.DO UNITE 0,100) ' [A] MATRIX VALUES’ no Him-1,3 no JCOL=1,3 UNITE (*,*) IROU,Jc0L,A(IROU,JcOL) ENDDO ENDDO CONTINUE RETURN END SUBROUTINE MAKESHAPES(MASTERK,NUMNODES) DDUDLE PRECISION NASTERKINUNNDDES,NUNNDDES) DOUBLE PRECISION K(8,8) COMMON/NUM/NUMELS INTEGER EIC325).EJ(325),EK(325),EM(325) COMMON/ELNODES/EI,EJ,EK,EM 100 FORMAT (A) IEIGHT-B C -This Subroutine through each of the elements and (using the proper c subroutine) develops an elemental [k] matrix to be added into the [K] via c the MERGER subroutine. C * Put the two together and route to appropriate calculation of IkJ for C each element then add the value for that element into the global [K] C But first, the MASTERK must be cleared from the last tile-step. Do II-1,NUNNDDES Do JJ=1,NUMNODES NASTERKIII,JJ)=0.DDD ENDDO ENDDO C Proceed to assemble next MASTERK DO 30,1Element81,NUMELS IF (EM(IElement).E0.0) TNEN C O I CALL TRIANGLELEM(IElement,K) GO TO 27 ENDIF C O I CALL SOUARELEM(IELEMENT,K) 27 CONTINUE C O 0 CALL MERGEK(K,IElement,MASTERK,NUMNODES) C Note that MERGER merges the element's [K] value into the COMMON MASTERK() C and does not return any value. 30 CONTINUE RETURN END 597 c AAAAMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA C I O SUBROUTINE TRIANGLELEM(IEIement,K) DOUBLE PRECISION K(B,B) DOUBLE PRECISION NODX(200),NODY(200) COMMON/NODEINODX,NODY INTEGER EI(325),EJ(325),EK(325),EM(325) COMMON/ELNODES/EI,EJ,EK,EM DOUBLE PRECISION A(3,3) COMMON/ABLOCK/A DDUDLE PRECISION XI,XJ,XK,XM,YI,YJ,YK,YM,BI,BJ,BK,CI,CJ,CK DOUBLE PRECISION x<3>,v<3>.a<3,6),C(6,3),AR2,suI C This subroutine uses the brute-force calculations in TRIANGLECALC C to produce a [kl for the selected element. 100 FORMAT (A) c 09...... THII.ODO c tiflflfiiii XIINODX I:’,I URITE (*,100)’ '[K] ' (DU)' NRITE (FILEDELU,70) 'DELUFILE’,I JINT'I UNITE (FILEKNUM,75) 'KNUMBER',J UNITE (*,*) J,I FORMAT(A,IZ) FORMAT(A,I2) OPEN (IZ,FILE=FILEDELU, STATUS=’OLD’) OPEN (15,FILE=FILEKNUM, STATUS=’OLD') RENIND 12 RENIND 15 DO L-1,NUNNOOES DD M81,NUMNODES READ (15,*)KPAST(L,M) ENDDO ENDDO READ (12,*)(DU(KK),KK=1,NUMNODES) CLOSE (12) CLOSE (15) C C C C C C C C 0000 000 00 000 00000 000 0000 11353 DIAGNOSTIC OF FILE-READ KPAST AND DU NRITE (*,100) ' OUTPUT OF READ VALUES OF [KPAST] ' NRITE (*,100) ' ID JD KPAST (ID,JD)' DO ID-1,NUNNOOES DO JD=I,NUMNODES UNITE (*,*) ID,JD,KPAST(ID,JD) ENDDO ENDDO DO ID81,NUMBDY URITE (*,101) ’ BOUNDARY AT ROU:',IBDY(ID) ENDDO NRITE (*,100) ' ID DU(ID) AS READ' DO IDSI,NUMNODES URITE (*,*) ID,DU(ID) ENDDO Collapse the [KPAST] and (DU) matrices according to the boundary conditions associated with the (DU)'s timestep index. Do (DU) first, producing (DUINTER): Then This Rous Then IFLAG-O DO II=1,NUNNOOES IF (II .50. IBDY(IFLAG+1)) THEN UNITE (*,101) ' SKIPPING Du ROU',II IFLAG=IFLAG+1 ELSE DUINTER(II-IFLAG)=DU(II) UNITE (*,101) ' COPYING DU ROU’,II UNITE (*,101) ' T0 DUINTER ROU',lI-IFLAG UNITE (*.*) DU(II) NRITE (*,100) ’ ID DUINTER(ID)' DO ID-1,IOIN UNITE (*,*> IO,OUINTER(IO) ENDDO [KPAST] producing [KINTER]. is a nested sort similar to the BYVAL() subroutine. first: IFLAG-o DO II-1,NUMNOOES IF (II .50. IBDY(IFLAG+1)) THEN IFLAG=IFLAG+1 ELSE columns: UNITE (*,100) I MYSTERY GLITCH, NUMNODES IS' UNITE (*.*) NUMNODES UNITE (*,100) JJ, JFLAG, IBDY(JFLAG+1) JFLAG=0 DO JJ=1,NUMNODES 109 C URITE (*,*) JJ,JFLAG,IBDY(JFLAG+1) IF (JJ .EO. IBDY(JFLAG+1)) THEN JFLAG=JFLAG+1 ELSE KINTER(II'IFLAG,JJ-JFLAG)=KPAST(II,JJ) ENDIF ENDDO ENDIF ENDDO C DIAGNOSTIC OF COLLAPSED KPAST AND DU MATRICES (KINTER AND DUINTER) C UNITE (*,100) ' OUTPUT or VALUES or [KINTER] AND (DUINTER)' C UNITE (*,100) ' ID JD KINTER (ID,JD)' C DO ID-I,IDIM C DO JD=1,IDIM C UNITE (*.*) ID,JD,KINTER(ID,JD) c ENDDO C ENDDO C UNITE (*,100) . ID DUINTER(ID)' C DO IDsI,IDIN c UNITE (*.*) ID,DUINTER(ID) C ENDDO C Multiply IKINTER] by (DUINTER) to produce a "coupressfi" (RINTER) CALL SOUARBYCOL(IDIM,KINTER,DUINTER,RINTER) DIACNOSTIC OF RINTER URITE (“.100) ' AFTER SOUARBYCOL:’ URITE (*,100) ' ID RINTER(ID)' DD Io-1,IOIN URITE (*,*) ID,RINTER(ID) ENDDO 000000 C "Unpack“ the (RINTER) values into the appropriate rOUs of (RESID) IFLAC=0 DO II=1,NUMNOOES IF (II .50. IBDY(IFLAG+1)) THEN If this is a “deleted" row, regenerate the (RESID(II)) value from the appropriate row and column of [KPAST] and (DU) URITE (*,100) ' REGENERATION CHECK OF KPAST VALUES’ URITE (*,100) ’ ID JD KPAST (ID,JD)' DO ID81,NUMNODES DO JD=1,NUMNODES URITE (*,*) ID,JD,KPAST(ID,JD) ENDDO ENDDO 00000 00 00 SUM=0.0DO C URITE (*,100) ' IBDY (IFLAG+1) KK SUM KPAST DU’ DO KK=1,NUMNODES SUM=SUM+KPAST(IBDY(IFLAG+1),KK)*DU(KK) 000 0 000 000 110 URITE (*.*) IBDY(IFLAG+1),KK,SUM,KPAST(IBDY(IFLAGfl),KK),DU(KK) ENDDO RESID(II)=SUM URITE (*,101) ' REGENERATED R VALUE FOR R:',II URITE ('.*) RESID(II) IFLAG=IFLAG+1 ELSE If not, copy the value from (RINTER) and subtract the necessary coefficients (Uhich are carried through from the solution of the system of compressed equations). RESID(II)=RINTER(II-IFLAG) DO LL=1,NUMBDY RESID(II)=RESID(II)+KPAST(II,IBDY(LL))*DU(IBDT(LL)) ENDDO Finally, subtract the residual terms generated for this series of arrays from the original force vector (R) ENDIF ENDDO UNITE (*,100) ' (R) VECTOR JUST BEFORE SUBTRACTION' 00 I081, NUMNODES UNITE (*.*> ID,R(ID) ENDDO URITE (*,100) . N R(N) RESID(N)’ DO N-1,NUNNOOEs R(N)=R(N)-RESID(N) URITE (*,*> N, R(N), RESID(N) ENDDO CONTINUE C Copy whatever is left into (FINALR) for return from the subroutine. DO N81,NUMNODES FINALRIN)=R(N) ENDDO RETURN END 111 c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA SUBROUTINE STOREDELU(DELU,NT,NUMNODES) DOUBLE PRECISION DELU(200) CNARACTER*IS FILEDELU C This subroutine stores the calculated values of (/\U) in individual files CONMON/NUM/NUMELS INTEGER IBDY(200) COMMON/IBINDX/IBDY DOUBLE PRECISION BVAL(200) CONMON/BOUNDVAL/BVAL DOUBLE PRECISION DUMMY(200) INTEGER NODEHALF,NNUM Reinsert known boundary values in the (DELU) array. NOTE the (DUMMY) vector is the "unpacked" solution vector which also contains the known boundary conditions at the particular time-step. 0000 IFLAG=0 DO ll=1,NUMNODEs IF (II.EO.IBDY(IFLAG+1)) THEN DUMMY(lI)=BVAL(IFLAG+l) IFLAG=IFLAG+1 ELSE DUNMY(II)=DELU(Il-IFLAG) C Truncate very small values to avoid underflow errors: 00 KK=1,NUMNODES IF (DABS(DUMMY(KK)) .LE. 1.00-24) THEN DUMMY(KK)=0.0DO ENDIF ENDDO IF (NT.NE.0) THEN c URITE (*,100) ' STOREDELU RUN ' 100 FORMAT (A) URITE (FILEDELU,85)’PARADELUFILE',NT 85 FORMAT(A,IZ) OPEN(22,FlLE=FILEDELU,STATUS=’NEU') C Nrite to appropriate file. 00 I=1,NUMNODES URITE(22,*)DUMMY(I) ENDDO CLOSE(22) 135 FORMAT (I2,014.6,014.6) IF ((NFE .E0. 1) .OR. (NPRNT .GE. 2)) THEN C Screen output for instant gratification. URITE (*,100)' ' URITE (*,85) ' DISPLACEMENT VALUES FOR TIMESTEP',NT UNITE (*,100) ' VALUES LESS THAN 1.00-25 ARE STORED AS 0.000' URITE (*,100) ' NODE DX 112 Do ID-1,NUMNODES,2 JOIID+1 NNUM=(JD/2) URITE (',*) NNUM,DUMMY(ID),DUMMY(JD) ENDDO ENDIF ENDIF URITE(*,100)' **' RETURN END ILJL3 c AAAAAAMAMAAAAA A A A A A AAA AA AA AA A A A AAAAA AAAAAAAAAAAAAAAAAAAAAAAAMAAAAAAAAAAAAA 000000 100 SUBROUTINE BTVALINUNBDT,NUHNOOES,NASTENK, FINALR,IDIM,LASTR,FINALK,IDPLU51,LSIG) DOUBLE PRECISION MASTERK(NUMNOOES,NUMNODES) DOUBLE PRECISION FINALK(IDIM,IDPLUSI) INTEGER IBDY(200) COMMON/IBINDX/IBDY DOUBLE PRECISION BVAL(200) CONNON/BOUNDVAL/BVAL DOUBLE PRECISION FINALR(200) DOUBLE PRECISION LASTR(200) DOUBLE PRECISION TEMP This subroutine takes the appropriate boundary conditions from the IBDY (indexed list Of ROUS which have known boundary conditions) and BVAL (the values associated with the known boundary conditions) and substitute the proper values into the MASTER! array at each time step. Note that the boundary conditions are not time- dependent in this version. FORMAT (A) c ZERO OUT THE APPROPRIATE Rows IN [MASTERK] AND (FINALR) Do III,NUMBDY DO J-1,NUMNODES MASTERK if the C. Itoh printer is not working ' 8eb162 _imagen for NO printout show queue eb...... gives printer stack-up. FORTRAN/LIST (FILENAME) LINK LINK .obj,IMSL/LIB To dump output to a file: DEFINE SYSSOUTPUT DEASSION SYSSGJTPUT To MERGE two files: COPY , lLOC 0 01560000000600 flflflflflflflflfifl 118 Table 4 (Cont’d). Note that the I/O convention applies to each subroutine as an independent entity. From this, in a subroutine, values come in (I) and others are returned (0). Subroutines which call another subroutine: Hhen the call is made, values are sent out (0) and others return to the caller (I). This is similar to dolble-entry bookkeeping in that there should be a match between I and O designations throughout the program. AAMAAAAAAAAAAAAAAAAAAMAAAAAAMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA DECLARE comou STATEMENTS INCLLDE ' KSTWEJM' DCIJCLE PRECISIGI K(8,8) MIMI/mu MINNBLxK/NWNQES DOJBLE PRECISIN RlZOO) MIRDLWR DMLE PRECISIGI CGFHOI MIMFFDLOCK/CNFF DQJBLE PRECI SIOI TN WITNICKBLOCKITN COMMON/IDIMBLOCK/IDIM COMMON lNTS/NUMBEROFTIMESTEPS DOUBLE PRECISION TIMSTART, TIMINCR DWBLE PRECISICNI NCDXQOO) .NWYlZOO) WINwE/me,NwY INTEGER EI (325).EJ(325).EK(325),EM(325) CGDiON/ELNwES/EI ,EJ ,EK,EM DWBLE PRECISIQI MASTERK(200,200) COHON/BIGK/MASTERK DWBLE PRECISICNI “3.3) CORN/“LOCK/A INTEGER NIHSPYK WINS/”START , NlIISTOP INTEGER IDDY(200) COIN/IBINDX/IBDY DMLE PRECISICNI BVALIZOO) CGINNUBGJNDVAL/BVAL INTEGER NLHDDY CGIKNUBDY/NlMBDY INTEGER IDPLUS‘I WI IDP/ IDPLUSI 119 DMLE PRECI SIN KPASTQOO, 200) WIWASTMWKPAST INTEGER IRW(ZDO),JCOL(200),JWD(ZOO) WIIJJ/IRWJCOLJMD DMLE PRECI SIN KZERO( 200 , 200) WIHIKZERO DMLE PRECI SIN “200) MAME/Y MILE PRECISIN FINALK(200,200) , LASTR(200) WIENDO/FINALK,LASTR INTEGER IBANDNIDTH WIIBU/IBANDNIDTH INTEGER ISTMFLAG MIISF/ISTORFLAG SAVE 100 FNMAHA) “11E (i'100)l eeeeaem Vlsco 2 W I URITE (*,100)' v1.7 ' RITE (*,100)' SCOTT IDRRIS ' RITE (*,‘IOO)' Michigan State University 1991 ' URITE (*,100)' ' MAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA WEN F I LES This block of con-Dands opens, labels, and mmbers the appropriate files for use by VISCOi with the exception of the series of files needed for (KR)! storage, as those are created as needed. nonnnnn WEN (3, FILE"GENERAL_DATA', STATUS. 'OLD') RENIND 3 (PEN (9, FILEII'DCAINDARIES', STATUS8 'OLD') RENIND 9 MAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA c C C READ IN INITAL DATA C This reads in some of the necessary parameters to operate some of c the arrays used in this program. READ (3,*)NuIELs,Nunbeu-0fuodes NllinESINLIIberOfNodes*2 NUNN-NuiNODEs READ (3,*)(COEFF(I),I-1,6) READ (3.*)NLHBEROFTIMESTEPS READ <3,*)NIms1ARt,NuISIOP READ (3,*)TIMSTART,TIMINCR READ <3,*)ISIORFLAC CALL NAREARRAvs c mmmeeeeeeeeaea C (INSTANT THICKNESS TERM: JJZC) THII .1130 c Win...” URITE (',IOO) ' COEFF 1 --> 6' DO J-1,6 URITE ('.')J.COEFF(J) ENDDO CLOSE (3) C Read in the known boundary conditions from the 'boundaries.dat file RITE (*,100)' ' RITE 0,100) ' HUBER OF KNOJN DISPLACEMENT VALUES:' READ (9)) MDT URITE (*.‘) NUMDDV RITE (*,100)' ' RITE (',100) ' DIRECTION INDICATGI: X81 Y-O' RITE 0,100) ' ME DIRECTIGI VALUE' Do I-I,NUNst READ (9.*) IsNDx, IDIR, BVALlI) wRIIE (*,*) I8NDX,IDIR,BVAL(I) wRIIE (*,100)' . C IDIR: XII YIO FOR 2'0 PROBLEMS I30“ I )‘I IBNDX‘Z)‘ IDIR EUDO CLOSE (9) C * Once arrays are ready, construct the series of [Kttil values. C for all of the timesteps in the problem. RITE (*,100)' STORING no MATRICES FOR . DO 5,NI-O,NUNEEROFIINESIEPS CALL NAREA) C(oi-TN*((A(3,3)*8)/(6.ODO*AA)) C O I CALL RECTANGLECALC (C,k) RETURN END 14259 c AAAAAAAAAAAAAAAAAA A A A A A A A A A AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA SUBROUTINE MERGEK(K,IElement,MASTERK,NUMNODES) DOUBLE PRECISION MASTERK(NUMNODES,NUMNODES) COMMON/NUM/NUMELS DOUBLE PRECISION K(8,8) INTEGER S!(8) INTEGER EI(325),EJ(325),EK(325),EM(325) COMMON/ELNOOES/EI,EJ,EK,EM 100 FORMAT(A) -This Subroutine adds the [K] for some element into the global I!) for some time-step. -Master! is Global [K] for a given time step. -IEL Contains the node-list for the element. 'NEL Number of elements. nnnnn S!(1)=2*EI(IELEMENT)-I S!(2)-2*EI(IELEMENT) S!(3)=2*EJ(IELEMENT)-1 S!(£)-2*EJ(IELEMENT) S!(S)-2*E!(IELEMENT)-I S!(6)-2*E!(IELENENT) C This skips the EM for the triangular element to avoid C S!(7) and SK(8) = -1 and MERGEs the SOUARE element. IF (EM(IELEMENT).NE.0) THEN S!(7)-2*EMIIELEMENT)-I S!(8)=2*EMIIELEMENT) C eeeeeeeeeeeeee Merge into [K(t)] seeseeeeaeeeeeeeeeeeeeeeee DO 15,I=1,8 Do ID,J=:,S MASTERK(+c<4) !(2,I)-k(1.2> !(2,2)-2.ODO*(C(2)+C(5)) K(2,3)I-1.ODO*C(3)+C(4) !(2,b)8C(2)-2.ODO*C(S) N<2,5)=-I.OCO~ !(3.2)-K(2.3> !(3,3)-2.00“*(C(1)+C(6)) !(3,A)--1.0F”"C(3)+C(4)) !(3.5)8C(1)-;.‘OO*C(6) !(3.6)=C(3)-"3> K(3,7)l-1.OFI'(C(T)*C(6)) !(3.8)=C(3)*24«) !(k,i)=!(1,4‘ !(6,2)=!(2,Ai K“,3)=K(3,4; x<(.4)s2.oo‘~'¢<2)+c<5)) !(4,5)=-1.0F‘*7(3)+C(4) !(4,6)--2.0P'*C(2)+C(5) !(6,7)8C(3)+C(4) !(6,8)=-1.0C“' !(7,3)-!(3,7I !(7,4)-!(4,7) !(7,5)=!(S.7> !(7,6)-!(6.TE !(7,7)-2.008'6C<1)+c<6)) !(7,8)--T.OO"CCI3)+C(4)) !(D,1)I!(1.8) K(8.2)-!(2.?3 !(8,3)=!(3,“f ‘(83‘)'K(le , P) !(8.5)'K(5.9) K(3.6)'K(6.8) !(8,7)=!(7,8> !(8,8)82.ODR'C(2)+2.ODO*C(5) C >THICKNESS i“CTOR HEREIIIIIIIIIIIIIIIII 131 1332 c AAAMWAAAAAAAA A A A A A A A A A A AAAAAAAAAAAAAAAAMAMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA C C 100 O nor-Inn 150 50 55 SUBROUTINE STOREMASTERK(NT,MASTER!,NUMNODES) STORE the values for [K(t)) in a unique file, once the values have been calculated. INCLUDE ' KSTORE.FOR' DOUBLE PRECISION MASTERK(NUMNODES,NUMNODES) COMMON/NUM/NUMELS INTEGER IBANDUIDTH COMMON/IBU/IBANDUIDTH CHARACTER*IS FILENAME Note that this only stores the diagonal and upper values of the I!) matrix. IF IT AIN'T SYMMETRICAL IT'S GONNA BEII FORMAT (A) URITE (*,100)’ TEST DUMP BEFORE STORAGE IN STOREMASTERK' DO IT81,NUMNOFFS DO JT81,NH“V"DES URITE (*.7") IT,JT,MASTER!(IT,JT) ENDDO ENDDO forllt (I3,13,c‘2.4) '*********** Column Matrix Storage Conversion ***‘** ICOUNTER=1 DO 55,IRON=T,NHMNODES IF ((IRON‘IRANDUIDTH).GE.NUMNOOES) THEN MAXCOL=”””HODES ELSE MAXCOL='"’T~IRANDUIDTH'I ENDIF DO 50,ICOLrTROu,MAXCOL RITE <*,*> m, ICOUNTER, IRON, ICOL STOREDK(NT,Y"O”NTER)=MASTERK(IROU,ICOL) ICOUNTER=iCO”;TER+1 CONTINUE CONTINUE C ******* END OF COLUMN-MATRIX STORAGE ROUTINE ******** RETURN END 1J33 c AAAMAAAAAAAAAA A A a A A A A A A AAAAAAAAAAAAAAAAAAMAMAMAAAAAAAWMAMAAAAAA 100 000 000 SUBROUTINE SOLUTION (NUMNODES,NNPLUSI,KZERO,IDIN,IDPLUSI,FINALN, LASTR) INCLUDE ' KSTORE.FOR' DOUBLE PRECISION !INTER(200,200) DOUBLE PRECISION KZERO(NUMNODES,NNPLUSI) DOUBLE PRECISION FINALK(IDIM,IDPLU51) INTEGER NRC COMMON/NUM/N“”FLS COMMON/RBLOCR/R COMMON/BDY/NU“RDY COMMON/IBU/IRANDUIDTH COMMON INTS/NUNREROFTIMESTEPS DOUBLE PRECISION RIZOO),FINALR(200) DOUBLE PRECISION DELU(200) DOUBLE PRECISICN LASTR(200) DOUBLE PRECISION Y(200) FORMAT (A) This subrou':na retrieves the appropriate values for I '7)I, [dU] etc. and steps through the appropriate sets of solutions. Bu: first the [!(0)I values must be retieved. NTIO Minfitfi - s s v t ttfitiiiiittitfififiififltimimmm This unpacks '“O [K] matrix from the appropriate column of the STO“"""IME,COLUMN) array. NRITE (*,*: ““NODES, IBANOUIDTH ICOUNTER=1 DO IROU=1,NHMNOOES IF ((IRO”"iiflDNIDTH).GE.NUMNODE$) THEN MAXCOL?“ “2OOES ELSE MAXCOL=Z°TU+IBANOUIDTH-1 ENDIF URITE (*,1°“\ ' MAXCOL' URITE (*,*T “ VCOL DO ICOL=' ".WAXCOL “It: (t'1r‘"~ . tettttifi’il KZERO('”",ICOL)=STOREDK(NT,ICOUNTER) !ZERO(?” ,I?CN)=KZERO(IRON,ICOL) ICOUNTV":;‘CHNTER+1 URITE (*,1"‘ TRON,ICOL,KZERO(IROU,ICOL) ENDDO ENDDO c . f) f! flfoIfIfl fD-O (if? (1 f) CDC! 1134: *‘******** END OF UNPACKING ROUTINE ***““******* for-t (I3,i3,e12.b) write (*,100) ' test (Alp after retrieval routine in SOLUTION DO II-1,NUMNODES DO JJ-1,NlliNmES RITE 0,150) II,.IJ,!ZERO(II,JJ) flittiffiiittCitififitfifiiiiifitttfiittfiti0..titttitittfi Retrieve the (R) vector's value NOTE: Because of a glitch in VMS-FORTRAN (R) is passed via a comson statement. 00 NT-1,NI.MBEROFTIMESTEPS O 0 CALL FINDR(NT,NUMNODES) subtract the stresses from the previous timesteps I O O O O 0 CALL SUBRESIDUAL(FINALR,NT,!INTER) Accomt for the know bouIdary conditions 0 0 0 O 0 I I CALL BYVAL(NUMBDY,NUMNODES,!ZERO,FINALR,IDIM,LASTR,FINAL!, C IDPLUST,0) NOTE: FINAL! is the auusnted matrix containing the remaining siwltaneous «nations to be solved. Attsapt a solution: URITE (*,100) ' SIMULT CALLED, DEL-(U) VALUES FOLLON' CALL SIMULT(IDIM,FINALK,DELU,1.0D-ZS,1,IDPLUS1,Y) Store the results, AFTER re-inserting know bouIdary conditions. CALL STOREDELU(DELU,NT,NUMNODES) Print the results CALL PRINTDELU( ) ENDDO RETURN END 135 c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA SUBRGJTINE FINDR(NT,WES) DMLE PRECISIW "200) DCIIBLE PRECISIN VAL INTEGER INDEX,IDIR MIM/MLS WINS/NWTART ,MSTW CGIDN/RBLOCK/R CHARACTER*15 FILENAIE C This sLbroutine develops the value of (R) for timestep m:- NT (actual C time vslue of TIMEVAL(NT) seconds) and readies it for the slbtraction of the C resithal stress values. 100 FWMAT (A) IF (NT.GE.NlMSTART .AND. NT.LE.NIMSTW) THEN NNI‘I ELSE NN'O ENDIF This applies a stepped imut between NRSTART and NuIISTw time intervals IdIich rstpires the availability of RNRBERI and RNlNlBERO files. A Dirac spike has the same value for NUTSTART and NRSTW. Other irput shapes my be created by modifying the values of the anction and letting NN-NT daring the relevant time-frlae. flflflflfl RITE (FILENAME,50)'RNlliBER' ,NN 50 FWMAT (A,I2) WEN (20,FILEIFILENAIE, STATUSI'OLD') RENIND 20 c RITE (*,‘|00) ' FILENAME:' ' C RITE(",TOD) FILENAME C Set the IdIole vector to 0.000 00 I'1,NlliNCOES R(I)=0.0DO ENDDO Read the non-zero values of R from the appropriate file Remeaber: x-i, y-O on READ (20,")NlR4RVAL RITE 0,100) 00 J81,NlliRVAL READ(20,*) INDEX,IDIR,VAL NENI NDEX8( INDEX‘Z)’ IDIR R(NENINDEX)'VAL C ' RITE ('3') "DEX, IDIR, VAL, NENIIDEX ENDDO IF (NT .LE. 2) THEN RITE (*,100)' ' RITE (*,50) ' (R) VECTW FOR TIMESTEP',NT RITE 0,100) ’ NmE N Y’ RITE (*,‘|00)' ' ENDIF 136 DO IDI‘I,IAIIN(DES,2 JD-IDH NNlRiIJD/Z IF (NT .LE. 2) THEN RITE (’.') NNUI,R(ID),R(JD) ENDIF EIDDO CLOSE(20) RETURN END c AAAAAAAAAAAAAAAAAAA C 101 100 137 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA O I I SUBRGITINE SlflRESIDUAL(FINALR,NT,KINTER) This stbroutine slbtrscts the resichal force, from preview timesteps from the current (R) value. INCLIDE ' NSTWE . FW ' W IIDIIBLxK/IDIM WINNBLOCKINMES W/NW/NUMELS cm INTS/NLRIBEROFTIESTEPS MI I B I NDX/ I BOY MINT/HWY WI RBLOCK/ R WIKPASTBLWK/KPAST MI I Bill I BADNIDTH INTEGER IBANDHIDTH INTEGER IBDY(200) , NLHNWES DMLE PRECISIW DU(ZOO),DIHIY(200) DMLE PRECISIW R(200),RESID(200),FINALR(200) DGJBLE PRECISION RINTER(200),DUINTER(ZDO) DGIBLE PRECISION KPAST(200,200) DwBLE PRECISIN KINTER(IDIM,IDIM) DMLE PRECISIOI SUI CHARACTER'15 F I LEDELU CHARACTER'15 F I LEW FRMAT (A, I2) FWMAT (A) RITE 0,100) ' IDIM IN SUBRESIDUAL' RITE (',') IDIM RITE 0,101) ' TIMESTEP NulBER',NT NTPLUSISNT+1 IEND-NT-i eeeeeeeeeeeemmeeeeee MCEI‘IWCI Low W DO 80,I=O,IEND JsNT-I ”m This extracts [“0] from the appropriate column of the STmED!(TIIE, Element) array. ICGJNTER-i DO IRw-1,NuiNmES IF ((IRMIBANDNIDTH).GE.IAMNG)ES) THEN MAXCOLOWES ELSE MAXCOL-IRWIBANDUIDTHSI ENDIF DO ICOL=IRGI,MAXCOL 50 fifififif) fD-fi (if) CDCDCD 138 KPAST(ICOL,IROU)'STOREDK(J,ICOUNTER) KPAST(IRON,ICOL)8KPAST(ICOL,IRON) ICOUNTERIICOUNTER01 .fifitififififlflflfiiififl End of [‘(t)] Extr.ction *fitiitttti format (I3,i3,e12.b) write 0,100) ' test (up after retrieval routine IN sLbresithel' DO II81,NUNNODES DO JJ-1,NUMNODES URITE (*,150) II,JJ,kpeBt(II,JJ) ENDDO ENDDO eeeeee (dU(t)) Extraction eeeeeeeeeeeeeeeeeeeeeeee URITE (*,100) ' I KN DELUSTORED(I,KK)' DO KK!1,NUMNODES DU(KK)'DELUSTWED (1.x!) URITE (*,.) I,KK,DELUSTORED(I,KK) ENDDO ‘0'... End of .xtr.cti°n fflfifififtfitfifiitti Collapse the [!PASTI and (DU) matrices according to the boundary conditions associated with the (001's timestep index. Do (00) first, producing (DUINTER): Than This Rows Then IFLAGIO DO II81,NUMNODES IF (II .EO. IBDY(IFLAG+1)) THEN IFLAG=IFLAG+I ELSE DUINTER(II'IFLAG)-DU(II) ENDIF ENDDO I!PASTI prodacing IKINTERI . is a nested sort similar to the BYVAL() subroutine. first: IFLAG-O DO II81,NUMNODES IF (II .ED. IBDY(IFLAG+1)) THEN IFLAG'IFLAG+1 ELSE columns: JPLACno Do .IJ-I,NUNNmEs IF (JJ .Eo. IBDY(JFLAG+1)) THEN .IPLACaJPLACoI ELSE RITE 0,100) ' II IFLAG JJ JFLAG' RITE (an) II,IFLAG,JJ,JFLAG !INTER(II-IFLAG,JJ-JFLAG)8!PAST(II,JJ) ENDIF ENDDO 0000000000 0000 00“ 1J3!) ENDIF EIDDO RITE 0,100) ' BEFWE SGBYCOL’ RITE (*,101) ’ IDIM 8 ',IDIM DO JJ81,IDIM RITE 0.") JJ,DUINTER(JJ) DO KX'IJDIM RITE (',*) JJ,KX,KINTER(JJ,KK) ENDDO ENDDO Multiply (KINTERI by (DUINTER) to produce a “compressed“ (RINTER) CALL SOUARBYCOUIOIM,!INTER,DUINTER,RINTER) RITE 0,100) ' AFTER SDAREBYCOL CALL -C(IiPRESSED-' DO NN-I,NUNNODES RITE (*.*) NN,RINTER(NN) ENDDO “Unpack“ the (RINTER) values into the appropriate rows of (RESID) IFLAGsO DO II81,Nl.NiNCI)ES IF (II .EO. IBDY(IFLAG+1)) THEN If this is a "deleted'I row, regenerate the (RESID(II)) value from the appropriate row and calm of IKPASTJ and (DU) Slll=0.m0 DO K!!1,NlliNmES IF ( DU(KK) .LE. 1.w'25) DU(KK)30.WO SlflISUHKPASTUBDNIFLAGH).KK)"DU(KK) ENDDO RESID( I I )‘SLM IFLAGSIFLAG‘H ELSE If not, copy the value from (RINTER) and subtract the necessary coefficients (which are carried through from the solution of the system of compressed equations). RESID(I I )8RINTER(II-IFLAG) DO LL'1,NLMBDY RESID( I I )8RESID( I I )fKPASTU I , IBDY(LL))*DU( IBDY(LL)) ENDDO Finally, subtract the residual terms generated for this series of arrays from the original force vector (R) ENDIF ENDDO write (‘,100) ' in subresidual' write (*,100) ' n r(n) resid(n)' 140 Do N-I,NImNODES RITE ('2') N,R(N).RESID(N) RINI-RINI-RESIDIN) ENDDO CONTINUE C Copy whatever is left into (FINALR) for return from the subroutine. nonnn write (*,100) ' at end of sub resid. ' write (*,100) ' n r(n) finalr(n)' DO NII,NUMNODES write (*.*) N,R(N),FINALR(N) ENDDO DO NII,NUMNODES FINALR(N)8R(N) ENDDO RETURN END 141 c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA SIBRGIT INE STWEDELU(DELU,NT ,MES) INCLlDE ' KSTRE . FW ' DMLE PRECISIGI DELU(200) CHARACTER‘1S FILEDELU C This stbroutine stores the calculated values of (AU) in individasl files Min/“LS INTEGER IBDY(200) WIIBIRX/IBDY DGJBLE PRECISIN BVAL(200) COIN/MilleVAL INTEGER ISTWFLAG WIISF/ISTWFLAG DGJBLE PRECISIW DWNZDO) INTEGER Nd)EHALF,NNlNi Reinsert known bomdary values in the (DELU) array. NOTE the (DWI) vector is the 'mpscked" solution vector IATich also contains the known boundary conditions at the particular time-step. 0000 I FLAGSO DO II'1,NlHNWES IF (II.E0.IBDY(IFLAG+1)) THEN DllDiY(II)-BVAL(IFLAG+1) IFLAGIIIFLAG+1 ELSE DW(II)IDELU(II‘IFLAG) ENDIF ENDDO C Trmcate very small values to avoid mderflow errors: 00 !!=i,NlNinES IF (DABS(DINIIY(!!)) .LE. LOO-15) THEN DIMY(!!)=0.(!)0 ENDIF ENDDO IF (NT.NE.0) THEN C RITE (',100) ' STWEDELU RUN ' 100 FWMAT (A) C Nrite to appropriate array colum. 00 I81,NINAN(I)ES DELUSTRED(NT,I) I m0) ENDDO C Screen output for instant gratification. 85 FWMAT (A,I2) 142 RITE (',100)' ' RITE (',85) ' DISPLACEIENT VALLES FW TIMESTEP',NT RITE 0,100) ' VALLES LESS THAN 1.w-15 ARE STWED AS o.mo' RITE 0,100) ' ME D! DT' 00 IDIi,NlliN¢ES,2 JD-ID+1 NNlNi-(JDIZ) RITE (‘5') NM,DLOIIY(ID).DLIDIY(JD) EIDDO C This stores the (GI) vectors to disk for use with perneter C estimation software. Note that there are two format statements C since sue coapilers won't add in the leading 0 to the appended C filename call. 86 FmMAT (A,Ii) IF (ISTWFLAG .ED. 1) THEN IF (NT .LT. 10) THEN RITE (FILEDELU,86)'DELUFILE’,NT ELSE RITE (FILEDELU,‘)'DELUFILE',NT ERIF WEN(22, FILEIFILEDELU,STATUSI'NEU') DO I'LWES RITE(ZZ.')DWY( I) ENDDO ENDIF END I F c weeenummueeewm RITE(*,100)’ ”' RETURN END 143 AAAAAWAAAAWAAAAAAAAAAAAAAAAAAAAAAAAAAAAMAAAAAAAAA cAA SUBRGITINE BYVAL(MBDY,MES,MASTERK, C FINALR,IDIM,LASTR,FINALK,IDPLUSI,LSIG) DOBLE PRECISIN MASTERX(WES,WES) DWBLE PRECISIN FINALK(IDIM, IDPLUS1) INTEGER IBDY(2W) MI I B I RX] I BDY DMLE PRECI SIN BVAL(200) WIBGJNDVAL/BVAL DNBLE PRECISIN FINALR(200) DMLE PRECISIW LASTR(200) DMLE PRECISIW TER This subroutine takes the appropriate boundary conditions from the IBDY (indexed list of ROUS which have known boundary conditions) and BVAL (the values associated with the known boundary conditions) and substitute the proper values into the NASTER! array at each time step. Note that the boundary conditions are NOT time- dependent in this version. 100 FRMAT (A) nonnnn C ZERO GIT THE APPRWRIATE RWS IN [MASTERKJ AND (FINALR) 00 M ,NIAIBDII Do «I ,NUNNODES NASTERKIIBDTtI).J)-D.ODO ENDDO FINALRITBDV(I))-0.0DD ENDDO C ”TRACT THE APPRWRIATE VALUES FRM THE (FINALR) VECTW AID ZERO GIT THE C APPRWRIATE COLIMNS IN (MASTERKI DO I81,Nl.liBDY DO J81,NlliNmES FINALR(J )8FINALR(J)°MASTERK(J , IBDY( I ) )‘BVAL( I) MASTERKU, IBDY(I))'0.M0 ENDDO ENDDO NW, RESTRUCTURE THE ARRAY SO THAT THE ”ZERO“ RM AID COLLINS ARE ELIMINATED, AND THE FINAL ARRAYS OF IN] AID (R) ARE PRWERLY DI- MENSIOIED. n OOH FIRST, (R): IFLAGIO DO HILWES IF (II.EO.IBDY(IFLAG+I)) THEN IFLAGSIFLAG'PI ELSE LASTR(II-IFLAG)'FINALR(II) ENDIF EIDDO c WORWWOOM C THEN (MASTER!) : C FIRST RWS. . . 1.44 IFLAGIO DO II‘TJARANCDES IF (II.EO.IBDY(IFLAG+1)) THEN IFLAGSIFLAGH ELSE m AAA----AAAAAAAA-A---AAA c "ESIED mm ”l““'"""""""'"" -— JFLAGIO DO NILWES IF (JJ.E0.IBDY(JFLAG+1)) THEN JFLAGIJFLAGO‘I ELSE FINALKUI'IFLAG,JJ-JFLAG)IMASTERK(II,JJ) ENDIF ENDDO ceemeeeeeeeeeeemeeemeeeeeeemee END I F EIDDO C ...A)D FINALLY AUGIENT [FINALE] UITH (LASTR) FW RETRN AID SOLUTIW. Do JJ-1,IDIM FINAL!(JJ,IDPLUS1)8LASTR(JJ) ENDDO RETRN EDD 145 c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA SUBRGITINE SINLT(N,A,X,EPS,IIDIC,NRC,Y) DOUBLE PRECISION Y(200).A(N,NRC),X(200) DOUBLE PRECISION EPS,PIVOT,DETER,AIJCK,SIMUL INTEGER INDIC,NRC,N INTEGER IRONIZDO)..IOOL(200).JORD<200) C FRO! “APPLIED NUNERICAL NETNODS-, R.CARNANAN, H.A.LUTHER, C 4.0. NILKES. .I.UILEII t. SONS,NEU you, 1969 CHAPTER 5, p.276 iflfiflfifitfifitfiiitii. DIRECVORY or VARIABLES .fiiflflflttiflfififiifiifi N NUMBER OF ROUS IN A A AUGMENTED MATRIX OF COEFFICIENTS X VECTOR OF SOLUTIONS EPS MINIMUM ALLONABLE MAGNITUDE FOR A PIVOT ELEMENT INDIC COMPUTATIONAL SNITCH FOR SOLUTION TYPE (+1 FOR NO INVERSE RETURNED) NRC COLUMN DIMENSIONS FOR THE MATRIX [A] (NPI) ifififlflifififififitfltiit.tOtttfittfifittfi.ttififiitfitiOititfiiitifitfififit flflflflflflflflflfl 100 FORMAT (A) MIT NPLUSMINRC C BEGIN ELIMINATION PROCEDURE DETER81.000 DO 9 K81,N C CHECK FOR A TOO-SMALL PIVOT ELEMENT C URITE (*,*) N,K,A(!,K),EPS IF (DABS(A(K,K)).GT.EPS) GO TO 5 NRITE (*,100) ' PIVOT TOO SMALL...I OUITI' URITE (*,100) ' ERROR TRAPPED IN SUBROUTINE *SIMULT" STOP C NORMALIZE THE PIVOT ROM 5 Ian-m DO 6 J8!Pi ,NPLUSM 6 A(!,J)=A(!.J)/A(!.!) A(!,!)-1.0D0 C ELIMINATE THE K(TH) COLUMN ELEMENTS EXCEPT FOR THE PIVOT DO 9 I81,N IF (I.E0.X .OR. A(I,K).E0.0.000) GO TO 9 00 B JIKP1,NPLUSM 8 A(I.J)'A(I.J)'A(I.K)*A(K.J) A(I,K)ID.000 9 CONTINUE C URITE THE SOLUTION VECTOR INTO (X) FOR RETURN DD II'1,N X(II)'A(II,NRC) ENDDO C THAT'S ALL FOLKS... RETURN END 146 147 c AAAAAAMAAAAAAAAAAAAAAMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA SLBRGITINE SWARBYCOH ISPEC,SQIARE,COL,RESULT) C Multiplies ISOUARE] by (COL) and prochces (RESULT). C Note that there is no safety check on dimensions here. DGJBLE PRECISIN SGIAREUSPECJSPEC) DMLE PRECISIW COL(ISPEC),RESULT(ISPEC) 100 FWMAT (A) C MFRTRANIIIIIIIIHII C This stuff gets arm a sunspot factor in VMS FRTRAN CBISPEC A-COL(1) B-SOUARE(1,1) C Clean House 00 I81, ISPEC RESULT( I )80.m0 EIDDO C Do the nasty DO IRMIJSPEC DO ICOL81,ISPEC IF (DABS(SDUARE(IRW,ICOL)*COL(ICOL)) .LE. 1.m-20) GO TO 20 RESULT( IRW)IRESULT( IRW)*SOUARE( IRW, ICOL)*COL( ICOL) 20 CWTINUE E1000 E1000 RETURN END c .......................................................................... 148 Table 5. P21.FOR PROGRAM P21. FW C C PROGRAM FW 'SMART' GLWAL WTIMIZATICNII c C INITALIZE VARIABLES: C perincl.for here DOUBLE PRECISION CL(50), CU(50), FF(90), PPL(20) DOUBLE PRECISION RR(90, 30). NPEN(50), XC(30), xx<90, :0) DOUBLE PRECISION xx0LDISD) DwBLE PRECISIN ALPHAP, BETA, DELTA DNBLE PRECISIW Z, ZXC INTEGER IC, ICM, ICMI, IEV1 INTEGER IEVZ, IEV3, IIS, IOPT, IGANNA INTEGER IT, ITNAx, IZ, IZRO INTEGER szc, JC, JI, JJ INTEGER JZ, !OUNT,!!1, NALT, NC INTEGER NCNPLx, NFE, NINPS, NLEG INTEGER NRUN, NS, NSEG, NST! INTEGER NUMTIMSTEPS, NV, NPRNT WIN/ALPHAP, BETA, DELTA, IGAABIA WIBBB/IC, ICM, ICMI, IEVI COHON/CCC/IEVZ, IEV3, IIS, IWT WIDDD/IT, ITMAX, IZ, IZRO COHON/EEE/IZXC, JC, JI, JJ COHON/FFF/JZ, KGJNT,KK1, NALT, NC CONN/GGG/NCMPLX, NFE, NINPS, NLEG WIHHH/NRUN, NS, NSEG, NSTK WIIII/NV, Z, ZXC, NPRNT, NLHTIMSTEPS MUN/CL, CU, FF, PPL, RR, REN, XC, XX, XXOLD DGJBLE PRECISICNI TH CNN/THICKBLOCK/TH INTEGER IARG NFEID c .fifiitfiiiiiittfitfifiiiiitititiAtflttO C CUISTANT THICKNESS TERM: THII.m0 c *fififittifiiiifi*RtfiflfiflitfifitifiittQR.O 12 NCMPLX I 1 100 FORMAT (A) 110 FORMAT (A,IZ) 115 FORMAT (12,012.3,012.3,012.3) NRITE (*,100) ' PARAMZ ' URITE (*,100) ' Parameter Estimation Program ' 149 Table 5 (Cont'd). URITE (‘,100) ' Copyright 1991' URITE (“.100) ' Scott Morris ' URITE (*,100) ' Michigan State University ' 15000 CONTINUE 15001 IOPT I 1 NRUN I 1 ICMI I 0 NALT I 0 HST! I 0 C C NPEN(I) IS HEIGHT GIVEN TO PENALTY I IN FUNCTION 15028 KOUNT I 0 C +e++ee+e+++++++++++e+++++++++++++ C READ IN THE PARAMETRIC DATA C ee++++e++++++++++++++++++++++++++ OPEN (16,FILEI'par_est',STATUSI'Old’) RENIND 16 READ (16,*) NC,NS,NV,ITMAX,BETA,IGAMMA READ (16,*) NUMPAR READ (14,*) NUMTIMSTEPS READ (16,*) NPRNT URITE (*,100) ' Parameter indeces and limits:' URITE (*,100) ' NI CL(I) CU(I) !X(1,NI)' write 0,110) ' rquarI', nlmpsr DO III,NUMPAR write (*,110) ' II', I READ (16,*) NI,CL(NI).CU(NI),XX(1,NI) URITE (*,115) NI,CL(NI),CU(NI),XX(1,NI) END 00 c .......................................... DELTAI1.OD°A C .......................................... C NOTE THAT DELTA VALUE IS FIXED c .......................................... CLOSE (1‘) c MitiMQWAQAWQMAW CALL LINE15050 c Weeeeeumnumummememe 15200 CONTINUE IC I NC - NS 115C) IARGI1199999 do “1,10 Wranarg) enddo 15600 DO 12 I 2,NV 15610 DO 42 I 1,NS write (*.100) ' line 15610' 15420 RR(IZ, J2) I ran(iarg) write (*.100) ' line 15(20' END DO 15630 continue END 00 write (*,100) ' line 15430' 15‘50 URITE (*,100)' PARAMETER VALUES FOR THIS RUN' URITE (*.100)' NS NC NV ITMAX' 15460 URITE (*,*)NS,NC,Nv,ITNAx NRITE (*,100)' ALPHAP BETA IGANNA DELTA' 120 FORMAT (012.2,3X,012.2.17,012.2) 15470 URITE t*.120) ALPHAP, BETA, IGANNA, DELTA RITE ‘e.e) I I URITE (*,100)' Parameter Estimation Routine' NRITE (',100)' Subroutine trace:' URITE (*,100)' ' 15520 GO TO 15535 c tittififiittitifliiiiIfitfiflttfiiit c ifliiiiiififlifltiifliitttiififliti C CALL MAIN SUBROUTINE 15535 CALL LINE15700 c 0....*flifiitiiiflflfliitfliiifliflii c iifitittittitiiflifififltittflfitit 15560 IF ((IT - ITMAX) .LE. 0) THEN GO TO 15550 ELSE GO TO 15660 END IF 15550 00 J2 I 1.NS XX(IEV1.J2) ' XC(J2) END 00 15552 IIS I IEV1 c eeeeeeeeemeeeeeeeeeemeeaeemeeaeee CALL FUNC c eemeemeteeeeeeemaeeaeeeeeeeeeeeeee URITE (*.100) ' ' 151 15553 URITE (*,100) ' VALUE OF THE FUNCTION AT THE CENTROIDI' URITE (‘.') FF(IEV1) URITE (',100) ' ' 15555 NRITE (*.100) ' COODINATES OF THE CENTROID' 15556 00 J2 I 1,NS URITE (*,100) ' INDEX VALUE' 15557 URITE ('.*) J2. XC(JZ) END 00 URITE (*.100) I I 15559 URITE ('.100) I BEST NON-CENTROID VALUE OF THE FUNCTION -I URITE (*.*) FF(IEVZ) URITE (*,100) ' I 15560 URITE (',100) ' BEST NON-CENTROID X VALUES' 15590 DO JZ I 1,NS URITE (*,100) ' INDEX VALUE' 15600 URITE (*.*) JZ.XX(IEV2, J2) 15620 continue END DO 15630 RITE ('.175)' “BER OF ITERATIGISI'," 175 FORMAT (A,I3) URITE ('.175)' FUNCTION EVALUATIONSI'.NFE 15650 IF (FF(IEV1) .GE. FF(IEV2)) THEN DO J2 I 1,NS XX(1,JZ) I XC(JZ) END 00 GOTO 15691 ELSE 00 J2 I 1,NS XX(1,JZ) I XX(IEV2,JZ) END 00 GOTO 15691 END IF 15660 URITE (*,175) I THE NUNBER or ITERATIONS HAS EXCEEDED',ITMAX 15665 URITE (*,100) I PROGRAN TERMINATED PREMATURELY. Evaluations:' URITE (',*) NFE 15690 IF (IT .GE. 0) THEN IEV3 I 1 DO ICM I 2,NV IF ((FF(IEV3) ' FF(ICM)) .LE. 0.0) THEN IEV3 I ICM END IF END 00 URITE ('.100)' THE BEST FUNCTION VALUE YET IS' URITE ('.') FF(IEV3) URITE (*,100) ' Parameters associated with this best point:' DO JC I 1,NS URITE ('.100)' IEV3 JC XX( )' URITE (*.')IEV3.JC,XX(IEV3, JC) END 00 00 JJ I 1,NS XX(1.JJ) - XX(IEV3,JJ) END 00 152 END IF 15691 RITE (i'100’ I MW! URITE (*,100) ' END OF ESTIMATION RUN' “‘11; (e'1oo) o eeeeeeeeeeeeeeeeeeeee a 99999 END eeeeeeeeaeeeeeeeeeeeseeeaeeeeeeeeeeeeaaeeeeeeeee XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SUBROUTINE LINE157OO XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX maeeeeenemmemmmmee MAIN SUBRGJTINE STARTS HERE 000 on C parincl.for here DOUBLE PRECISION CL(50). CU(SO). FF(90). PPL(20) DOUBLE PRECISION RR(90. 30), NPEN(50). XC(30). xx<90, 30) DOUBLE PRECISION xx0L0(30) DOUBLE PRECISION ALPHAP, BETA, DELTA DOUBLE PRECISION Z, ZXC INTEGER IC. ICM, ICMI, IEV1 INTEGER IEvz, IEV3, IIS. IOPT, IGAIeIA INTEGER IT, max, 12, Tue INTEGER szc, JC. JI. JJ INTEGER J2. !OUNT.!!1, NALT, NC INTEGER NCNPLx, NFE, NINPS, NLEG INTEGER NRUN, Ns, NSEG, NST! INTEGER NUNTINSTEPS, NV, NPRNT COMMON/AAA/ALPHAP. BETA. DELTA. IGAMMA COMMON/BBBIIC, ICN, ICMI, IEV1 OOIDION/CCC/IEvz, IEv3, IIS, IOPT COMMON/DDD/IT, ITNAx, IZ, IZRO COMMON/EEE/IZXC, JC, JI, JJ COMMON/FFF/JZ, !OUNT,!!1, NALT, NC COMMON/GGG/NCMPLX, NFE, NINPS, NLEG COMMON/HHH/NRUN, NS, NSEG, NST! COMMON/III/NV. z, zxc, NPRNT, NUNTINSTEPS COMMON/JJJ/CL. CU, FF, PPL, RR, NPEN, xc, xx, xx0LD 100 FORMAT (A) URITE (“.100) ' (15700)' 15700 CONTINUE 15715 KSTK I 0 NFE I 0 15730 IT I 1 15760 KODE I 0 15750 IF ((NC - NS) .LE. 0) THEN GO TO 15770 ELSE GO TO 15760 END IF 15760 KODE I 1 15770 CONTINUE 15780 00 II I 2,NV 15790 00 IV I 1,NS 153 15300 xxm,m - 0.000 E1!) 00 ETD 00 15820 00 II I 2,NV 15830 00 IV I 1,NS 15860 IIS I II 15850 CONTINUE 15860 XX(II. IV) I CL(IV) + RR(II, IV) * (CU(IV) - CL(IV)) 15870 continue END 00 15880 KK1 I II 15890 CONTINUE c memmeeeaaemmeeemee C CHECK CONSTRAINTS 15895 CALL LINE167OO CMWWW 15900 IF ((II - 2) .LE. 0) THEN GO TO 15910 ELSE GO TO 15970 END IF 15910 IF (IPRINT .NE. 0) THEN GO TO 15920 ELSE GO TO 15990 END IF 15920 URITE (*,100) 'COORDINATES OF INITIAL COMPLEX' 15960 ID I 1 15965 DO IV I 1,NS 15950 URITE (“.100) XX(IO, IV) 15955 continue END DO 15970 IF (IPRINT .NE. 0) THEN GO TO 15975 ELSE GO TO 15990 END IF URITE (',100) ' II IV XX(II,IV)' 15975 00 IV I 1,NS 15980 URITE (‘.*) II, IV, RR(II, IV) 15985 continue END 00 15990 continue END 00 16000 !!1 I NV 16010 00 IIS I 1,NV 16020 CONTINUE 16025 NFE I NFE + 1 C eeeeeeeeeeeeeeaeeeeeeeeeeases CALL FUNC c eeeeeeeeeeeeeeeeeeeeeeeeeeeee 154 16030 continue END DO 16060 KDUNT I 1 16050 IA I 0 16060 IF (IPRINT .NE. 0) THEN 00 TO 16070 ELSE 00 TO 16110 END IF 16070 HRIIE (*,100) 'VALUES OF THE FUNCIION' URITE (',100) ' IV FF(IV)’ 16080 00 IV I 1,NV 16090 URITE (*.*) IV,FF(IV) 16100 contInuo END DO 16110 IEV1 I 1 16120 DO ICN I 2,NV 16130 IF ((FF(IEV1) ' FF(ICM)) .LE. 0.000) THEN GO TO 16150 ELSE GO TO 16140 END IF 16140 IEV1 I ICN 16150 continu- END 00 16160 IEV2 I 1 16170 D0 ICN I 2,NV 16180 IF ((FFIIEVZ) ' FF(ICN)) .LE. 0.0) THEN 00 TO 16190 ELSE 68 TO 16200 END IF 16190 IEV2 I ICN 16200 continue END 00 16210 IF ((FF(IEV2)'(FF(IEV1)+8ETA)) .LT. 0.0) THEN 00 TO 16240 ELSE GO TO 16220 END IF 16220 KOUNT I 1 16230 GDTO 16260 16250 KOUNT 8 KOUNT I I 16250 IF ((KOUNT ° IGANHA) .LT. 0) THEN GO TO 16260 ELSE GO TO 16600 END IF 16260 CONTINUE C Qtittittfittiititttitttittflirt! C CDNPUTE CENTHOID 16265 CALL LINE17000 11555 c Ofififitfififififfii.i***..fiflfitfiiitfitfi 16267 IF (IT .LE. (2 * NV)) THEN ALPHA I 1.6 ELSEIF (((2 * NV) .LT. IT) .AND. (2 ‘ NV) .LE. (4 * NV)) THEN ALPHA I 1.3 ELSEIF ((4 ‘ NV) .LT. IT) THEN ALPHA I 1 END IF 16268 IF (KOUNT .GT. 0) THEN ALPHA I .8 16269 IF (IALPH .E0. 0) THEN ALPHA I ALPHAP 16271 IF (NSTK .LT. 1) 00 TO 16275 16272 00 JJ . 1,us 16273 XXOLDIJJ) I XX(IEV1, JJ) END DO 16274 KSIKZ I 0 16275 00 J.) I 1,NS 16280 XX(IEV1, JJ) I (1.0D0 + ALPHA) * (XC(JJ)) - ALPHA * (XX(IEV1, JJ)) 16285 continue END DO 16290 IIS I IEV1 16300 CGITINUE c cuscx CONSTRAINTS******************* 16305 CALL LINE16700 c .O...‘*C***ififitifi*fitttitittiifiifltflflfi 16310 CONTINUE 16315 NFE I NFE + 1 C COMPUIE FUflCTlouttttttttttfittttfitfitfi CALL FUNC C tittiittfittttitttitttttiifittiitittt. 16320 IEV2 I 1 16330 00 1c» . 2,uv 16340 IF ((FF(IEV2) - FF(ICM)) .LE. 0.0) IHEN 00 10 16360 ELSE 00 10 16350 END It 16350 IEV2 I ICN 16360 continue END DO 16370 IF ((IEV2 - IEV1) .NE. 0) THEN GO TO 16450 ELSE GO TO 16380 END IF 16380 KSTK I KSTK + 1 110 FORMAT (A,12) IF (KSTK .GE. 8) THEN URITE (I,100) ' HELP! 1~u STUCK! -ERROR TRAP AT 16380-' URITE 1*,100) ' xsrx >6' 156 STOP ENDIF NSTK I NSTK I 1 16381 IF ((NSTK .LT. 3) .AND. (KSTK .GE. 6)) THEN URITE (“.100) ' JUNPINO OUT AT LINE 16381' GO TO 16600 END IF 16382 IF ((NSTK .GE. 3) .AND. (KSTK .GE. 6) .AND. (KOUNT .LT. 2)) THEN KSTKZ I KSTKZ I 1 IF (KSTK2 .GE. 2) THEN URITE (*,100)' JUNPINC OUT AT LINE 16382' GO TO 16600 END IF IEV3 I 1 DO ICM I 2,NV IF ((FF(IEV3) - FF(ICN)) .LE. 0.0) THEN IEV3 I ICM END IF END D0 D0 JJ I 1,NS XC(JJ) I XX(IEV3, JJ) URITE (*,*) XC(JJ) XX(IEV1, JJ) I XXOLD(JJ) END 00 KSTK I 0 URITE (',100) ' REPLACING CENTROID DY BEST POINT' GOTO 16275 END IF 16385 00 JJ I 1,NS 16390 XX data file: NUMELS NUMNODES C1C2C3¢405¢6 NUMBEROFTIMESTEPS NUMSTART NUMSTOP TIMSTART TIMINCR IIII For the Data File: 1 wafi) NGWU) < ----- These are the x,y coordinates for the 2 NODX(2) NODY(2) nodes. n NmXIn) NcDfln) 1 EI(1) EJ(1) EK(1) EHII) <--- These are the nodes INIich mite q) 2 EI(2) EJ(2) EI((Z) EH(2) ' each elenent. For a triangular . . . . . element, the EM ) value should be . . . . . 0 (Integer). n EI(n) EJ(n) EKIn) EN(n) IIII For the <80u1daries.dat> data file: MBDY IBNDX IDIR BVAL(IBNDX) \ IBNDX IDIR BVAL(IBNDX) 1"” NWT WI“ of times IBNDX IDIR BVAL(IBNDX) / IBNDX The node at which the bomdary condition is know) >>>>>>>THESE HIST BE IN SEOUENTIAL MDERlll<<<<< IDIR Direction of Displacement: XI1 YI0 BVAL(IBNDX) Displacement Value Note that the index systee obviates the need for semential ordering of the bomday conditions EXCEPT that d). to a bin in the program, the "1" IDIR "IIISTII precede the '0' IDIR for any node INIiCh has two known bomdary conditions. IIII for the file: II Number )1 I 1 and Y I 0 NulRVAL 171 Table 8 (Cont'd). INDEX IDIR VAL INDEX Node nulber IDIR Direction code (see above) VAL Force value 172 The F118 PAR_EST.DAT required for using P21.F is: NO NS NV ITMAX BETA IGAMMA NUMPAR NUNTINSTEPS NPRNT NIl CL(NIl) CU(NIl) XX(1,NIl) N12 CL(N12) CU(N12) XX(1,NI2) Where NPRNT Print level during function evaluation 0 - Prints nothing after first evaluation 1 - Prints new coefficients after first run 2 - (1) plus displacement values NUMPAR Number of parameters to be estimated NC Number of Constraints NS Number of Search Variables NV Number of Vertices ITMAX Maximum Number of Iterations BETA Model Parameter (Convergence Criteria) IGAMMA Model Parameter NIn Numerical Index for Parameter n CL(NIn) Lower Constraint for Parameter n CU(NIn) Upper Constraint for Parameter n XX(1,NIn) Best Guesstimate of Parameter n 173 The File Elastic.DAT required for usingg p21a.F is: NV ITMAX BETA IGAMMA NPRNT N11 CL(NIl) CU(NIl) XX(1,NIl) Where NPRNT Print level during function evaluation 0 - Prints nothing after first evaluation 1 - Prints new coefficients after first run 2 - (1) plus displacement values NV Number of Vertices ITMAX Maximum Number of Iterations BETA Model Parameter (Convergence Criteria) IGAMMA Model Parameter NIn Numerical Index for Parameter n (n31 for all elastic data) CL(NIn) Lower Constraint for Parameter n CU(NIn) Upper Constraint for Parameter n XX(1,NIn) Best Guesstimate of Parameter n APPENDIX B Material Deflection Data 174 Note that the symmetry of the sample was exploited to simplify the parameter estimation process. The 1/4 grid (Figure 31) deflection points are listed below. 175 (6) (5) (4) (3) (2) 1 3 (1) 2 4 Figure 30. Finite Eleoent Grid Used in Parameter Estimation Braluation Model . 10 176 Table 9. Change in Position of Irdeces at 100% loading. Tinestep Point xy: 0 1 2 3 4 S 6 1 x 2 0 0.0656 0.1517 0.1662 0.2667 0,3173 0.3563 1 y i 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 2 x : 0 0.0662 0.1243 0.2197 0.2629 0.3129 0.3563 2 y l 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 3 x ' 0 0.0796 0.0606 0.0967 0.1435 0.1695 0.1613 3 y 0 0.0409 0.0944 0.0947 0.1159 0.1432 0.1375 4 x 0 0.0437 0.0799 0.1113 0.1631 0.1636 0.1613 4 y 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 5 x 0 0.0596 0.1020 0.1267 0.1637 0.1661 0.2013 5 y I 0 0.0569 0.1100 0.1393 0.1919 0.2452 0.2516 6 x i 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 6 y i 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 7 x 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 7 y 0 0.0701 0.1521 0.2013 0.2791 0.3210 0.3563 6 x 0 0.0243 0.0675 0.0902 0.1260 0.1167 0.1375 6 y . 0 0.0570 0.0776 0.1167 0.1606 0.1959 0.1613 9 x 2 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 9 y : 0 0.0952 0.1556 0.2273 0.2412 0.3106 0.3563 10 x ' 0 0.0000 0.0000 0.0000 0.0000 0.0000 0 10 y 0 0.0730 0.0750 0.1033 0.1592 0.1917 0.1613 Position of Irxiex Mark at 100% Loading. Values in Inches. pp. oomomoqqmmmmobwwMNuw K x‘