# Full text of "DTIC ADB015343: HELP - A Multimaterial Eulerian Program in Two Space Dimensions and Time"

## See other formats

UNCLASSIFIED AD NUMBER ADB015343 LIMITATION CHANGES TO: Approved for public release; distribution is unlimited. FROM: Distribution authorized to U.S. Gov't, agencies only; Test and Evaluation; APR 1976. Other requests shall be referred to Air Force Armament Lab.^ Eglin AFB^ FL. AUTHORITY ADTC Itr 3 Jul 1979 THIS PAGE IS UNCLASSIFIED THIS REPORT HAS BEEN DELIMaED AND CLEARED FOR PUBLIC RELEASE UNDER OOD DIRECTIVE 5200.20 AND NO RESTRICTIONS ARE IKPOSED UPON ITS USE AND DISCLOSURE. DISTRIBUTION STATEMENT A APPROVEO FOR PUBLIC RELEASE; DISTRIBUTION UNLIHITED. ^ELP • 4 ^ULTIMATERIAL JULERI AN JROGRAM InTwO IPACE (DIMENSIONS AND JIME, j ' SYSTEMS, SCIENCE AND SOFTWARE F. 0. BOX 1620 LA JOLLA, CAUFORNIA 92037 Laura. “S Vnrt\ ff, KlpeYt T Distribution limited to U. S. Government agencies only; this report documents test and evaluation; distribution limitation applied April 1976 . Other requests for this document must be referred to the Air Force Armamenc Florida 32542. Laboratory (DUW) , Eglin Air Force Base AIR FORCE ARMAMENT LABORATORY '-Udlhlkbdk: r. ■ .. v'-lU' .-y^iH.-/.^ •: t W w«ila>tew<iifcita.»o ...i..w,v.t; .•■ ■■ ^i »W security classification of this page (tWi«n Palm Enlttud ) READ INSTRUCTIONS BEFORE COMPLETING FORM REPORT DOCUMENTATION PAGE . report NUMBER 2 GOVT ACCESSION NO. 3 RECIPIENT'S CATALOG NUMBER AFATL-TR-76-45 4 . title I'and Submit ) HELP - A MULTIMATERIAL EULERIAN PROGRAM IN TWO SPACE DIMENSIONS AND TIME, BOOKS 1 AND 2 S. TYPE OF REPORT A PERIOD COVERED Final Report July 1974 to June 1975 6 PERFORMING ORG. REPORT NUMBER 7. AuTMORfu) Laura J. Hageman Jerry L. Waddell Darin E. Wilkins Robert T. Sedgwick 9. PERrORMING ORGANIZATION NAME AND ADDRESS Systems, .Science and Software •/ P 0 Box 1620 i, ^ La Jolla, California 92037 II. CONTROLLING OFFICE NAME AND ADDRESS Air Force Armament Laboratory (DUW) Armament Development and Test Center Eglin Air Force Base, Florida 32542 8. CONTRACT OR GRANT NUMBERCaJ F08635-75-C-0044 / 10. PROGRAM element, PROJECT, TASK AREA a WORK UNIT NUMBERS 25600232 12 REPORT DATE April 1976 13. NUMBER OF PAGES 437 14 MONITORING AGENCY NAME » bOOfltSS(ll dllletml Irom Cnnirnlllne Ollirt) IS. SEC U Rl T Y CL ASS. (III* ri-porfj UNCLASSIFIED 'fSA“TE.CL~ASSiFlC ATION ' DOWNGRADING schedule 16. DISTRIBUTION ST ATEMEN T fol (h(» R»porO Distribution limited to U. S. Government agencies only; this report documents test and evaluation; distribution limitation applied April 1976. Other requests for this document must be referred to the Air Force Armament Laboratory (DUW), Eglin Air Force Base, Florida 32542. I *7. distribution statement (ot fh# 9bHtrmct •nferetl In Btfyrh 20, It dlttBrent from B 0 port) 18. supplementary notes Available in DDC. 19. KEY WORDS (Confinue on revorao aide it n9C0B9mry anti Idantifv by block nowbar) HELP Multimaterial Interaction Elastic-Plastic Massless Tracer Particles Continuum Mechanics Hydrodynamic Codes RPM, OIL, and PIC ABSTRACT (Continua on ravarga aide If neceeeary and identity hv block niimhar) HELP is a two-dimensional, multimaterial, Eulerian Code for solving material flow problems in the hydrodynamic and elastic-plastic regimes. Although the code is basically Eulerian, material interfaces and free surfaces are propagate through the Calculation Mesh, in a Lagrangian manner, as discrete interfaces across which the materials are not allowed to interdiffuse. This interface treatment gives HELP a distinct advantage over other pure Eulerian Codes and allows it to be applied to the solution of a variety of complex, multimaterial problems. The basic HELP techniques have been exercised over a wide range of^ W 1 JAN 73 EDITION OF I NOV 65 IS OBSOLETE UNCLASSIFIED l SECURITY CLASSIFICATION OF THIS PAGE (Whtn Data Enitttd) L 1 UNCLASSIFIED 8ICURITY CLAI»iriCATION Of THIS f AOSf)t7i«n Dal* Bnlan^ ITEM 20 (CONCLUDED): solutions including hypervelocity and ballistic impact, fragmentation munitions, shaped charge jet formation, and shaped charge jet penetration. Considerable confidence has been gained in the basic numerical model . UNCLASSIFIED SPCURJTY classification of this PA0EflW?«n Bnffd) msmB^mansssEsa PREFACE This report documents work performed during the period July 1974 through June 1975 by Systems, Science and Software, P. 0. Box 1620, La Jolla, California 92036, under contract F08635-75-C-0044 with the Air Force Armament Laboratory, Armament Development and Test Center, Eglin Air Force Base, Florida 32542. Mr Leonard L. Wilson (DLDG) managed the pregram for the Armament Laboratory. This technical report has been reviewed and is approved for publication. FOR THE C0M4ANDER: GERALD P. D'ARCY, Colonel, USAF, Chief, Guns, Rockets 6 l2xplosiv< Division i .J- (The reverse of this page is blank) CONTENTS Chapter Page INTRODUCTION GENERAL DESCRIPTION OF THE NUMERICAL METHOD 2.1 CONSERVATION EQUATIONS 2.2 DIVISION INTO PHASES HI 2.2.1 SPHASE - The Effects of Deviatoric Stresses . . . , 2. 2. 1.1 Continuity Equation 2. 2. 1.2 Equation of Motion , 2. 2. 1.3 Energy Equation 2.2.2 HPHASE - The Effects of Pressure 2. 2. 2.1 Continuity Equation 2. 2. 2. 2 Equation of Motion , 2-11 2-11 2-11 2. 2. 2. 3 Energy Equation 2. 2. 2. 4 Treatment of Free Surfaces and Large Density Discontinuities 2-12 2-13 2, 2. 2. 5 Artificial Viscosity 2,2.3 TPHASE - The Effects of Transport . 2-16 2-17 2. 2. 3.1 Continuity Equation 2-18 2. 2. 3. 2 Equation of Motion 2. 2. 3.3 Energy Equation 2.3 MATERIAL MODEL 2-20 2-22 2-22 2.3.1 Equation of State 2-23 2.3. 1.1 High Explosives 2-23 CONTENTS Chapter 2. 3. 1.2 Inert Materials 2. 3. 1.3 Ideal Gas Page 2-26 2-27 2.3.2 Elastic-Plastic Constitutive Relation. 2-29 2. 3. 2.1 Strain Rate Deviators 2-29 i 2.3. 2.2 Stress Deviators 2. 3. 2. 3 Yield Criterion , 2.3.3 Material Failure in Tension 2.3.4, Detonation of High Explosives 2. 3. 4.1 Detonation Procedure 2. 3. 4. 2 Definition of the Detonation Front . . . , 2-32 2-33 2-34 2-34 2-34 2-39 GRID BOUNDARY CONDITIONS 3.1 BASIC ASSUMPTIONS 3.2 STRENGTH PHASE (SPHASE) 3.2.1 Definition of Strain Rate Derivatives for Cells at a Grid Boundary 3.2.2 Definition of Interpolated Strain Rates and Stresses for Cells at a Grid Boundary 3.2.3 Definition of Velocities and Deviator Stresses at Grid Boundaries 3.2.4 Correction to the Theoretical Energy for Work Done at Grid Boundaries in SPHASE 3.3 HYDRODYNAMIC PHASE (HPHASE) 3.3.1 Definition of Velocities and Pressures at Transmittive Grid Boundaries CONTENTS 3.3.2 Definition of Velocities and Pressures at Reflective Grid Boundaries 3.3.3 Correction to Theoretical Energy for Work Done at Grid Boundaries in HPHASE 3.4 TRANSPORT PHASE (TPHASE) 3.4.1 Transport of Masse, Momentun and Energy Across Transmittive Grid Boundaries 3.5 TRACER PARTICLE MOTION NEAR GRID BOUNDARIES 3.6 EDITING OF FLUXES AT GRID BOUNDARIES IV MATERIAL INTERFACES AND INTERFACE CELLS . 4.1 DEFINITION OF MATERIAL INTERFACES .. 4.1.1 Use of Tracer Particles 4.1.2 Movement of Tracer Particles 4.2 CREATION OF INTERFACE CELLS 4.2.1 Use of MFLAG Array 4.2.2 Definition of Material Arrays 4.2.3 Accounting for Subcycles ..... 4.3 CALCULATION OF CELL FACE AREAS FOR TRANSPORT ACROSS INTERFACE CELL BOUNDARIES 4.3.1 Defining Fractional Areas of Intersected Cell Faces ...... Page 3.4.2 Change of Momentum for Cells at Reflective Grid Boundaries in TPHASE.. 3.4.3 Correction to Theoretical Energy for Energy Transported Across Grid Boundaries in TPHASE 3-9 3-10 3-12 3-12 3- 12 4- 1 4-1 4-1 4-2 4-4 4-4 4-6 4-8 1 4-9 1 1 : o t ■ M CONTENTS Chapter Page 4.3.2 Presetting Non-Intersected Cell Face Areas 4-11 4.3.3 Resetting Non-Intersected Cell Face Areas 4-15 4.4 ADJUSTMENT OF INTERFACE CELL MASS TRANSPORT TERMS 4-17 4.4.1 Evacuation of a Material 4-18 4.4.2 Overemptying of a Material 4-19 4.5 TYPES OF INTERFACE CELLS 4-20 4.5.1 Multimaterial Cells 4-20 4.5.2 Free Surface Cells 4-20 4.5.3 Slipline Cells 4-21 4.6 DETERMINATION OF MATERIAL PROPERTIES IN MULTIMATERIAL CELLS 4-22 4.6.1 Pressure and Density Iteration for Multimaterial Cells 4-22 4. 6. 1.1 Pre- iteration Calculation .... 4-23 4.6. 1.2 Iteration Calculation 4-25 4. 6. 1.3 Equation of State Modifications for the Iteration Procedure 4-28 4.6.2 Shear Yield Strength of Multimaterral Cells 4-29 4.6.3 Partitioning of Internal Energy Increments in Multimaterial Cells 4-30 V SLIPLINES 5-1 5.1 GENERAL COMMENTS 5-1 5.2 RELATIVE ACCELERATION OF MASTER AND SLAVE MATERIALS (UVCALC) 5-2 vi BtSSwEg!g3!TO'‘'IBWjatl'lJ-lil 1 1 III im iinw — I ■■ ■-v»?r*nv .-^f, -«|A ' CONTENTS Chapter Page 5.3 TRANSPORT OP MASTER AND SLAVE MATERIALS 5-3 5.4 MODIFICATION OF POST-TPHASE MASTER AND SLAVE VELOCITIES (UVMOD) 5-4 5.5 USE OF SLIPLINES 5-5 VI PLUGGING FAILURE 6-1 6.1 GENERATION OF A PLUGGING CALCULATION 6-1 6.1.1 Designation of Material Packages 6-1 6.1.2 Placement of Plug Tracer Particles .... 6-1 6.1.3 Placement of Passive Tracer Particles.. 6-3 6.1.4 Slipline Specification 6-3 6.2 METHOD OF EVOLVING THE PLUG 6-4 6.2.1 Plugging Failure Criterion 6-4 6.2.2 Conventions Governing Plug Tracer Particles 6-7 6. 2. 2.1 Partially Formed Plug 6-7 6. 2. 2. 2 Completely Formed Plug 6-8 6.2.3 Movement of Plug Tracer Particles 6-10 6.2.4 Conversion of Target Material into Plug Material 6-12 6. 2. 4.1 Redefinition of Tracer Particles 6-12 6. 2. 4. 2 Redefinition of Target and Plug Masses 6-15 VII DESCRIPTION OF INPUT VARIABLES AND RESTART PROCEDURES 7-1 7.1 GENERAL CAPABILITIES AND LIMITATIONS 7-1 7.2 INSTRUCTIONS FOR GENERATING A PROBLEM 7-3 m vii Chapter CONTENTS Page 7.2.1 Z'Block Variables 7-3 7.2.2 Cell Dimensions 7-14 7 .1.1> Initial Density, Velocity, and Specific Internal Energy of Each Material Package 7-14 1 .1 A Strength Constants for Each Material Package 7-16 7.2.5 Material Tracer Particles 7-16 7.2.6 Slipline Endpoints 7-18 7.2.7 Definition of High Explosive Detonation Points 7-19 7.3 RESTART PROCEDURES 7-20 7.4 REDEFINING TRACER PARTICLES WHEN RESTARTING A CALCULATION 7-23 VIII INTERACTIVE CAPABILITIES 8-1 8.1 REZONING THE GRID 8-1 8.1.1 Combining Cells 8-1 8.1.2 Adding Material 8-2 8.1.3 Activating the Rezone 8-3 8. 1.3.1 Automatic Rezone 8-4 8. 1.3. 2 User-activated Rezone 8-4 8.2 AUTOMATIC ADDITION OF MATERIAL TRACER PARTICLES 8-5 8.3 REDEFINING SLIPLINES 8-5 8.4 CLOSING A VOID 8-6 8.4.1 Identifying Void Closing Region 8-7 8.4.2 Criteria for Closing a Void 8-7 viii CONTENTS Chapter Page 8.4.3 Method for Closing the Void 8-7 IX ERROR CONDITIONS 9-1 9.1 ALPHABETIC LIST OF ERROR MESSAGES 9-1 9.2 INTER DIAGNOSTIC PRINTS 9-6 X GUIDE TO THE FORTRAN PROGRAM 10-1 10.1 A GENERAL FLOW DIAGRAM 10-1 10.2 DESCRIPTION OF THE SUBROUTINES 10-1 10.3 DICTIONARY OF IMPORTANT VARIABLES 10-20 10.4 FLAGS AND CONVENTIONS 10-93 10.4.1 Flags Governing Interface Cells .... 10-93 10.4.2 Definition of Transport Variables .. 10-95 10.4.3 Radial and Axial Terras 10-97 XI SAMPLE HELP INPUT AND OUTPUT 11-1 11.1 SHAPED CHARGE LINER COLLAPSE 11-1 11.1.1 Input 11-7 11.1.2 Cycle 0 Output 11-11 11.2 PERFORATION OF A THIN TARGET 11-38 11.2.1 Input 11-43 11 2.2 Cycle 0 Output 11-47 11.3 IMPACT INTO A THICK TARGET 11-70 11.3.1 Input 11-73 11.3.2 Cycle 0 Output 11-77 XII REPRESENTATIVE APPLICATIONS OF THE HELP CODE 12-1 12.1 HYPERVELOCITY IMPACT 12-2 IX ^ CONTENTS /hapter Page 12.2 PLUGGING FAILURE 12-10 12.3 LONG ROD PENETRATION AT OBLIQUE INCIDENCE... 12-16 12.4 MULTIPLE IMPACT 12-18 12.5 FRAGMENTATION MUNITIONS 12-22 12.6 SHAPED CHARGE JET FORMATION 12-27 12.7 SHAPED CHARGE JET PENETRATION 12-32 12.8 WAVE PROPAGATION STUDY 12-36 REFERENCES 12-41 APPENDIX A - Dimensioning the Arrays A-1 APPENDIX B - Abbreviated HELP Input Instructions B-1 APPENDIX C - Segmenting the HELP Code C-1 CHAPTER I INTRODUCTION This is the second complete documentation of the HELP^^^ code since it was developed five years ago. For those un- familiar with HELP, the following remarks will give some historical perspective. The HELP code evolved from three major hydrodynamic codes developed over the past twenty years: RPM, OIL, and PIC. Starting from the one-material code, RPM^^^, the HELP code added the capability of describing multimaterial interactions and of treating material strength as elastic-plastic rather than rigid-plastic. RPM, in turn, was an extension of OIL^^^, a compressible fluid code, and included material strength in the rigid-perfectly plastic approximation. The OIL code was derived from the PIC^^^ code and replaced the discrete PIC particles by a continuum. This latter transition to a continuum calculation was made in order to permit the treatment of very small compressions without having to use an extremely large number of particles for ade- quate mass resolution. The disadvantage of this transition, however, was the loss of certain Lagrangian- type features of PIC, and specifically the ability to treat flows containing more than one mate'll al. The m; ;or purpose of developing the HELP code was to restore che multimaterial capability of PiC while retaining the ease of treating small compressions and the computing economy of a continuous Eulerian model. The accurate treatment of material interface or free surface motion, however, requires a detailed knowledge of surface location, which is not contained in a purely Eulerian model. The developers of the HELP code, therefore, deviated from the pure Eulerian model by employing massless tracer particles to track material interfaces and free surfaces. 1-1 These particles are initially placed along the surface bound- aries and are moved thereafter with a local material velocity, thus creating a Lagrangian-type definition of moving surfaces without sacrificing the capability of easily treating extreme distortions. Furthermore, diffusion of materials across these surfaces is prevented in HELP by taking into account the surface positions when transporting material across the Eulerian cell boundaries. In the years since HELP was first developed and docu- mented it has been applied to a wide variety of multimaterial problems in continuum mechanics. In the process it has been significantly expanded to increase its generality and applic- ability. For example, sliplines were added in order to model more accurately the plugging failure in thin targets^^’^^ r 7 1 and the high-explosive metal interactions in shaped charges ^ and fragmenting munitions . Furthermore, a high explosive detonation model was added which would account for multiple explosives, multiple detonation points and wave shapers. As the applications of the HELP code increased the need for a more general problem generator grew. The original problem generator was written primarily for calculations involving spheres and cylinders impacting finite and semi-infinite plates. This original generator clearly was inadequate for generating shaped charge problems or even projectiles with various shaped noses. The new HELP gener- ator can, without modification, generate any configuration where the package boundaries can be decomposed into straight line segments or arcs of circles and ellipses. In addition to these improvements, some aspects of the code have been reformulated. For example, the hydrodynamic phase (HPHASE) has been simplified (the equations are solved in a single pass) in order to provide more accurate definitions of pressures and velocities near free surfaces. (See Section 2.2.2.) This has eliminated the need 1-2 to "glue" free surface cells, a method v/hich proved unsatis- factory in certain situations. The single pass formulation also simplified the addition of an artificial viscosity term to cell boundary pressures. Another change involved the ordering of the three phases. It has been found that more realistic strain rates are computed using the post - transport cell velocities. Therefore the strength phase now precedes the hydrodynamic phase in the computational cycle. This reordering of the phases also led to their renaming: PH3 has become SPHASE; PHI has become HPIIASE, and PH2 has become TPHASE. And finally, the INFACE subroutine has been broken down into five subroutines, thereby reducing the core memory requirements of the code when it is segmented, and, perhaps more importantly, clarifying the use and treatment of the material interfaces. A summary of the major chapters of the text follows. Chapter II gives a general description of the numerical method and the material models employed by HELP. The grid boundary conditions are discussed in detail in Chapter III as they apply to each phase of the calculation. Chapter IV offers a discussion of the material interfaces, as well as the determination of the pressure and shear yield strength for multimaterial cells. The equations which govern the slipline cells are given in Chapter V along with general guidelines on their use, and in Chapter VI the mechanisms added to the code for modelling plugging failure are described. Procedures for generating and restarting a HELP calculation are given in Chapter VII as well as some general statements about the code's capabilities and limitations. The inter- active options available to the user are described in Chapter VIII, and a list of the error conditions are discussed in Chapter IX. In Chapter X a general flow diagram of the code, a description of all the subroutines, and a dictionary of variables is given. Several sample input decks for the 1-3 various types of problems to which HELP is commonly applied are given in Chapter XI. The HELP output is also described in this chapter. The last chapter, Xll, describes several applications of the HELP code. The three appendices give instructions for dimensioning the arrays, an input form with abbreviated instructions for generating initial conditions, and a diagram indicating how the code might be segmented to save core storage. This second documentation of the HELP code is consider- ably more detailed than the original and is more oriented toward the new user of the code. As with any large, complex computer code, HELP is most effectively used by those persons who have a good grasp of its numerical methods, material models, conventions, assumptions, capabilities and limitations. It is hoped by the present authors that this document will provide this basic information. Any suggestions for improve- ments in the code or its documentation will be appreciated by the authors. 1-4 CHAPTER II GENERAL DESCRIPTION OF THE NUMERICAL METHOD In the present section the conservation equations and the finite difference analog to treat these equations are given. The material model, which includes strength effects, is also presented. I 2.1 CONSERVATION EQUATIONS Space is divided into fixed Eulerian cells through which the material moves. To arrive at expressions for the rate of change of total mass, momentum and energy within such a cell, it is convenient to start with the conservation equations governing the motions and interactions of contin- uous media in the form: f?-' ( H (Continuity Equation) ^ (pu.) 3t dx7 (2.1) 1 Du. - ! f (Equation of Motion) P Tt " -577 (2.2) II r > (Energy Equation) DEt g ^ UT " ^ • (2.3) Here is the stress tensor, which can be regarded as the sum of the hydrostatic stress, and a stress deviator tensor, s^j , i . e. , . . = s. . - 6. .P 11 11 11 and ^i^i ^i total energy (kinetic plus internal) per unit mass. Tensor notation is implied, so that repeated indices denote summations. 2-1 Expanding the convective derivatives in Eqs. (2.2) and (2.3), Df/Dt = 9f/3t + u. 3f/3x., then adding Eq. (2.1) times u. to Eq. (2.2), and F.q. (2.1) times to Eq. (2.3), and collecting terms, gives ^ ^ “ij ■ -sIt (2.5) ( 2 . 6 ) For the developments to follow it is desirable to re- place these differential equations by the analogous integral equations, obtained by integrating over the cell volume, V, and then converting the volume integral of divergences to surface integrals over the cell surfaces. Equations (2.1), (2.5), and (2.6) then become Jt f pdv - - / pUj^n^dS (2.7) 3 Jt f pu.dV = J Oijn.dS J puju.njds (2.8) 3 f pE.j,dV =J o^jU^n^dS -J pu^E.j.n^dS . (2.9) 2-2 li II 2.2 DIVISION INTO PHASLS It is convenient to express the integral conservation relations, Eqs. (2.7) through (2.9), as finite difference equations over the time step At and also to decompose the total stress, o^j, into its deviator and hydrostatic components, according to Eq . (2.4). This gives, for the increments of total mass (m) , momenta (muj) and energy (mE.p) within the cell, Am = A(mUj) A(mE.p) tf Sj.njdS ■'s tf Sj.u.njdS PnjdS Pu^n^dS t / pu.n.dS •'s (2 tf (pu.u.)nj S (2 tj (pu^E.p)n^ ( 2 . 10 ) ( 2 . 11 ) ( 2 . 12 ) Here, the terms on the right are divided into increments due to stress deviator forces on the cell surface (first column), those due to tlie pressure forces on the cell surface (second column), and those due to the transport of mass, momentum and energy through the surface of the cell (third column). These three types of increments are accounted for in distinct phases of the computation. Specifically, during each time step all cells are updated three times to account for the following: • Effects of the stress deviators (SPHASE, Strength Phase)^ • Effects of pressure (HPHASE, Hydrodynamic Phase)^ •Effects of transport (TPHASE, Transport Phase). The phases are calculated in the order listed. In the dis- cussion that follows, the calculation of the terms on the right of Eqs. (2.10) through (2.12) are described sequentially starting with SPHASE. The equations are written in cylindrical * See DT in Section 10.3. coordinates although tlie code models both the axisynunetric and plane strain cases. Some preliminary definitions will be useful. Super- script n on a variable refers to the value of the variable at the beginning of the time step and superscript (n+1) de- notes the value at the end. In this discussion a typical cell in the interior of the grid is considered; Chapter III has a discussion of the special conditions which describe the calculation at the grid boundaries or at the axis of symmetry. For a typical cell, denoted by a value of the index k, the dependent variables for that cell are written P(k), u(k), v(k) , Ej(k), m(k), representing respectively the pressure, radial and axial components of velocity, the specific internal energy and the mass for cell k. The adjacent cells above, below, to the right and left of k will be designated respectively as ka, kb, kr and kn . Here the terms above, below, right and left refer to a cross- section view of the cells, in which the left border is parallel to the axis of symmetry with z increasing upward (see Figure 2.1). Each cell is, then, the torus obtained by rotating the rectangle (since Ar f Az in general) about the axis of symmetry. 2.2.1 SPHASE - The Effects of Deviatoric Stresses This section discusses the modification of the cell velocity and internal energy due to the stress deviators. The derivation of the stress deviators is given in Section 2.3.2. The SPHASE subroutine should be executed only if at least one material in the problem has shear yield strength. 2. 2. 1.1 Continuity Equation. (2.10) No contribution in SPHASE. 2. 2. 1.2 Equation of Motion, (2.11' The SPHASn contribution to the equation of motion (2.11) f Ag(muj) = At f s. n 7 1) a dS - The axial and radial momentum increments will be discussed separately. Effect on Axial Motion - In this case the velocity Uj is the axial velocity, and the term s^^n^ is the axial component of the deviatoric stress on a surface. For a torus with rectang- ular section, s^jn^ on the various faces is s^^^ cell top - s^^^ cell bottom s^,-' outer cell surface rz - inner cell surface where zz and rz subscripts denote normal and shear stresses. The top, bottom, right and left surfaces are indicated by (t) » (b) , (r), and ()i), respectively. The equation for axial motion is therefore Ag(mv) V( t\ - r^lAt + s^^^ 27rr2 Az At - s^*'^ 27rri Az At where rj^ and r 2 are the inner and outer radii, respectively, of the cell and Az is its axial dimension. ..'ibV''" -'-i; ‘ ' i ■' I'iiiiftfttolft 2-6 In this equation, and in the remainder of this section, cell boundary stresses are determined from a simple average of the time n cell centered stresses, e.g.. (t) . (r) . Effect on Radial Motion In this case the velocity u^ is the radial velocity, and the radial components Sj^jn^ of the deviatoric stresses on the various faces of the mass element as shown in Figure 2.2 are s^ •' top of mass element zr ‘ bottom zr right left rr - Sgg(A0/2) front - Sfl„(A0/2) back . The mass of the element is mA0/2Tr, where m is the cell mass. 2-7 17 ^s A, (mu) 4V) (4 - 4 ^ At + T2 ^6 Az At - r,A0 Az At - Sq^AzAt A6 At rr 1 06 or (multiplying by 27 t/A 0) the equation for radial motion is Ag(mu) = 2tt At H - ^i) ♦ s(''> s;_-' r- Az rr L rj^ Az - SqqAz ArJ 2 . 2 . 1 , 3 Energy Equation, (2.12) In the SPHASE term from Eq. (2.12), Ag(mE.j,) = Aty s..u^n.dS s Si^UjUi is the work rate per unit surface area which, for the various faces of the torus, is (s(^)v^^) + cell top \ Z Z Z IT f cell bottom \ zz zr / /s(r)y^^^ + outer cell surface y rz rr / - inner cell surface \ rz rr / where (t) , (b) , (r) , (i) denotes stresses and velocities at the top, bottom, right and left cell faces, respectively. The cell face velocities are determined from a simple average of the time n velocities in the two cells, e.g., (t) , v"(k) ^ v"(k a) 2 , u"[k) ^ u"(klt) 2 The equation for the SPHASE change in total cell energy is Ag(mE, 3(t)^(t) zz zr (b) . 3(b)uCb: zr m - ^i) Ws(r)^(r) , 3(r)„(r)\ 3 , 4 ^ 4 ^ \ rz rr / z ■ ('rz^'^^^^ * ^rr^^^^^) ’^l The post-SPHASE cell specific kinetic energy, 1/2 u^u^, is determined by the post-SPHASE velocities obtained from the equation of motion. The post-SPHASE specific internal energy for the cell, Ej, is then calculated from the updated value of specific total energy and specific kinetic energy. Ej = E.J. - 1/2 u^u^ The partitioning of the internal energy increment among the various materials in a multimaterial cell is discussed in Section 4.6.3. This completes the discussion of the momentum and energy increments due to the deviator stresses as calculated in SPHASE. 2-10 2.2.2 HPHASE - The Effects of Pressure This section discusses the modification of the cell velocity and internal energy due to hydrostatic pressures. The discussion will be in terms of pressures acting on cell boundaries. Sections 2. 2. 2. 4 and 2.2. 2.5 will discuss the definitions of these cell boundary pressures. 2. 2. 2.1 Continuity Equation, (2. 10] No contribution in HPHASE. 2. 2. 2. 2 Equation of Motion, (2.11) The HPHASE contribution to the equation of motion (2.11) A„(mu.) = -At / ""j' The axial and radial n.omentum increments will be discussed separately. Effect on Axial Motion In this case u^ is the axial velocity and nj = -1, 0, +1 are the axial components of the unit normal to the bottom, sides and top of the cell, respectively. The equation of axial motion is therefore A|_j(mv) “ P^^ff|r2 - r^jAt - p'’'*V|: r‘ - rjAt where r., and r, are the radii of the outer and inner cell surfaces, respectively, and P^ , P^ ^are the pressures on the bottom and top cell faces. 2-11 Effect on Radial Motion In this case the velocity Uj is the radial velocity. Referring to Figure 2.2, the radial components of the unit normal take the values n^ = -1, 1, -A6/2, -A0/2, 0, 0, corresponding to the inner and outer surfaces, the two side surfaces and the top and bottom, respectively. The pressure- integral contribution in Equation (2.11) becomes ^ Aj^(mu) = P^^r^AS Az At - P^*'1'2A6 Az At I^^]\r Az A0 At where P^^^^and P^^^re the left, right and side face (s) pressures, respectively. The side pressure P is taken to be the average of the left and right boundary pressures f!^>. Since Ar * r 2 - r^ the equation for radial motion becomes Ap^(mu) = ' 2 ZiTAzAt . 2. 2. 2. 3 Energy Equation, (2.12) In the HPHASE term from Equation (2.12), Aj^(mE^) = -At J Pu^n^dS , S Pu^h^ is the work rate per unit surface area which, for the various faces of the torus, is 2-12 p(a)^(a) .pCb)^(b) pCr)u(r) . p ( 0^(0 cell top cell bottom outer cell surface inner cell surface where the superscripts a, b, r, t. denote pressures and veloc- ities at the top, bottom, right and left cell faces, respectively. Multiplying these terms by the time step and by the surface area they are acting upon gives the HPIIASE equation for the change in cell total energy as r I A^(mE^) = - AtKP^^^v^^^ - 7T (r^ - r^) + (p(r)u(r)r^ - P^*"^u^^^rj) 2 ttAz] . The post-HPHASE cell specific kinetic energy, 1/2 u^u^, is determined by the post-HPHASE velocities obtained from the equation of motion. The post-HPHASE specific internal energy, Ej, is then calculated from the updated value of specific total energy and specific kinetic energy: F = F u.Ui The partitioning of the updated internal energy among materials in multimaterial cells is discussed in Section 4,6.3. 2 . 2 . 2 . 4 Treatment of Free Surfaces and Large bensity biscbntinuities Thus far, the HPHASE equations for momentum and energy increments have been derived in terms of cell boundary pres- sures and velocities. Normally, the cell boundary values are defined to be the simple average of the post-SPHASE cell centered quantities on either side of the boundary, e.g., P(k) P(kb) P(k) + P(ka) ,(r)^ u(k) ^ u(kr) ,U)^ u(k) •»• u(k&) However, these definitions can lead to excessive acceleration of partially-filled, free-surface cells or of cells containing relatively low density material. In the case of a full cell adjacent to a partially- filled, free- surface cell, it would seem that the condition at the boundary between them should relate more closely to the void than the normal averaging procedure would allow. The over-acceleration of the partially filled cells would be prevented by reducing the pressure gradient acting on them. One method of requiring a boundary pressure to reflect the presence of voids and low-density materials more closely is to define the cell boundary pressure to be an inverse compression weighted average of the adjacent cells, i.e., for cells k and kb the boundary pressure between them is defined to be C(kb)P(k) + CCk)P(kb) ^ c(ki)T'+ rcki — where C(k) and C(kb) are defined (for pure cells and free- surface cells having only one material) as follows: C(k) = [m(k)/VOL(k)]/p^ C(kb) = [m(kb)/VOL(kb)]/p„ 2-14 whore m(k) , m(kb) , VOL(k), VOL(kb), are the total masses and total volumes of their respective cells, and p is the reference density. The above compression definitions are meant to be used as weighting factors, and, in the case of one material free-surface cells, assume the material is smeared over the entire volume of the cell. For multimaterial cells the com- pression term is a simple average of the compressions of the materials in the cell. The effect of inverse compression weighting is to weigh more heavily the pressure associated with the underdense cell when defining the cell boundary pressure. This lowers the pressure gradient acting on the underdense cell. It should be noted that if the neighboring cells have the same compression, inverse compression weighting reduces to the simple averaging procedure used before. For boundaries where inverse compression weighting is used to determine the cell boundary pressures, the cell boundary velocities used in the energy equation are defined to be a compression weighted average, e.g., „(b)_ cmufk) + £(kb)u(kb) „(b)_ C(k)v(k) + C(kb)vCkb) " C(k) +”C(kb) ' • This definition weighs more heavily the velocity of the denser cell and assumes its velocity is probably more realistic than the velocity of the less dense cell. The process of inverse compression weighting for boundary pressures and compression weighting for boundary velocities is used at all boundaries of a free-surface cell. However, since the over-acceleration problem can occur when- ever a dense, high-pressure cell is adjacent to an underdense. 2-15 low-pressure cell, an input variable, CRATIO, has been added to the code in order to allow the user to control the use of the compression weighting. Whenever the compression ratio between two cells is greater than CRATIO, the pressure at the boundary between them is taken to be the inverse compression weighted average of the cell-centered values, and the boundary velocity is taken to be the compression weighted average. 4 CRATIO has a default value of 10 which essentially restricts the use of the weighted averages to the free surface. 2.2.2. 5 Artificial Viscosity The stability of a HELP calculation depends, in part, on the effective viscosity which results from converting kinetic energy into internal energy during mass transport. If velocities are small, there is little mass flux between cells, and as a result, little effective viscosity. Since these conditions have arisen in a variety of calculations, an option, LVISC, has been added which, when set equal to one, will signal the code to add an artificial viscosity term to the cell boundary pressures in subroutine HPHASE. The discussion which follows presents a form of artificial viscosity which has been useful in several calculations. How- ever, when the code is applied to a new class of problems, it may be that other artificial viscosity formulations will give improved results. The user can make such changes if the need arises . fa) The definition of the artificial viscosity, Q' , which can be added to the pressure at the top cell boundary is ( v(k) - v(ka)). ^ [ P (k) + P (ka) . P(k) V. P(ka) i 2' 2 2-16 where p is material density. An analogous term, added to the pressure at the bottom cell boundary, can be defined by substituting kb for ka in the above equation. For the right cell boundary, the artificial viscosity * 1 11 (!<] ■ 11 (kr)j . P fkl + P (kr; p(k1 + P(kr: (j,) and an analogous term, Q , is defined for the left cell boundary by substituting kil for kr. If one of the cell boundaries is also a transmittive grid boundary, the artificial viscosity term for that cell boundary is zero. If, instead, one of the cell boundaries is a reflec- tive grid boundary, the artificial viscosity for that cell boundary is -2 . o)(k) Vl P(JOl • P (kT , where co is the velocity component normal to that boundary. 2.2.3 TPHASE - The Effects of Transport The purpose of this section is to describe how the trans- port of mass, momentum and energy from cell to cell is accounted for in the code. This is done by calculating the integrals in the last terms of Eqs. (2.10) through (2.12). In the dis- cussion below, attention is given to the case of transport between pure cells. The special provisions required to treat interface cells are described in Chapter IV. 2-17 2*2. 3.1 Continuity Equation. C2.10] From Eq. (2.10) the mass being transported is A^(m) = - At y* pu^n^dS and is determined for each cell face from 6m = - u^A^At where p° is the density of the cell from which the mass moves (donor cell), is the area of the cell face and u^ is an inter- polated value of the velocity component normal to the cell face and represents the velocity at the interface at the end of the time step. For example, considering transport from cell k into the cell above it, ka, the cell boundary velocity, v^ \ is defined to be a simple average of the velocities of the two cells. However, we want the mass transport velocity, v, to be the cell boundary velocity at the end of the time step. At. This transport velocity therefore is the velocity of the material located a depth Az = vAt into the donor cell k (see Figure 2.3). However, since the transport velocity v, when adjusted by the axial velocity gradient, must equal the cell boundary velocity, we have V (1 At) = V 2-18 ka f a1 Figure 2.3 --The relationship between v'‘ , the cell boundary velocity at the beginning of TPHASE, and v, the transport velocity. [i > v.LKaj - vjjO At1 Az J . Similarly, the radial mass transport velocity between the donor cell k and its neighbor on the right, kr, is defined to be j [u(k) <• u(kr)] “ ■ |~i ^ ■'iua At] Calculated transport masses are subtracted from the donor cell mass and added to the acceptor cell mass, and the cell momenta and internal energy are updated accordingly. Since all transport terms are computed using the post-HPHASE quantities, a given cell is updated only after the transport terms have been calculated for all four of its boundaries. 2 . 1 . 1.1 Equation of Motion, (2.11) The TPHASE contribution from the equation of motion (2.11) is Aj(mUj) The axial and radial momentum separately. (pu.Uj)n.dS . increments will be discussed I Effect on Axial Motion At each of the four cell faces comprising the surface of the cell, the specific axial momentun, u^ , being transported Is the axial velocity of the cell from which the mass moves, i.e.. 2-20 77 I 1 1, i I I Uj “ v(kcl) where kd is the index of the donor cell. Since this velocity is constant over the cell face it can be removed from the integral sign and we obtain A.j.(mv) = v(kd) The quantity in brackets is seen to be the mass transported through the cell face as calculated in Section 2. 2. 3.1. Thus for each cell face. ^-At J pu^n^dsj A.j.(mv) = v[kd) 6m The total axial momentum increment for the cell is the sum of the increments obtained for each of the four cell faces. Effect on Radial Motion The specific radial momentum, Uj , being transported is the radial velocity of the cell from which the mass moves, i.e. , Uj = u(kd) where kd is the index of the donor cell. By an analysis similar to that for the axial momentum increment, the radial momentum increment is A.p(mu) = u(kd) 6m for each of the cell faces. The total radial momentum incre- ment for the cell is the sum of the increments obtained for each of the four cell faces. i«rT Tii-v; Vk. • '• 2-21 2 . 2 . 3 . 3 Ene luation The expression for the transport of energy from Eq. (2.12) i^) = J (pu^E^) n.dS . A,p(mE To evaluate this integral, the transported specific energy^ E^^ across each cell face is taken to be that of the donor cell, kd, i . e . , - EjCkd) * \ [(u(kd))^ * (v(kd))^] and the total energy which is transported across a given inter- face is therefore the product of this specific energy and the associated transport mass, which was computed above in the continuity equation, i.e. , d Aj(mEj) » E^6m The total energy increment for the cell is the sum of the increments obtained for each of the four cell faces. The internal energy is obtained by subtracting the new "kinetic energy (obtained from the updated mass and momenta) from the updated total energy. The partitioning of the updated internal energy among materials in multimaterial cells is discussed in Section 4.6.3. 2.3 MATERIAL MODEL As is well known, tho conservation equations (Eqs. 2.1 through 2.3), which consist (in two dimensions) of four equations, are insufficient to allow determination of the eight unknowns (p, u, v, P, E™, s , s , s ) contained in them and must be supplemented by equations defining the 2-22 irirf'ff'! ‘'WMrfrfy.' material model to be used. The material model employed in the HELP code can represent inert materials or high explosives (HE) and provides a fairly sophisticated HE detonation routine allowing for the detonation of multiple HE packages and for detonation around inert obstacles. The code also includes strength effects, modeling them in the elastic-plastic regime, and a simple tensile failure model. By bypassing the strength effects subroutine, SPHASE, the code solves problems in com- pressible fluid mechanics. A more complete discussion of these points follows. 2.3.1 Equation of State The HELP code employs three different equations of state: the Jones-Wilkins-Lee (JWL) equation of state high explosives, the Tillotson equation of state for for for inert materials, and the ideal gas equation of state. The array dimensions of the HELP code allow up to 30 materials in any problem. The following convention is used to relate the 30 material numbers to one of the three equations of state: materials 1-19 are assumed to use the Tillotson equation of state; material 20 is assumed to be an ideal gas; and materials 21-30 are assumed to be high explosives and use the JWL equation of state. 2. 3.1.1 High Explosives Once an explosive cell is detonated, its pressure is computed using the Jones-Wilkins-Lee (JWL) equation of state. ^ The equation is P(Ej,p) = A (l - + B 1 - + apEj L 2-23 where specific internal energy in ergs/g density of detonation products in g/cm^ density of undetonated explosive in g/cm^. The constants p^, A, B, a, a and 3 for four explosives are listed in Table 2.1. The *onstants fur these four explosives are defined via DATA statements in EQST, and the material numbers of these explosives are: 21-Comp B; 22-TNT; 23-OCTOL; 24-PBX 9404. The JWL equation of state has a pressure maximum which is potentially hazardous to the pressure iteration (see Section 4.6.1). The pressure iteration adjusts the densities of the materials in a multimaterial cell to achieve pressure equilibrium. During this iteration process, it is possible for the code to reference the JWL equation of state with a density corresponding to a pressure beyond this maximum, thereby preventing convergence of the iteration. In order to preclude this possibility, a variable, RHOMAX, is defined via DATA statements in subroutine CDT, This is the maximum value of density for each HE to be used in the pressure iteration. The values of RHOMAX for the four specified explosives are listed in Table 2.1. The sound speed, C, in the detonated explosive is de- fined by the equation f m - Pp' c = Vr P/p • The value of r (defined by DATA statements in subroutine CDT) as well as the detonation velocity, D, (defined by DATA state- ments in subroutine DETIME) and detonation energy, E^ , (de- fined by DATA statements in ADDENG) for the four named explo- sives are listed in Table 2.1. 2-24 m'"Trr I I'li'Triil i " iiniiiai|iiiniiiiiiiiiiiMiiiiin.i 2. 3. 1.2 Inert Materials The equation of state used in HELP, for initially con- densed materials, is that due to J. H. Tillotson, modified to give a smooth transition between condensed and expanded states. For the condensed states, i.e., when p/p^ > 1, or for any cold states, E^<Eg, the equation has the form P = = a + V ^ » 1 Ejp + Ay + By' J where , ! Ej * specific internal energy in ergs/g p * material density in g/cm^ n * p/pQ y " p/Pq - 1 For expanded hot states, i.e,, when p/p^ < 1 and Ej > E' the equation of state has the form ff K'. » aEjP + l)E,P ^ ^ + Aye -B(Pq/p-1) LEoH -a(p^/p-l) A smooth transition between the condensed and expanded states Is insured by a transition equation for the intermediate region defined by E < E < E" and p/p < 1. This blended portion of sis o the equation of state has the form „ CK,-E,) Pj ♦ (E'-Ej) P, 2-26 i In these equations, a, b, a, 3, E^, E^, E^, A, B, and are constants for the particular material. The values of these constants for nineteen materials are given in Table 2.2. For multimaterial cells, an iteration is used to deter- mine the cell pressure and the densities of the various materials within the cell. This iteration is discussed in Section 4.6.1. 2. 3. 1.3 Ideal Gas As noted earlier, material number 20 is assumed to be an ideal gas and the pressure is calculated from the ideal gas equation of state P = (Y-1) P Ej and the sound speed given by C » with p being the density, Ej the specific internal energy, and Y the ideal gas constant. The normal density, p^, of the ideal gas is defined in a DATA statement to be that of air at standard temperature and pressure (0.00129 g/cm^) and is used only for editing purposes. Since y is stored in a non-dimensioned variable, GAMMA, the code has the capability of handling only one ideal gas. Note, however, that the JWL equation of state reduces to the ideal gas equation of state if A * B ■ 0 and a = Y-1* Also, in order for a material using the JWL equation of state to have the ideal gas sound speed, it is necessary to have T = y* 2-27 n* </) 00 « 09 kf^oom«foor^o o ooooooo UJ M • • O i/)'0O4/)iOfNOt9tSOOt«O 00 00 O <C 9 9 i/> O ■ ^C>«of^ooi/>a»r'ih'»oa» r^or'-Kioor-f>j Ci, u ^ Aoor«r4«-i^ooo»H^O fN rj im « in r* rj (A 06 -4 00 « sO jj ^ ^K)«toOi/)ooooor4n» m t/> <n <n m i/) o • u <u| Oi/)ini/)i/)i/)L/)i/)Oi/)i/) i/) 4/> lA u> i/) m u) m a I 0l/)br)4/>l/li/>t/)i/)0)»00 lA (A lA lA l/> i/> l/> 4A O 0« •H *>s H </t 4A iA lA U1 (A iA M riMOiAr^r^aj^r^rsifN 1h m»aoo»^oooooo ^ oddododdddo o ooooooo >v o «/) «H if X lAiAiA lA iAf4 p lA«-)OOiAiA«A'0^<-40 fNrH^OO0-«^0O0 OOOOOOO tf) M O ^ O OOOsJOlAr'K)— i-^K)\0tv ^ otor^r^'-'Oor^iA^o T3 ►o^^O*-Hi-*-HrNOOO ooooooo •H <-• i-H ^ O O ^ • 4AlAinLAlAtAiAiA^n<&> iA iAiAiAiAiAIAtA OOOOOOOOOOO O OOOOOOO iA iA 04 44 >-4 4; c if 44 44 c« V o *J 3 p £ -*4 *4 «> •r4 tf)HHtABV>«4 C 41 O 04 *4 jrt rt •04^Xf-i^e^ rj rt "U *-> X >-•’-• E —« o 34>^B-HHO,C^X ^4 CVh-4p.H(0 in uu,<«hs!:3EPq.o o < m OOP -j x «-* irNfO^LAvOr OOOO^ *N fO^tiAOl OOCTl 2.3.2 Clastic-Plastic Constitutive Relation Section 2.2.1 discussed the modification of the material velocity and internal energy fields due to the effects of the stress deviators in subroutine SPHASE. This section will dis- cuss the elastic-plastic material model incorporated in the HELP code to calculate the strass deviators. Subroutine SPHASE has two major tasks: first, it computes instantaneous strain rate deviators and updates the stress deviators, and second, it updates the cell momenta and energy to account for the effects of the stress deviators (as discussed in Section 2.2.1). This section will discuss the first task of SPHASE: the computation of the strain rate deviators and the updating of the stress deviators, including a description of the von Mises yield condition. 2. 3.2.1 Strain Rate Deviators As a first step in computing stress deviators, the instantaneous strain rate deviators are computed from the velocity field, i.e., for cylindrical coordinates • e * V 1 (u + V + -) zz y O X y Sr = S 1 ■ T (u + V X y ^ -) X-' Sz ■ (“y * . where the x and y subscripts denote partial differentiation with respect to those variables. Here the necessary space derivatives of the velocity field are computed (centered at cell centers) in a straightforward manner, i.e., 2-29 The above calculation gives only the strain rate deviator for that material which is at cell center at the beginning of the time step, and it is necessary to recognize that constitutive equations are to be applied only along a particle path. To this end, one must take a convective derivative. Consider the time n position q” of a particle which, at time n+1, will be at a cell center. This point will be offset from the cell center by 6x and 6y where 6x ■ -u(k) At 6y ■ -v(k) At Ay 2-30 Now we wish to determine, by interpolation using the value of the strain rate deviators at cell centers, the strain rate deviators at the offset point q’^. This is done by weighting the contribution of neighboring cells in proportion to their overlap areas with a rectangle of cell dimensions, centered at Q^. For example, referring to the above sketch, the weighting factors are for cell k I WK *= (Ax - |6x|) (Ay - l6y|) for cell in horizontal direction WKH = |6x| (Ay - |6y|) for cell in vertical direction WKV = |6y I (Ax - |6x j) for cell in diagonal direction WKD = |6x| |6y| so that the resulting interpolated value of a strain rate deviator component, say is given by (WK) ^^^(k) + (WKH)e^^(kh) + (WKV)E^^(kv) + (WKD)ey^(kd) Having thus determined , e' , i' for the mass which will zz rr rz be at the center of cell k at the end of the time step, we can now calculate the strain deviator increments. For example , Ae^_ = i' At rz rz 2-31 2. 3. 2.2 Stress Deviators These strain deviator increments, ^re next used in the material constitutive equation to update the stress deviator tensor. To do this one needs the time n stress deviators for the =ij particle in question and these s'^^j are determined (from stored cell centered values of the s^^j at time n) by an identical area weighting procedure as described above for strain rate deviators, For an elastic-plastic material the following calculation then (provisionally) updates the s^j according to n+1 ij ■ =ij * 2G lej. where G is the material rigidity modulus; i.e., the three components of interest are given by s^'*'^ = s' + 2G Ae zz z_ zz = s' + 2G Ac rr rr rr = s'^^ + 2G Ae^^ rz rz rz The provision (above) is that the yield condition for the material not be violated. 2-32 2. 3. 2. 3 Yield Criterion HELP employs the von Mises yield condition which requires that s. .s. . < 2Y‘- ij 11 - where s s- • is twice the second invariant of the stress tensor^ ij ij 1 • C • f ®lj =ij " =rr ^ * *63 " where s-, oO (s + s ) is the hoop stress and Y is the '• rr zz shear yield strength. If the yield condition is violated during a time step, the deviatoric stress components are proportionately reduced normal to the yield surface. If the yield condition is not violated, the deformation is purely elastic and the stresses are not reduced. A variable yield strength Y = t^o ^1 ^ * ^2 ^ 5 - Ej/y is defined in the code (in subroutine STRNG) to account for the increase in strength at high compressions and the decrease in strength at elevated energies. In this equation Y^, Yj^ and Y 2 are input constants, p = p/p^ - 1 where p is the density and the reference density, is the specific internal energy and is the melt energy. This expression is never allowed to be negative, however, so that once the melting point energy is exceeded the material yield strength is set to zero. It siiould oe noted tiiat cells (including multi- material cells) containing material which has failed in tension (see the following subsection) are assumed to have no shear yield strength. Further modifications to the yield criterion for multimaterial cells are discussed in Section 4.6.2, 2-33 2.3.3 Material Failure in Tension The tensile failure criterion for pure cells is a simple one based on relative volume. When the relative volume (Vq/V = p/Pq) of material in a pure cell reaches a value which is greater than a specified maximum dilitation (an input con'stant for each material, AMDM-), the material in the cell is considered to have failed. A cell which has failed is not allowed to sustain tensile hydrostatic pressures nor, as noted in the preceding subsection, to have shear yield strength. The tensile failure criterion is more stringent for multi- material cells in that they are not allowed to sustain negative pressures, regardless of the relative volumes of the materials they contain, since, presumably, the negative pressure would cause the materials to separate at the interface, thus immediately relieving the tension. 2.3.4 Detonation of High Explosives To model high explosives a detonation procedure has been incorporated into HELP as well as special prescriptions which maintain a sharp boundary between the detonated and undetonated explosive . 2. 3. 4.1 Detonation Procedure The detonation procedure incorporated into HELP makes no assumptions about the geometric shape of the explosive package, the number of explosives in the problem, or the number of initiation points, although certain variable dimensions allow a maximum of ten primary and ten secondary initiation points. 2-34 ?v t Primary initiation points may be located in any con- figuration within an explosive package, and may have a time delay associated with them. A primary initiation point can be detonated at any time on or after time zero (the start of the calculation). For example, a string of five detonators may be initiated sequentially or simultaneously, depending upon the desired results. Secondary initiation points which are used to detonate around obstacles may also be located in any configuration throughout the high explosive packages. The location of secondary initiation points should be selected carefully if they are to perform the intended function, that of detonating around obstacles. The initiation times for each of the secondary initiation points are computed by the same method used to compute the detonation time for a cell. Briefly, the numerical method for handling the detonation of a high explosive is to calculate the detonation time of each explosive cell in the grid based upon the cell's location and the constant detonation velocities of each high explosive through which the detonation front passes between the primary sources and the cell to be detonated. The high explosive detonation routine employed in the HELP code allows multi - initiation of multiple, adjacent, high explosive packages. For the case of adjacent HE packages the current method is restricted to explosives whose detonation velocities and C-J pressures are within ten percent of one another so that the degree to which one of the detonations is overdriven is minimal and can be neglected. Figure 2.4 illustrates the method employed to compute the detonation times of cells in adjacent high explosive packages. In this example the primary initiation point, P^, is in the HE. package, and cell K1 is a pure cell in the HE, A, <4^ package, cell K2 is a pure cell in the IIE 2 package, and cell K3 is an interface cell containing both HEj^ and HE 2 . In all cases the time to detonate a cell is the sum of the time in- crements, t^, associated with the time for the detonation front to travel from the points to which lie on the line from the primary initiation point to the center of the cell. If this line is at an angle which is less than or equal to 45° with respect to tht' x-rixis, the P. arc the inter- cepts of that line with the centers of the grid columns (e.g., the line to cell K3) . If the angle is greater tlian 45°, the are the intercepts of that line with the centers of the grid rows (e.g., the lines to cells K1 and K2) . If both points, P^ j^ and P^, are in pure HEj^ cells, the detonation time, t^, is based on the IlE^^ detonation velocity, Dj. If the second point P^, is in an interface cell containing both HEj^ and HE 2 explosives, the detonation time, t^, is based on the maximum of the HE^ and ME 2 detonation velocities: ^°1» °2^max* Referring to Figure 2.4, all of the time increments for cell K1 are computed using . For cell K2, however, the first three increments are based on , the next two on (0i>^2^max’ and the last one on D 2 . For cell K3, all but the last are based on and the last is a function of C^i>^2^max’ The HELP code has been applied to the jet formation process associated with lined shaped charges in which the detonation wave has been "shaped" by the insertion of inert material into the high explosive package In order to accomplish this it was necessary to incorporate the capability of detonating around obstacles into the code. Figure 2.5 shows a disc of inert material situated in a package of HE in which detonation is initiated at a point on the axis of symmetry. The primary initiation point is shown as Point 1. The detonation time of each cell contained in 2-37 area I is based on its distance from Point 1 and the detona- tion velocity of the HE. Similarly, the cells in areas II and III will be detonated at times based on their distances from the secondary initiation points, 2 and 3, and the detonation velocity. Points 2 and 3 in Figure 2.5 act as secondary initiation points for their respective areas as soon as the HE in the cells containing these points has been detonated . Both the primary and secondary detonation points, as well as their corresponding detonation regions, are part of th.e initial conditions and are specified by the input cards. At the beginning of the problem the time to detonate an explosive cell is computed and stored in the DETIM array (see Section 7.2.7). On the cycle for which the value of T in the calculation equals or exceeds the detonation time associated with a given cell, the chemical energy of detonation of the explosive (E^ in Table 2.1) is added to the cell. If the cell contains more than one explosive, the appropriate chemical energy of detonation of each high explosive is added to the cell. Each cycle, ADDENG prints the row (I) and column (J) of every cell detonated on that cycle. The total energy released on that cycle is also printed, and the theoretical energy of the grid (ETK) is updated accordingly. 2. 3.4.2 Definition of the Detonation Front Because the detonation procedure used by HELP detonates all of the explosive in a cell at one time, the detonation front is approximated by the boundaries between detonated and undetonated cells. In order to prevent small numeri- cal signals from propagating through the grid ahead ' 1 %/, ■■ , I H 2-39 of the detonation front, it is necessary to place certain restrictions on cells containing undetonated explosive dur- ing the three phases of the code. There are two restrictions necessary to prevent the pressure effects from propagating undesired numerical sig- nals. The first of these is implemented in CDT and requires that any cell containing undetonated explosive have zero pressure. For these cells, the equation of state routine, EQST, is not called by CDT. This restriction is necessary since the explosive equation of state applies only to the detonated explosive. The second restriction, implemented in the pressure effects subroutine, HPHASE, requires the pressure and velocities to be zero at the four boundaries of an interface or pure cell containing undetonated explosive, even if one or more of the neighbor cells contains detonated explosive. In the strength effects phase, SPHASE, any mixed or pure cell which contains explosive is assumed to have no strength and its deviator stresses are set to zero. This assumption applies to cells containing detonated as well as undetonated explosive. In subroutine TPHASE, the transport phase of the code, no mass, momenta, or energy is transported into or out of a cell containing undetonated explosive. These transport restrictions are coded into DMCALC as well, which computes the mass of each material to be transported across the four boundaries of an interface cell. If an interface cell contains undetonated explosive, no mass of any material is transported into or out of the cell. 2-40 ■t?favy-icj>-»i!fe>-^.<itMu^ffaw-«^»->vi--'--iJ»'^ ' ***■■— — i^smus^ CHAPTER III GRID BOUNDARY CONDITIONS 3.1 BASIC ASSUMPTIONS Additional specifications are required for cells which border the grid boundaries because necessary quantities are not defined for neighbor cells which would be outside the grid. The code provides for a transmittive or reflective boundary condition at the bottom and a reflective boundary condition at the left (see Figure 3.1). Boundary conditions for grid bound- ary cells are then derived by assuming fictitious neighbor cells outside the grid. For transmittive boundaries the flow variables are the same in the fictional cell as in the border cell, and for reflective boundaries the state is assumed to be the same in the fictional cell except that the velocity component normal to the boundary has the opposite sign. These assumptions, as they apply to subroutines ‘ oE, HPHASE, and TPHASE, are described in more detail in the sections that follow. Special prescriptions for moving tracer particles on and near the grid boundaries as well as special edited quantities which indicate the work done and the mass, momentum, and energy fluxes at the grid boundaries are described as well. 3.2 STRENGTH PHASE (SPHASE) The SPHASE grid boundary conditions will be discussed as they pertain to the three subtasks of SPHASE: 1. the calculation of cell centered strain rates; 2. the interpolation of cell centered stresses and strain rates to account for convection; 3. the application of deviator stresses to update cell velocities and internal energies. 3-1 Within each subtask the transmittive and reflective boundary cases will be discussed. 3.2.1 Definition of Strain Rate Derivatives for Cells at a Grid boundary The strain rates are computed using space derivatives centered at cell centers. For example the finite difference analogs of 8v/9y and 3u/3y for cell k are 3v _ v(ka) - v(kb) 3y “ -SAy. ^ + Ay . + .SAy.^^ 3u ^ u(ka) - u(kb1 3y “ .SAy..^ ♦ Ay^ + -SAy^^^ ’ with ka and kb being the indices of the cells above and below cell k, respectively, and Ay being the axial dimension of the corresponding cell. However, if cell k is at the top trans- mittive grid boundary, these space derivatives are approxi- mated by 3y -5Ay^._^ + v(kl - v(kb) Ay/, + l.SAy, u(k) - u(kb] •5Ayj_;^ + •SAy^ which assumes that the fictitious cell above cell k has the same velocity components and axial dimension as cell k. If, on the other hand, cell k is at the bottom grid boundary and that boundary is reflective, these space derivatives become 3-3 3v _ y(ka) r (-v(k)) 8y ■ l. SAyV + . SAy/^j 3u _ u(ka) - u(k) 3y " l.SAy^ + .SAyV^^j which assumes that a fictitious cell below the grid has the same velocity components and axial dimension as cell k, except that the velocity component normal to that boundary has the opposite sign Similar approximations are made for 8u/3x and 3v/3x for cells at the right transmittive grid boundary and for cells at the left reflective grid boundary. 3.2.2 Definition of Interpolated Strain Rates an d Stresses for Cells at a Grid Boundary SPHASE computes the deviator stresses of a particle that will be at the center of cell k at the end of the time step. The strain rates and time n stresses of such a particle are computed by taking an area weighted average of cell centered strain rates and stresses. (See Section 2. 3. 2.1.) If cell k is at a grid boundary, and the velocity field is such that the material is flowing from outside the grid, then it is assumed that the fictitious cells outside the grid have the same stresses and strain rates as the contiguous cells in- side the grid. 3.2.3 Definition of Velocities and Deviator Stresses at Grid Boundaries The cell velocities and internal energies are updated in SPHASE using cell boundary stresses and velocities which 3-4 i. . are averages of the cell centered quantities. For example, the stresses and velocities at the right boundary of cell k» which is not at a grid boundary, are -r . u(r) * u(kr: -r . vfk) * v(tr: _ s^^(k) * s^^(kr) S “ ■ ' ' ' ' • rz 2 If, however, cell k is at the right transmittive grid boundary then the stresses and velocities of the fictitious cell outside the grid are assumed to be the same as those of cell k, and the boundary stresses and velocities are r.r _ ,, flr-J u" = u(k) v^ * v(k) ^rr ^rr^^^ And if cell k is at the left reflective grid boundary, the following equations define the stresses and velocities at the left boundary: vCk) = 'rrC') 3-5 i.e., the shear stress and the velocity component normal to the boundary are set to zero, and the normal stress and tangential velocity component are assumed to be equal to the cell k values. Likewise, the stresses and velocities at the bottom grid boundary, when it is reflective, are u^ = u(k) .-.b = 0 = 0 rz s„00 3.2.4 Correction to the Theoretical Ener one at for Work The theoretical energy of the grid, ETH, is updated only for the work done at the transmittive grid boundaries. At the right boundary <5E^ = [s^^(k).u(k) + s^^(k).vCk) )-Aj.-At where s^^, s^^, u, v, A^ are the deviator stress normal to the right boundary, the shear stress, the velocity component normal to the right boundary, the velocity component tangent to the right boundary, and the area of the right cell face, respectively Similarly, at the top transmittive grid boundary ^t “ [s^^(k)*v(k) + s^^(k)-u(k)] *A^-At and at the bottom grid boundary, when it is transmittive 'vH-ff / iwjjr !■*» *•*»' '"•‘.t " ■ “ ts^^(k).v(k) + s^j,(k) -u(k)] .A^j.At When the bottom boundary is reflective, = 0, since in that case the normal velocity, v, and the shear stress, s^„, are zero, rz 3.3 HYDRODYNAMIC PHASE (HPHASEl The effects of pressure are computed by using pressures and velocities defined at cell boundaries. If a cell boundary is part of the grid boundary, the definition of the pressure and velocity at that boundary depends on whether or not that boundary is reflective or transmittive . 3.3.1 Definition of Velocities and Pressures at transmittive Grid Boundaries If the grid boundary is transmittive, the cell boundary pressure is assumed to be equal to the cell centered pressure. Likewise, the velocity component normal to the boundary is assumed to be equal to the cell centered velocity component. For example, if cell k is at the top transmittive boundary P(k) V® = v(k) and if cell k is at the right transmittive boundary = P(k) u u(k) Because the transmittive boundary condition assumes all material outside the grid is at the same state as the material in the border cells of the grid, the user should, by rezoning, prevent a strong shock from moving across a trans- mittive boundary. In such an instance a numerical signal can be generated by the boundary condition, and the stresses in the material behind the shock can become noisy and meaningless. 3.3.2 Definition of Velocities and Pressures at deflective tjrid boundaries If the grid boundary is reflective, the cell boundary pressure is assumed to be equal to the cell centered pressure. However, the velocity component normal to the boundary is assumed to be zero. For example, if cell k is at the left reflective grid boundary (an axis of symmetry in cylindrical coordinates) = P(k) and if cell k is at the bottom grid boundary which is acting as a reflective boundary (CVIS = 0), P^ = P(k) 3.3.3 Correction to Theoretical Energy for Work Done at ffrid boundaries in HpHAsIE The theoretical energy of the grid, ETH, is updated in HPHASE for the work done at the transmittive grid boundaries. At the top transmittive boundary 3-8 6E ^ = -P (k>v(k>A^-At , and at the right transmittive boundary *= -P (k>u(k>A* At , where A^ and A^ are the areas of the top and right cell faces, respectively. If the bottom boundary is transmittive = P(k>v(k>AgAt , and if it is reflective, dE^^ * 0, since -b 0 . 3.4 TRANSPORT PHASE (TPHASE) Material is transported through the grid by computing mass, momentum and energy fluxes at each of the four boundaries of every cell. If one of these boundaries is a grid boundary, these transport terms are determined by the following con- ventions for transmittive and reflective boundaries, and the theoretical energy is adjusted accordingly. 3.4.1 Transport of Mass, Momentum and Energy Across Transmittive Grid Boundaries The velocity of the material at and beyond a grid boundary is assumed to be the same as that of the boundary cell. The transport velocity used to compute the mass transported across a grid boundary is therefore the cell centered velocity. How- ever, if the cell centered velocity component normal to the grid boundary is such that the material is moving into the grid, the mass transport term at that boundary is set to zero. For example, if cell k is at the top transmittive grid boundary and v(k) > 0, then 3-9 6nia = v(k) . At. A^. p(k) On the other hand, if v(k) $ 0 then = 0. Likewise, if a cell k is at the right transmittive boundary and u(k) > 0, then 6m^ = u(k) .At. A^. p(k) If u(k) £ 0, however, = 0. For the case where 6m > 0, the momentum and energy trans- ported across grid boundaries are debtrmined by the donor cell method, i.e . , 6 (mu) = 6m>u(k) 6(mE) = 6m.E(k) The code does not have provisions for feeding material into the grid across any of the three transmittive grid boundaries; however, TPHASE can be easily modified by the knowledgeable user to do so. 3.4.2 Change of Momentum jgt Cells at Reflective Grid Boundaries in TPHASE Even though no mass is transported across a reflective boundary, the momentum of a cell at a reflective boundary is modified in TPHASE. A useful way to think about reflective boundaries, such as the axis or grid boundaries, is to take them to be planes (or a line) of symmetry through which equal and opposite transports occur in both directions. 3-10 Assume, for example, that the mass transport calculation at the axis indicates that a mass 6m* would pass through the axis (i.e., u < 0) . Then an equal and opposite transport is assumed from the left. The net effect on cell mass is zero. However, there is an effect on the cell momentum. At the axis, momentum 6mu is lost and a momentum 6m(-u) is gained. The net result is a momentum change of -26mu. Thus the momentum equation has to have this correction term; specifically A^mu ^ 6 (m^^Uj^) - 26mu , i where the subscript i denotes the four sides of the cell. Similarly, it can be shown that the correction term at the bottom grid boundary, when it is reflective, is -26mv. In the case of multimaterial cells, the energy equation must also be modified. For example, consider a cell on the axis of symmetry, where u is its pre-TPHASE radial velocity and u* is its post-TPHASE radial velocity. The correction to the cell's thermalized kinetic energy (see Section 4.6.3) is as follows: - ^ (u-u')^ + (-u-u')^ * 26muu' where the first term corresponds to the mass that "left" with velocity u, and the second term corresponds to the mass that "entered" with velocity -u. Similarly, the correction term for a multimaterial cell at the bottom grid boundary, when it is reflective, is 26mvv'. * The discussion above is independent of how one calculates dm. However, in HELP, at the axis 6m ■ pu(iTArAz) At if u < 0, and 6m ■ 0 if u > TT7 The area term, irArAz, is the area of the surface whicK cuts through the center of the cell. 3-11 3.4.3 Correction to Theoretical Energy for Energy l^ransported Across Grid Boundaries in TPHASE In TPHASE the theoretical energy of the grid, ETH, is adjusted to account for the total energy of the mass that leaves the grid across the transmittive grid boundaries. For example, if cell k is at the top transmittive grid boundary, and 6m^ > 0, then 6Et = 6m^-(E(k) + 1/2 (u(k)^ + v(k)^)) where E(k) is the specific internal energy of cell k. There is no mass loss and therefore no energy loss at a reflective boundary. 3.5 TRACER PARTICLE MOTION NEAR GRID BOUNDARIES As elsewhere in HELP, subroutine MOVTCR assumes the velocity of material outside a transmittive grid boundary is equal to that of the material in the cells bordering the grid boundary. Therefore, if the overlap area (described in Section 4.1.2) used to compute the velocity of a given particle extends into a fictitious cell outside the grid, the velocity of that fictitious cell is assumed to be equal to its neigh- bor inside the grid. For example, as illustrated in the following figure, the radial velocity of the particle is defined as follows: Aj^u(k) + A 2 u(k) + AjU(ka) + A^u(ka) u = — — — Ai + A2 + A3 + A4 3-12 ka ! 1 1 1 r^4 5ii. i lA, 1 ^ • A2| j 1 ^ i 1 • i. .J 1 F ic V Riglit Transmitt Boundary Cells ive Grid I£ the grid boundary is reflective, however, the velocity component normal to that boundary in the fictitious cell is assumed to have the opposite sign. If the particle is near the axis of symmetry, as illustrated in the following figure, the radial velocity of the particle is defined as follows: A^u(k) A2 (-u(k)) * Ag (-u(ka)) + A^u(ka) Aj ♦ ♦ Aj . A^ Fictitious Cel C J r — 1 1 1 1 1 • ka 1 -Ti KZ 1 }^2 1 1 1 1 1 1 ! I . • k ! 1 1 1 i Axis of Symmetry An analogous procedure is followed for computing the axial component of velocity for a particle near the bottom grid boundary when it is reflective (CVIS = 0) . The coordinates of the tracer particles are stored in grid units. If the new position of a particle is computed to be beyond the grid, its x-coordinate and/or y-coordinate is set equal to the grid boundary line it has exceeded, i.e., 0, IMAX, or JMAX, since it is not possible to track a particle outside the grid. Furthermore, once a particle is moved onto a grid boundary, it is not moved off of that boundary, only along it. 3.6 EDITING OF FLUXES AT GRID BOUNDARIES In SPHASE and HPHASE the total work done at each of the transmittive grid boundaries is accumulated, as well as the internal energy lost by each material package due to work done at these boundaries (calculated in subroutine EOUT) . The total work is printed by subroutine EDIT under the heading WORK DONE and is included in the total printed under the heading IE OUT. In TPUASE, the mass that leaves the grid and its asso- ciated total enfjrgy and momenta are accumulated and are printed by subroutine EDIT under the headings MASS OUT, ENERGY OUT, MU OUT, and MV OUT, respectively. The kinetic and internal energy lost by each material package due to trans- port across grid boundaries is also accumulated and is in- cluded in the totals printed by EDIT under the headings KE OUT and IE OUT. . 3-14 CHAPTER IV MATERIAL INTERFACES AND INTERFACE CELLS The following sections describe the HELP method of de- fining and propagating material interfaces through an Eulerian grid and discuss the special treatment given to multi- material cells. The use of tracer particles to define the in- terfaces and the manner in which these particles are moved are discussed in Section 4.1. The method for storing information on several materials in an interface cell and for converting a pure cell into an interface cell is described in Section 4.2. The procedures for calculating fractional cell face areas (used in the transport of materials across interface cell boundaries) are given in Section 4,3. The adjustment of mass transport terms to exactly evacuate a material when an inter- face leaves a cell and to prevent overemptying of a material is discussed in Section 4.4. In Section 4.5 the character- istics of the three types of interface cells - multimaterial, free surface and slipline - are given. Finally, the special treatment necessary to define various cell and material prop- erties for multimaterial cells is given in Section 4,6. 4.1 DEFINITION OF MATERIAL INTERFACES 4.1.1 Use of Tracer Particles A set of massless tracer particles are initially posi- tioned along the boundary of each material package. These tracers are numbered such that the package is on the left as one proceeds between any two consecutive tracers. The material interface therefore is defined by the line segments connecting these tracers. The initial spacing of the tracers is specified when they are generated by TSETUP and can be maintained in a specified region as the calculation progresses by periodic calls to ADDTCR. (See description of ADDTCR in Section 8.2.) 4-1 To prevent the boundaries of adjacent packages from cros- [ sing, the tracers along the common boundary coincide exactly. j Because these identical points are moved with exactly the same velocity, they remain superimposed during the entire " calculation. | To allow for spatially disconnected subpackages, the j code employs the following convention. The last point of a material package or subpackage is a dummy point, its x- coordinate is -1000, its y-coordinate is 0. A given material package can be divided into any number of spatially discon- nected subpackages. 4.1.2 Movement of Tracer Particles Each material package is circumscribed by a series of massless tracer particles, which are propagated each time step a distance'^ Ax^ = u^ ^^INFACE (4.1) where u. is a local, density weighted, average velocity vector for the continuum, determined by an area overlap method which- gives weight to velocities in the surrounding cells. Specifi- cally, a rectangle of cell dimensions is superimposed on the particle to be moved and then u, 2 w(L) (4.2) Since the tracer coordinates are in grid line units, this distance is converted from centimeters to cell units. 4-2 1 I j i < ^ ■f where w(L) is the cell-centered velocity* of the overlapped cell L, is the overlap area and Pj^ is the total cell mass divided by the total cell volume. This procedure gives a spatially continuous velocity field for particle propagation. The density weighting causes the particles to follow the motion of the denser material and prevents the inter- face from becoming unstable in regions where there are large density discontinuities. Infrequently, the sum of the area terms, A^, equals zero, which indicates all of the overlap cells associated with a given tracer particle are empty and it cannot be moved in the usual manner. When this occurs, the code will move the isolated tracer to the position of one of its immediate neighbors. However, if the isolated tracer is on a grid bound- ary, it will not be moved off of that boundary, only along it. If INFACE is subcycled (CYCMX > 1), then the time step. At INFACE' used to move the tracers on a given subcycle, is the time step computed by CDT divided by the number of INFACE subcycles, i.e. At At CDT INFACE CYCMX (4.3) *These cell-centered velocities have been updated by HPHASE (pressure effects) and SPHASE (strength effects) for this time step. 3 4-3 The passive tracers (XP, YP) , used only for tracking the mass within the package boundaries, are moved by the same method, except that if INFACE is subcycled they are moved only on the last subcycle using the time step 4.2 CREATION OF INTERFACE CELLS In the discussion that follows an interface cell is any cell whose boundary is intersected by a material interface, whereas a pure cell is completely within one of the material package boundaries. An interface cell may contain one or more distinct materials, it may contain a void and therefore be a free surface cell, and it may contain a slipline. The characteristics of the various types of interface cells are given in Section 4.5. The way in which interface cells are flagged and the additional information stored for each of them are discussed in the sections below. I 4.2,1 Use of MFLAG Array In all cases, the value of MFLAG (k) of an interface cell k is greater than 100. The value of MFLAG (k) of a pure cell k, however, is less than or equal to the number of material packages in the grid (0, 1, ..., NMAT) . For interface cells, the MFLAG value has a twofold function. First, as indicated above, it indicates that the cell is intersected by an interface. Second, it links the k index of the cell to an index in those arrays which store the mass, density, specific internal energy, and velocity components of each material in the interface cell. Given m = MFLAG(k) - 100, the mass of material package n in cell k is stored in XMASS (n,m) ; similarly the density, specific internal energy, and velocity components of material package n are stored in RIIO (n,m), SIE (n,m), US (n,m) and VS (n,m), respectively. In this way the k index of the cell 4-4 -V arrays are linked to the m index of the material arrays. Furthermore, the cell arrays are defined by summing over the materials in the interface cell as follows: AMX(k) XMASS(n,m) U(k) = yuS(n,m)‘XMASS( / n (4.4) XMASS(n,m) (4.5) V(k) = J]vS(n,m)-XMASS(n,m)/^ XMASS(n,m) (4.6) n / n AlX(k) = ^SIE(n,m) •XMASS(n,m)/j^XMASS(n,m). C4.7) n / n Listed in the table below are the arrays which store cell quantities and the arrays which store material quantities within the interface cells. Throughout this document these will be referred to as cell arrays and material arrays, re- spectively. QUANTITY CELL ARRAY MATERIAL ARRAY Mass AMX XMASS Specific Internal Energy AIX SIE Radial Velocity U US Axial Velocity V VS Density (Not stored) RHO Pressure P (Not stored) Shear Stress STRSRZ (Not stored) Radial Stress Deviator STRSRR (Not stored) Axial Stress Deviator STRSZZ (Not stored) 4.2.2 Definition of Material Arrays Three subroutines, NEWFLG, NEWMIX, and NEWRHO are called to define the material arrays for a cell k that is being con- verted from a pure to an interface cell, which occurs when a material interface initially enters a pure cell. First NEWFLG is called to find an unused location, m, in the material arrays. The convention which denotes that a material array location m is not being used is that RH0(l,m)=-l. Having found an unused location (HELP performs an error exit if no such location is found), NEWFLG redefines the value of MFLAG(k) to be m + 100. In the case when an interface cell becomes a pure cell, which occurs when all material interfaces leave the cell, it is necessary to release the members of the material arrays which have stored the material properties for that interface cell. This process, which is implemented in TPHASE, consists of initializing the particular material variables to zero and setting RH0(l,m) = - 1, which denotes that the material variables whose second index is m are unused and can be assigned to another interface cell by NEWFLG. The process of converting an interface cell to a pure cell is completed by setting the cell's MFLAG value to the package number of the material that now fills the cell, or to zero if the cell has become completely void. NEWMIX, after calling NEWFLG, defines the material arrays for the material which is currently in the cell. If mo was the value of MFLAG (k) on the previous cycle, when cell k was pure, and m is the location in the material arrays assigned to cell k by NEWFLG [m + 100 = MFLAG(k)], then NEIVMIX defines the material variables for the material in the cell as follows (providing cell k was not a pure void cell, i.e., mo f 0) : ,.<I XMASS (mo, m) = AMX(k) SIE (mo, m) = AlX(k) US (mo, m) = U(k) VS (mo, m) - V(k) RHO (mo, m) = AMX(k)yCELL VOLUME . If mo = 0, the cell was previously a void cell and has be- come a free surface cell, since it contains a material interface. In that case RHO (NVOID, m) is set equal to 1, signifying that the cell is a free surface interface cell. Any interface, of course, is a boundary between two materials, or one material and a void. Therefore, when a cell becomes an interface cell, both materials or the void have to be identified. The identity of one material, or void, is indicated by the flag of the cell when it was pure (mo in the example above) . The other material or void is identified when the code senses which package tracer particles are in- tersecting the cell's boundaries. If that other package is the void, then RH0(NV0ID,m) is set to 1. Otherwise sub- routine NEWRHO is called to assign a density to the other material. NEWRHO picks the density from the neighbor cells nearest the intercept of the interface with the cell boundary. This density value is used, in most cases, as a flag in- dicating which other material interface has entered the cell. Since the densities of multimaterial interface cells are re- defined in the pressure iteration (see Section 4.6.1) at the beginning of the next cycle, the value chosen need only be an approximation of the density of the material about to be transported into the cell. Only if the interface cell had been a pure void cell (mo = 0) and will contain only the new material which is entering the cell does the density defined by NEWRHO remain unmodified by the pressure iteration. This is the only case in which the density assigned to the material by NEWRHO is used as a real measure of the density of the material in the interface cell. 4-7 I It is important to note that the material density array, RHO, is used as a flag as well as a material property. RH0(n,m)>0 indicates that the interface of material package n intersects the m interface cell. This flag must be set before any mass of package n can be transported into inter- face cell m. However, it is possible for the density variable, RIIO (n,m), to be non-zero and the corresponding mass variable, XMASS (n,m), to be zero. In that case, the non-zero density only indicates that the interface of material package n is intersecting the cell. On the other hand, when the material n interface leaves the m interface cell, the density variable, RHO (n,m), is set to zero to indicate that all the mass of material n should be "evacuated" from the cell on that cycle. In this case the density variable, RHO (n,m), being zero and the corres- ponding mass variable, XMASS(n,m) , being nonzero signifies that any influxes of that material into the cell should be set to zero and the outfluxes adjusted so that their sum ex- actly equals the remaining mass. This adjustment of the trans- port terms is done by subroutine DMADJ and is discussed in more detail in Section 4.4. 4.2.3 Accounting for Subcycles INFACE usually is subcycled to minimize the transport noise which occurs when an interface leaves a cell. The routine is executed two or more times each cycle using a fraction of the time step on each subcycle. It is therefore possible to create an interface cell after one or more subcycles of INFACE have been completed. In that case the mass transport variables (SAMMP, SAMMY, SGAMC, SAMPY) for that cell have not been accumulated for the completed subcycles and they need to be computed before adding in the transport 4-8 M I terms for the current and subsequent subcycles. Therefore, when INFACE has completed at least one subcycle and storage for the new interface cell has been set up, NEWMIX calls DMCALC, where the computation of transport terms is performed. 4.3 CALCULATION OF CELL FACE AREAS FOR TRANSPI INTERFACE CELL BOUNDARIES The computation of fractional cell boundary areas is performed in order to use the material interface positions to define the transport of a given material across an inter- face cell boundary. For pure cells the mass transport be- tween two cells is V A At, (4.8) where is the donor cell density, v is an average of the cell-centered velocities of the two cells, A is the area of the cell boundary across which the mass is being transported and Atpj,™ is the time step* computed by CDT. For interface cells this equation becomes ''"'n ■ On % % “iNFACE (4.9) where the density, p^, and velocity, v^^, are based on material quantities (RHO, US, VS) rather than cell quantities (AMX, U, V) . The area term, A , is associated with material n and is n determined by the position of the material n interface, while the time step, is the time step for the cycle computed by CDT, divided by the number of passes through subroutine INFACE (see Eq.4.3). Thus the code computes a mass transport term for each material which intersects any boundary of an interface cell. * See DT in Section 10.3, 4-9 This section describes the procedures used to compute the cell face area terms which are used to transport material across the boundaries of interface cells. 4.3.1 Defining Fractional Areas of Intersected Cell Faces The cell face area terms associated with each material in an interface cell are computed in subroutine FRACS. The method given below for computing these areas puts no restrictions on how many interfaces are in a cell nor on how many times a single interface crosses a cell boundary. It does, however, assume that a material interface never intersects itself, and experience has shown that the logic of the code breaks down if it does. This is not a limitation of the method since the boundary will not cross itself unless the tracers become too sparse, and this can be prevented by periodi- cally adding new tracers in subroutine ADDTCR. (See Section 8.2) The method for computing the area terms will be better understood if the following assumptions and conventions are kept in mind. First of all, FRACS computes and stores only the area terms for the right and top boundaries of each in- terface cell, since the cells below and on the left provide the area terms for the other two boundaries. FRACS stores these area terms for interface cells in two arrays, FRACRT and FRACTP, for the right and top cell boundaries, respectively. These arrays are doubly dimensioned^ given FRACRT (n,m) or FRACTP(n,m), the first dimension, n, gives the material pack- age number, and the second dimension, m, is the index that links the interface cell to the material arrays by the relation m = MFLAG(k) - 100. Also, because the area term of a given material package appears in the mass trans- port equation, it must be defined for the right and top boundaries of every interface cell, even for those boundaries not intersected by a material interface. 4-10 Let us first consider how FRACS computes the area terms for the cell boundaries that are intersected by a material interface. The set of tracer particles which circumscribe a given material are considered consecutively (two at a time) by FRACS. (TX and TY are the FORTRAN variable names for the particles' x and y coordinates, respectively. The TX and TY arrays are doubly dimensioned,; the first dimension gives the mat- erial package number, the second gives the sequential order- ing of the particles. For simplicity the coordinates of the tracers are in grid line units rather than centimeters, since in grid line units the integer part of the coordinates in- dicates whether a pair of tracers straddles one or more cell boundaries.) If a pair of material n tracers straddles a cell bound ary, FRACS computes the intercept of that boundary with the line between the two tracers. Given that point of intersection FRACS defines the area term for material n at that cell boundary to be the fractional cell boundary area on one side of the intercept. To determine which side of the intercept to associate with material n, FRACS uses the convention that material n is on the left as one proceeds from the first to the second tracer of the pair. For the cell boundary not intersected by a material interface the area term for a given material is either the total cell boundary area or is zero. The following discus- sion illustrates how the area terms of these non- intersected cell boundaries are correctly defined by "presetting" them when an interface enters a cell and "resetting" them when it leaves. 4.3.2 Presetting Non- Intersected Cell Face Area s The sequential order of the tracer particles determines whether an interface is entering a cell or leaving it. When the interface is entering a cell, and is crossing this cell boundary for the first time, the program processes the other boundaries of the cell in a clockwise order, presetting the corresponding area terms to the total area of that cell bound- ary, until it encounters a side that has already been crossed by the same material interface (i.e., a side which has an associated area term that is non- zero and is less than the total area). The first three c the presetting procedure. Case (1) ases diagrammed below illustrate Area terms preset: Bottom = IT - X? . (y. - y._j) zn Top » ifx? - Xj_j] Area terms preset: FottOm = TT |x^ - i-1 1 Case (3) Area terms preset: None - the bottom boundary has been intersected previously by the same interface. In order that there be no restriction on the number of times the material interface crosses a given cell boundary, the presetting procedure is generalized in the following manner. If the interface of material n crosses a cell boundary for the second time, the code does not preset the area terms of the other sides of the cell if material n lies between the two points of intersection. This is illustrated by Case 4. Case 4 Area terms preset: None - material n lies between the two intersection points. The pro , n senses that material n is between the points when 1 of the area term computed on the previous cross, as) and the one coinputed currently is greater than the total area of thn cell face. If, as in Case 5, the Kuni of the area term computed on the previous crossing (s) '’nd the one computed currently is less than the total area, material n is outside the two intercepts and the procedure described above for presetting the area terms of the other sides of the cell is followed. It should be noted that these conventions apply re- gardless of how many times the interface crosses the boundary, as illustrated by Case (6) and Case (7). Case (6) Area terms preset: / 2 2 \ Bottom = Left = (yj - Case (7) f/XrC di: Previous - X - — Current Area terms preset: None vrr. t 4.3.3 Resetting Non- Intersected Cell Face Areas The procedure for resetting the area terms when the in- terface leaves a cell is the following. If the interface is crossing the cell boundary for the first time, the program processes the other faces of the cell in a clockwise order and resets their area terms to zero, until it encounters a side that has been crossed previously by that same material interface. This procedure is illustrated by Cases (8), (9), and (10) below. Case (8) Area terms reset; Top ■ 0 Right = 0 4-15 Case (9) Area terms reset; Top = 0 Casa (10) Area terms reset None In Cases (11) and (12) below, the side where the inter- face leaves has been previously crossed by this interface. The program sums the area terms resulting from any previous crossings of that boundary and the current crossing. If, as in Case (11), the sum of these areas is less than the total cell boundary area, the material is outside the two inter- cepts and none of the area terms of the other sides are set to zero. 4-16 Case (11) Area terms reset; None , Previous ft#r/al n 'Current If, however, the sum is greater, as in Case (12), the material is between the intercepts and the program proceeds to reset to zero the areas of the other non- intersected sides of the cell . Case (12) Previous Area terms reset; Current xop = 0 Right = 0 Bottom = 0 4.4 ADJUSTMENT OF INTERFACE CELL MASS The four maS transport terms for interface ceils (one for each cell face) We computed by DMCALC using material densities and velocities and the cell face area terms associated with the material packages as computed by subroutine FRACS. However, sub- sequent to this calculation, DMADJ is called to (possibly) modify these transport terms. There are two reasons the 4-17 transport terms may need to be adjusted: one, to exactly empty a cell of a material when its interface leaves the cell, i.e., evacuation; and two, to prevent the outflux of material from exceeding the mass of the material in the cell, i.e., over- emptying. 4.4.1 Evacuation of a Material To signify that the material n interface has passed out of cell k the density of material n, RMOCn,m), is set to zero. This zero density is used as a signal to subroutine DMADJ to adjust the transport terms (one for each cell face) of material n so that cell k will be exactly evacuated of material n. This adjustment of the four transport terms usually occurs twice in the computational cycle; once at the beginning of the first INFACE subcycle (in case the material r interface left the cell on the previous cycle*) and once after all INFACE subcycles are completed (in case the mate- rial n interface moved out of cell k during or after the first, and before the last, INFACE subcycle). For the first case the four transport terms from the previous cycle are used to indicate the direction of the flow and the direction of the evacuation. In the second case, the four transport terms just computed by DMCALC are used to determine how material n is to be evacuated. If material n is to be evacuated from cell k, then all influx transport terms of material n are set to zero, and each outflux term is adjusted in proportion to its fraction of the total outflux computed by DMCALC. If the total out- flux computed by DMCALC is zero, the material is removed from the grid (evaporated) at the end of TPHASE. *The tracer particles are moved at the end of each INFACE sub- cycle. If the interface moves out of a cell on the last sub- cycle of INFACE, the code does not sense this until the first pass through INFACE on the next computational cycle. 4-18 rniMmiMimm nM^ These evacuation procedures are used for a cell that has lost one interface, but is still cut by another inter- face, as well as for cells that have become pure. After the area terms associated with each material package have been computed, FLGSET searches for cells that were interface cells on the previous cycle but are no longer cut by an interface. Such a cell has become pure and its flag, MFLAG(k), is made negative to signify that fact until the transport has been completed in TPHASE. To determine what material package a newly pure cell k belongs to, FLGSET first looks for a neighbor that is pure and assumes cell k will belong to the same package as one of its pure neighbors. If, however, all of its neighbors are interface cells, the top cell face area terms of the cell below? kb, are examined. If the cell face area for material n at the top of cell kb is equal to the total cell face area, then cell k must be in- side the material package n boundary. Therefore, all material densities for cell k except that for material n are set to zero, and at the end of TPIIASE cell k becomes a pure material package n cell with MFLAG(k) = n. 4.4.2 Overemptying of a Material After the subcycles of INFACE are completed, subroutine DMADJ is called to check that the outflux of each material in each interface cell does not exceed its current mass. If the sum of the outgoing transport terms of material n is greater than the current mass of material n, the outgoing transport terms of material n are reduced proportionately, as in the evacuation procedure described above. *The cell on the left of cell k is used analogously if cell k is in the bottom row of the grid. 4-19 4.5 There are several types of interface cells. In the following paragraphs they will be defined and then discussed in terms of special prescriptions that apply to them. 4.5.1 Multimaterial Cells Most interface cells contain more than one material, and the only practical restriction on the number of materials in any one cell is the array dimensions imposed by the user? A detailed description of the determination of material properties in multimaterial cells is given in Section 4.6. The following is a summary. The partial volumes of the materials in a multimaterial cell are determined when material pressures are equilibrated. The densities of the materials are adjusted in an iteration until the pressures of the materials equilibrate. (See Section 4.6.1) The shear yield strength of a multimaterial cell is a volume weighted average of the shear yield strengths of its constituents. The code assumes, however, that a multimaterial cell has no tensile strength, i.e., negative pressures are set to zero in multimaterial cells. (See Section 4.6.2) The internal energy increment of a multimaterial cell is partitioned by partial masses in SPHASE and HPHASE. The internal energy is updated in TPHASE to account separately for the transport of internal energy and the thermalization of kinetic energy of each material in the cell. (See Section 4.6.3) 4.5.2 Free Surface Cells Any cell that is intersected by the void interface is regarded as a free surface cell, and is flagged via the RHO array: RHO (NVOID,m) = 1, where NVOID is the number of material * The MFLAG convention, however, imposes a limit of 99 material packages in a given problem. 4-20 packages plus one (NMAT+1), and m = MFLAG(k) -100. A free surface cell can contain more than one material as well as the void, although frequently free surface cells do contain only one material. In CDT, if the equilibrated pressure of a multimaterial, free surface cell is less than PMIN, a user-controlled in- put number (See Section 7.2.1), the densities of the materials in the cell are not modified unless the sum of their partial volumes exceeds the volume of the cell. This rule prevents material densities in a free surface cell from being based on the mass being extended over the entire cell volume. How- ever, if the free surface cell has only one material, the pressure iteration is bypassed and the density, used in the equation of state, RHOW, is the total cell mass divided by the total cell volume. If the cell pressure based on this density is negative, it is set to zero and the density of the material, RHO, is redefined only if it is less than the assumed density, RHOW. (See also Section 4.2.2.) In SPHASE, the free surface cells are assumed to have no strength. Therefore the deviator and shear stresses of a free surface cell are zero. In DMCALC, when computing the mass transported from a free surface cell into a non-free surface cell, the acceptor cell density instead of the donor cell density is used, since, in general, the densities of the non-free surface cells are better determined. 4.5 . ?> Slipline Cells Sliplines in HELP are not restricted to grid lines or cell diagonals. They do coincide with material interfaces, however. Therefore a cell that contains a slipline is also an interface cell. 4-21 The materials in a slipline cell are accelerated in inverse proportion to their densities in HPHASE. (See Section V.) In SPHASE the assumption is made that a slip- line cell has no strengtli, therefore its deviator and shear stresses are set to zero. The velocity components of the "slave" and "master" materials in a slipline cell are different. The cell velocity components represent a mass weighted average of the material velocity components. These cell velocities are used in computing the work terms in HPHASE and SPHASE and in moving the material and passive tracers in MOVTCR. The material velocities are used to compute the mass and momentum transport terms in TPHASE. 4.6 DETERMINATION OF MATERIAL PROPERTIES IN MULTIMATERIAL CELLS Pressure and Density Iteration for Multimaterial Cells For one-material cells, it is possible to compute the pressure directly from the material density and specific internal energy. However, for cells containing more than one material, it is necessary to make additional assumptions or observations in order to determine the cell pressure. In the present method we determine the cell pressure for multimaterial cells by an iteration procedure. Specifically, the densities of the various materials within the cell are varied, subject to the condition that the cell be exactly filled by the masses therein, until the individual pressures converge to a common value, which is taken to be the cell pressure. Specific internal energies are held constant during the iteration. This process gives, in addition to the cell pressure, the densities of the individual materials for subsequent use in the DMCALC and TPHASE transport calculations. The procedure for determining the pressure and the densities in a multimaterial cell is actually divided into 4-22 two parts: a pre- iteration calculation to insure that the cell is filled exactly; and an iteration calculation to con- verge on the desired P*s and p's. 4. 6. 1.1 Pre- iteration Calculation In general, the masses and densities which exist at the beginning of the calculation do not exactly fill the cell, since cell masses are changed in the TPHASE calculation of the previous cycle. The pre- iteration calculation fills the cell exactly, taking into account the (computed) compressi- bilities of the various materials in the necessary alteration of densities. Stepwise, the code 1. Computes P^ from the old value of for each material within the cell (E^ is assumed constant throughout) . 2. Varies the by one percent to compute new P^. 3. Determines the constant-energy compres- 2 sibilities, called for each material C? = AP./Ap. (4.10) where the AP^ and Ap^ are given from the results of (1) and (2). 4. Normalizes (see proof below) to fill cell exactly: AV^ = -AP/h? (4.11) where 4-23 if i i li ’ I t- It- = ( Pj Ci ) (4.12) and AP (VOL-VCELL) (4.13) VOL (4.14) with "= mass of constituent i, and VCELL = volume of the cell. With the above increments in specific volume, the new specific volumes are given by " V. , + AV. 1 i,old 1 (4.15) Then, Pi . l/Vi (4.16) gives the desired nev; densities which cause the cell to be filled exactly. It is easy to see that the calculation defined above causes the cell to be filled exactly. To show this. 4-24 Summing , ’’i/Pi.old - "'i m. (VOL-VCELL) m . / D • n J “ — i' ,old = VOL ■E m , >'1 - (VOL-VCELL) (4.17) = VCELL (4.18) showing tliat the new volum^js m^V^ of the constituents add up to the total cell volume, as desired. 4 . 6 . 1 . 2 Iteration Calculation Given P. = £■ (v. , E.) for materials i= 1,...,M, the cell 1 1 \ 1 1 / pressure P = P- = ?2 = ... = Pj^ is found by varying the V^ until the P^ are equal within some specified accuracy, subject to the conditions that the cell remain exactly filled and that the are constants. The equation for volume conservation is J^m^ AV^ = 0 (4.19) and the equation that causes material i to undergo a change in specific volume AV^ , such that its pressure changes from 4-25 its current value = ^i(^i’^i) ^ value common to all materials in the cell at the end of the iteration cycle, is (4.20) These equations, for i = 1,...,M can be regarded as M+1 equations for the M+1 quantities AV^ and F" . The other quantities in the equations are either known constants (m^) or are updated each cycle of the iteration [the and the |3V^/9P^j ] and are taken to be constants while solving ^i — ^ . *1 these equations for AV^ and P . The solutions are where rn+1 AV. 1 p. w” LmJ 1 1 (p? - 1 1 FT 1 rn+1 (4.21) (4.22) (4.23) w. = m./h. 1 1 1 (4.24) That these equations are solutions to the given equations can be verified by direct substitution. The iteration would be exact (single cycle convergence) if the input coefficients (c'P./9V.) were constants, since the solution does not involve additonal approximations. In the iteration, these coefficients are refined each step of the iteration by making use of the P^, points which were generated during the previous step. Further, some pains are taken to^ determine realistic starting values of these coefficients, in the pre- iteration calculation, in the belief that this procedure will hasten convergence sufficiently to reduce total computing time. Stepwise, the iteration calculation is as follows: 1. Compute new and AV^ by the above formulas, using known values of and (av,/aP,)^^. n 4» 1 2. Compute new v” by adding the above AV^ to the old specific volumes. 3. Compute new P^ , using the updated V^. 4. Compute refined values of the coefficients |3V^/3P^j by using this last point [from (2) and (3) above] and that previous point (two such points are saved for each material during the iteration) which is The iteration is completed when the all equal P to some preset accuracy, PRCNT (See Section 7.2.1). In most exper- ience to date, convergence has been obtained in two or three steps, where a convergence criterion |P^ - F| < 10 P was used. 4. 6. 1.3 Equation o£ State Modifications for the Iteration roceoure The Tillotson equation of state (described in Section 2. 3. 1.2) has been used in most work to date with HELP and is sufficiently general to provide a satisfactory description of most media which is initially condensed and which may vaporize as a result of heating. However, for some materials which are sufficiently expanded and cold (E < E^), the slope of a con- stant energy curve (in the P-V plane) is positive. For pure cells, these meaningless states are avoided by replacing the 2 quantity (Ap + Bp ) with its minimum value (a negative pressure) if the cell density is less than the density corresponding to this minimum value. To assure convergence of the pressure iteration in multimaterial cells, it is necessary to have a continuous equation of state and for 8P/8V to be everywhere negative. To this end the term Ap + Bp‘ (4.26) is replaced by Ap + h(v)Bp‘ where 1. h(v) ■ +1 for V £ Vjj (compressed states) 2. h(v) decreases linearly to -1 in the interval ^0 < V < V^ ^ (4.28) 3. h(v) - -1 for V > V - ’'o ina The effect of these modifications is to change B to -B in a continuous manner and hence insure that 8P/3V is always nega- tive. Negative pressures that may result after the iteration is completed are set to zero, since tensile states in inter- face cells are assumed not to exist. In view of this latter assumption, the indicated equation of state changes do not normally affect the physical result, since they alter only states that would ordinarily be tensile. The purpose of the changes is to make the equation of state well-behaved, so that the iteration converges properly. 4.6.2 Shear Yield Strength of Multimaterial Cells The shear yield strength of a multimaterial cell is a volume weighted average, 7, of the shear yield strengths- of its constituents, Y„: n (4.29) If, however, one of the constituents is an ideal gas or a high explosive, (MAT(n) ^ 20), the shear yield strength of 2 .. 4-29 wmmmm the cell is set to zero. Likewise, if the cell contains the void interface (RHO(NVOID,m) ■ 1) or a slipline interface (THETA (m) >0), with m ■ MFLAG(k) - 100, the shear yield strength of the cell is set to zero. And finally, the shear yield strength of the cell is zero if the density of any one of its constituents is below the density corresponding to the maximum tensile strength of the constituent (RHO(n,m) <p *AMDM(n)). This presumes that a material that n has failed in tension has no shear yield strength. 4.6.3 Partitioning of Internal Energy Increments in Multimaterial Cells The internal energy increments, AEj, computed by HPHASF and by SPHASE are partitioned amongst the materials in a multimaterial cell in proportion to their masses. For ex- \ ample, the internal energy increment for material n is AEr »AE 'I n I ■vS rtii (4.30) Iii^ other words the HPHASE arid SPHASE gain in specific internal energy is the same for all the materials in the cell. The situation is somewhat more complex in TPHASE, since the iftternal energy increment is due not only to transport but also to the thermal ization of kinetic energy. It is natural to imagine that the interacting materials (all of the masses of material n within the cell at the end of the time step) have collided in a coordinate system in which their center of mass has zero velocity. First, we recall that the new velocity components for , material n are chbsen to conserve momentum 4-30 - 1,,^) . ^ /{<^ * ^«»l„) {♦•31) C" ' ? ♦M’n)] /("?’ * ?*"in) (♦•”) and that the new mass of material n is ■ bI|P + iC ^™4n n n j" in (4.33) Here the subscript i denotes the four sides of the cell; the subscript n denotes the material package, and the superscript (n) denotes pre-TPHASE quantities. The velocity associated with each of the transport masses is the donor cell material velocity; e.g., if a mass Am^^ entered the cell under consider- ation through side 1, the associated velocity u^^^ is the velocity of material n in the cell from whence the mass came. We note first that as determined to con- n ' n serve momentum by Eqs. (4.31) and (4.32) above, are the velocities of the center of mass of the material n in the cell at the end of the time step. This means that a component mass with pre-TPHASE velocities, following amount of kinetic energy thermalized ,(n+l) . ^(n) in the center of mass coordinate system. Thus the total ther- malized KE for material n is 4-31 •••' Tj) \ --r r,' ■■ ■ 7 (n) . ,,(n+l) Y * (vi"’ - vr^’)1 ,6niin ,(n+l) (n+1) Some interpretation of this equation clarifies its meaning: the first term, plus those terms from the summation for which 6min < 0 (material leaving the cell with donor cell velocities represents the ther- in n ’ in n malized KE for that fraction of the material which remains in the cell. The additional terms (6m. > 0 in the summation) in represent the thermalized KE associated with that material which came into the cell during the time step. For these terms the uf^^ v|”^ refer to the velocities of material in in in the neighbor cells which acted as donors. The change in internal energies is again given by EaEii„ . TKE„ (4.35) In other words, material n gains an amount of total internal energy which is that due to transport, plus the amount of thermalized kinetic energy for material n. The above discussion assumes that the material velocities are not equal to the cell velocities, i . e. , that the cell con- tains the slipline. If the cell does not contain a slipline, equations (4.33) through (4.35) still hold, only the material velocities are defined using cell quantities as follows: 4-32 CHAPTER V SLIPLINES S . 1 GENERAL COMMENTS Sliplines have existed in Lagrangian codes for several years ani some of our thinking and terminology comes from these earlier accomplishments. For example, in HELP one package, is considered the "master" along which another "slave" package slips. However, an important difference between this current effort using an liulerian code and the earlier Lagrangian models is that the sluipe of the slip surface is determined by th^ stress field and need not coincide with grid lines or node points. The slipline in HELP is also a material package inter- face; therefore, the slipline cells arc also interface cells. The angle of the slipline in a given slipline cell is measured from the x-axis of the grid, using *the intercepts of the slip- line interface with the cell boundaries. Tiiese intercepts are computed while processing the interfaces in TRACS and are stored by PTSAV. The stored intercepts are then used by THETAS to compute the angle of the slipline in each slipline cell. Since the slave and master materials in a slip cell develop different velocity components tangent to the slipline, it is necessary to store distinct velocity components for each material. (In an earlier version of HELP, it was assumed that all materials in an interface celJ would have the same velocity components.) Two material arrays, US, \^S, have been added to store the material velocities in interface cells. To describe the angle of the slipline across a given slip cell, the THETA array was added. When THETA (m)=-l (m = MFLAG(k) 100), the interface cell does not contain the slipline; when 5-1 r*i WpiWiiwi i ’irnYwmtfU'Wt. WT A ’IWt »l* I'smiO'IW W ’ BNl^f WWKf *T* f I T[!nTA(m) > 0 the slipline docs cross thr- interface cell at i an anslc (in radians) equal to THETA (m). I It should be noted that more than one material package I • can be designated as a slave or master. Therefore it is I possible for a sliplinc cell to contain more than one slave i ~ or master package. In that case, the velocities of all the • slave materials in the cell arc equilibrated as are the I velocities of all the master materials. The slave density, , Pg, used in the equations that follow is an average of the densities of all the slave materials in a given slip cell, likewise for the master density, 5.2 RELATIVE ACCELERATION OF MASTER AND SLAVE MATERIALS rUVCALO A cell that contains the sliplinc has, in a sense, failed. Therefore the code assumes the material in a slip cell has no shear yield strength and the deviator stresses of a slip cell are set to zero. The only stresses acting on a slip cell are simply the hydrostatic pressures. The acceleration of slip cells due to stresses is computed entirely in HPHASE. After updating the momenta of a slipline cell, IIPHASE calls UVCALC to establish the relative acceleration of the master and slave materials. Given the post-MPMASE cell centered momenta (m(k)u(k), m(k)v(k)) and the angle 9 of the slipline through cell k, the master and slave material velocities, uSj„» vSjj^, and us^, vs^., are updated by solving the following four simultaneous linear equations: vs^ + m„vs = (m^ + m„) v(k) (5.1) ssmm sm"^^' ni.. + m„us„ = (m^ + m„) u(k) (5.2) ssmm^sm'^^'^ ^ ' 5-2 s in 0 C5.3) uSg sin 0 - vSg cos 0 = u % vs cos 0 m Ps [(uSy-us^) cos 0 + (vs^ - vs^) sin 0] ( 5 . 4 ) where usl , vs'., us', vs' arc pre-llPHASl: slave and master material velocities. The updated material velocities must conserve the cell's axial and radial momenta (nqs. (5.1) and (5.2)). The m.aster and slave velocity components normal to the slipline must be eciiial fP.q. (5.3)), and the HI’HASE change to the slave and master velocity components tangent to the slipline must he inversely proportional to the slave and master material densities, p^., (Eq. (5.4)). Solving these four equations results in the thermalization of some of the slip cell's kinetic energy. This thermalized kinetic energy is converted to internal energy and is partitioned amongst the materials in the cell in proportion to their masses. 5.3 TRANSPORT OF MASTER AND SLAVE NIATERIALS The mass, momentum and^ energy transport terms are com- puted in DMCALC and TPHASE. Since some interface cells, namely those along the slipline, have different velocity components for the master and slave materials, the mass and momentum transport equations for interface cells must employ material velocities rather than cell velocities. For example, the mass of material n to be transported across the i cell face is 6mn At (5.5) 5-3 ♦ where p^ is the density of material n in the donor cell, u5 is the mass transport velocity based on the material n veloc- ities of the donor and acceptor cells, and is the i cell ^i face area associated with material n. Likewise, the momentum transport terms are based on material velocities rather than cell-centered velocities: «mun “ % * ^mVn - 6m^. vs^ . C5.7) Since the material velocities equal the cell velocities in cells that do not contain a slipline, the equations a re result in the same rnasc and momentum transport for non- slip interface cells as do equations that use cell-centered quan- tities . 5.4 MODIFICATION OF POST-TPHASE M^VSTER AND SLAVE VELOCITIES rUVMODl After the material velocities of a slip cell have been updated in TPHASE, they must be altered to insure that the material velocity components normal to the slipline are equal. Subroutine UVMOD is called from TPHASE after all cells have been updated for the effects of transport. UVMOD searches for interface cells that contain a slipline (THETA(m) > 0) and modifies the post-TPHASE master and slave material velocities of those cells so that their components normal to the slipline are equal, i.e., so that Eq. (5.3) in the discussion above is satisfied. The material velocities of non-slip interface cells are not modified by UVMOD. 5-4 5.5 USE OF SLIPLINES Sliplines in HELP have been employed in a variety of situations; for example, at the high explosive-metal inter- faces of shaped charge and fragmenting munitions calculations, and at the interface between the plug and target when modeling thin plate perforation. When there is a large density difference between the master and slave materials (as between an expanded HE and a metal casing) , the relative motion of the two packages is generated primarily by the HPHASE accelorations> which are density weighted and which satisfy Eqs. (5.1) through (5.4) in the discussion above. However, when the densities of the master and slave materials are close (as between a plug and a target), the relative motion of the packages results primarily from the different momentum transported to the master and slave materials in TPHASE. For example, given that the pro- jectile and plug are both slave packages, and the target is a master package, the material velocities of the plug and projectile in a slip cell are equilibrated, thereby trans- ferring the projectile momentum to the plug but not to the target. 5-5 WSW’KWIW*'''' iYVt'rVftViHWT ■’.nvviKT'X'^^T •W<lW'W«Wnw»<W'i'»»rr»^«W:W CflAPTER VI PLUGGING FAILURE A plugging model has been incorporated into the HELP code which the user activates if he or she anticipates plug- ging failure. As the following discussion will clarify, the decision to use the plugging failure model must be made when the calculation is generated. As the model now exists, its application is restricted to the impact of right circular cylinders or truncated cones into homogeneous targets of finite thickness. Furthermore, the user must supply an empirically derived material parameter (PLWMIN) which deter- mines the shape and mass of the plug that is formed In the sections that follow, the procedures for generating a plugging calculation as well as the code’s method for evolving the plug will be discussed. 6.1 GENERATION OF A PLUGGING CALCULATION 6.1.1 Designation of Material Packages The plugging model assumes the projectile is package 1, the plug is package 2, and the target is package 3. The pro- jectile can be made up of several m.itcTial packages as long as the part of the projectile which will interact with the target is homogeneous and is designated as package 1. Further- more, the projectile and plug packages must be "slave" packages and the target a "master" package. 6.1.2 Placement of Plug Tracer Particles Unlike the other material packages in a HELP calculation, the dimensions of the plug package are unknown at the beginning of the calculation. However, the user does specify in the input deck a set of plug tracer particles which later in the calculation are automatically redefined so as to circumscribe the plug package as it evolves. Initially these plug tracers coincide with the projectile and target tracers along the projectile-target interface. These initial tracer positions are based on the assumption that the vertical edge of the plug package begins at the top right edge of the projectile. In the example below, the relationship between the pro- jectile, target and plug tracer particles at time T = 0 is illustrated : Projectile - The last six projectile tracer particles (excluding the dummy end point*) are at positions F through A. * The last tracer of a package is a dummy point. It is used only to indicate the package has been completely circum- scribed. 6-2 I- '* Target - The first six target tracer particles are at positions A through F. Plug - The first six plug tracer particles are at positions A through F. The final six plug tracer particles (excluding the dummy point) are at positions F through A. w. There is a set of n duplicate plug and target tracer particles placed at position F in order to define the vertical edge of the plug package as it is extended. (See Section 6.2. 2.1.) These n duplicate particles are generated automatically by PLGADD, given the user- supplied definition of the plug slipline end- points (nbGSI)( 2), NENDSD(2)j. The number of duplicate particles generated is (NENDSD(2) - NBGSD(2) + 1) and is therefore controlled by the user. 6.1.3 Placement of Passive Tracer Particles The user must specify a region of the target in which the plug is expected to form. In each cell of that region, four passive tracer particles are automatically placed by SETUPA. These passive particles are used to track material within the target , to accumulate its plastic work, and to compute its angle of maximum shearing stress. This information is needed for the plugging failure criterion discussed in Section 6.2.1. 6.1.4 Slipline Specification To simulate the action of a plug being pushed out of a target plate, it is necessary to generate slip surfaces be- tween the plug and the target and between the deforming projectile head and target lip. These slip surfaces do not exist, of course, until the plug evolves. In a plugging calculation the tracer particles that define the endpoints of the slip- line for each package are all initially at the top right corner of the projectile, where the plug begins to form. 6.2 METHOD OF EVOLVING THE PLUG To simulate plugging failure, the HELP code essentially evolves a plug package the shape and mass of which are deter- mined by a plastic work failure criterion and by the assoc- iated angle of maximum shearing stress. 6.2.1 Plugging Failure Criterion The failure criterion used to propagate the vertical edge of the plug is based on the specific plastic work of the target material in the region near the endpoint of the plug edge. The plastic work of this material is approximated by using passive tracer points that move with the target material and by incrementing, by an area weighted average, the plastic work of the material associated with each point. Also assoc- iated with each point on a given cycle is an area weighted average angle of maximum shearing stress. Once the plastic work associated with one or more points in the region of the endpoint exceeds a predetermined cutoff value (PLWMIN) , the average, current angle of maximum shearing stress associated with those points is used to determine the direction in which the vertical edge of the plug is extended. Each cycle in SPHASE, a plastic work increment and an angle of maximum shearing stress are associated with the center of cells in the target. An area weighted average of these cell-centered quantities is used to define an angle of maximum shearing stress and to update the total plastic work 6-4 T*»» t T<»g- v « w grrvaW’^ for the material associated with each moving point. For example, consider point j in the figure below. 1 ^ ■ “aI 1 1 1 1 1 •j 1 \ 1 . _ 1_ _1 The plastic work, W, associated with point j at time n is defined as w- - ( 6 . 1 ) The average angle of maximum shearing stress, a, for the point is defined as: “j ’ • ( 6 . 2 ) 1 In both cases A. denotes the overlap area of the pseudo-cell, 1 11 n • centered at point j, with cell i. AW^ and are the time n plastic work increment and angle of maximum shearing stress, respectively, associated with the center of cell i. Each cycle, after the plastic work associated with the moving points has been updated, the code considers a region which, as illustrated below, is one cell in area and lies above the intercept between the topmost horizontal grid line and the plug's vertical surface. The radial position of the region is centered about this intercept. If one or more of the moving points which lie within that region has an accumulated plastic work total greatei than the specified empirical value, PLWMIN, then the vertical edge of the plug is extended to the next horizontal grid line. The angle at which the plug edge is extended is an average of the current angles of maximum shearing stress associated with those points within the specified region that have exceeded the plastic work cutoff. For example, in the case illustrated above, the code would consider the plastic work associated with passive tracer particles c, d, and f, only, since passive tracers a, b and e are outside the specified region. If, on a given cycle, the plastic work associated with all three points c, d and f exceeds the specified value, PLWMIN, the angle, a, at which the plug will be extended will be a = ^cl ^ “f^ ' point c satisfied the failure criterion, then u = a^. 6-6 I I When the criterion for extending the plug surface is met, the plug package itself must be expanded. The mechanisms for this expansion are described in the following sections. 6.2.2 Conventions Governing Plug Tracer Particles There are three stages of the plug formation. The con- ventions governing the initial placement of the plug tracer particles have been given in Section 6.2.1. The conventions for the intermediate and final stages are given below. 6. 2. 2.1 Partially Formed Plug Once the plug edge is extended, the projectile, target, and plug tracer particles follow the conventions illustrated by the example below, in which 5 of the set of n duplicate plug package tracer particles generated initially by PLGADD (see Section 6.1.2) now define the vertical edge of the plug. I Proj ectile - The last six projectile tracers (excluding the dummy tracer) are at positions F through A. Target - The first six target tracers are at P through K. The next (n-5) are also at K, where n is the number of duplicate points initially speci- fied by NBGSD(2) and NENDSD(2). The next five are at J through I'. Plug - The first six plug tracers are at A through F. The next four are at G through J. The next (n-5) are at K, and the last five (excluding the dummy tracer) are at L through P. Each time the plug is extended all but one of the dupli- cate particles is moved to the new endpoint of the vertical edge of the plug. The one particle left behind helps to de- fine the vertical plug interface. It should be noted here that the interface in the example above does not represent a discontinuity in the material as does the slipline interface KF. Since FF is a psuedo- interface, the material flow field is unaffected by its presence. 6. 2. 2. 2 Completely Formed Plug Once the plug edge is extended to the free surface back of the target and the projectile and target lips coincide, the projectile, target and plug tracer particles follow the con- ventions outlined in the example below. «.>**>*"* .If ««■— w«»V»Tiy»wr» - »' V ■wi'iw f Plun Target iPro j cct I le Projectile - The last six projectile tracers (excluding the dummy tracer) are at positions F through A. The target and projectile have joined between F and F' and the slipline extends from P to F'. Target - The target package no longer connects with the axis. The first eleven tracers are at P through F. The last target tracer (excluding the dummy tracer) is also at P in order to completely circumscribe the target package. Plug - The first six plug tracers are at A through F. The next ten are at G through P, and the last four (excluding the dummy tracer) are at P through S. Note that there are two duplicate plug points at P; one is moved with plug velocities and defines the top right corner of the plug as it is pushed out; the other is moved with target velocities and defines the top left corner of the target. The last three tracers of the plug are coincident with what were the last three tracers of the target just before the plug was completely formed. 6.2.3 Movement of Plug Tracer Particles As described in Section 4.1.2, tracer particles are moved with a density and area-weighted average velocity. However, to better simulate plugging and penetration mechanisms this method has been specialized for those plug, projectile and target tracer particles that define the plug package once it has begun to evolve. 1. The tracer particles along the projectile- plug interface are not moved with target material velocities, only plug and pro- jectile velocities. 2. The tracers along the slipline edge of the plug, except for the top endpoint, are also moved with only plug and projectile veloci- ties. 3. The tracers along the top surface of the plug, except the top endpoint, are moved normally using all velocities in the overlap region. 4. The top endpoint of the plug surface is moved with only target velocities. The following example illustrates the way in which a tracer in category 1 above is moved. If the area associated (2) \lARdliT 1.1) 1^1 1 kv\ kd ppolncTili; (1 with one of these tracers overlaps a cell which contains only target material, e.g., cell kd in the figure, the overlapped area and its associated velocity component are ignored in the de- finition of the tracer's velocity components. The axial velocity of the circled point, therefore, is % ' (^k'^VSi ♦ vS2)/2 * . A^|,-(vsj . VSj)/2)/(A^ * A|^„ * A^^^) (6.3) '4 where is the overlap area between cell i and the pseudo- cell, and vs^ is the axial velocity of material j in cell i. Notice that the target velocity, vs^, in the cells was not in- cluded in the average. The cell centered velocities are used in the averages only when the overlap cell is pure. The reasoning behind these special procedures is that the slipline rep- resents a discontinuity in the velocity field and a tracer on one side of this line should not be affected by the motion of the material on the other side. The tracers along the free surface of the target are moved with target velocities only. Tracers along the free surface of the plug are moved with plug velocities only. The 6-11 ! 1 tracers along the free surface of the projectile are moved with both plug and projectile velocities in the case where the pseudo cell overlaps an interface cell which contains plug material. Otherwise they are moved only with projectile material velocities. Once the plug surface has extended to the rear free surface of the target, the top endpoint of the slipline is considered a part of the target free surface; therefore, it is moved with target material velocities only. As the plug is pushed out, the slipline tracers below the top endpoint will, in general, have a larger axial component of velocity than the endpoint since they are moved with plug and pro- jectile material velocities. The plug slipline tracers are accumulated at the endpoint. The corresponding target slip- line tracers are removed and the number of target tracers is reduced accordingly. 6.2.4 Conversion of Target Material into Plug Materia l 6. 2.4.1 Redefinition of Tracer Particles When the plug surface is extended, the top edge of the plug is redefined, and the plug package boundary encompasses a larger volume as illustrated below. rARtlJ jPosition f orn er Position I'l.l G PRO.IBGTIU 6-12 The tracers at the top of the plug are moved a fraction above the next horizontal grid line, which is also the new y-cnordinate of the slipline endpoint. However, as illus- trated below, those tracers at the top edge of the plug that arc already above the new slipline endpoint are not moved. TAR( ;lt T icse not racei epos i s a re t iom d IVP ‘flSss' ♦ »- PRO JLCTI M V A \ Until the plug is completely formed, the upper endpoint of the slipline is on the top horizontal grid line of the last failed row. However, the top and vertical edges of the plug may extend above this horizontal line, as shown below, where • denotes the slipline endpoints. Last Failed Row TA” GLT PLUG \ --- X PROJI CTILP However, when the plug edge is extended into a cell that con- tains the free surface, the rear of the target has been reached, and the plug edge is extended beyond the grid line until it intersects the free surface, as illustrated below. The tracers at the top of the plug are not moved in this case. Instead, the target tracers to the left of the new plug surface endpoint replace the tracers at the top edge of the plug. At this time in the calculation the number of plug and target tracers is reduced. All but two of the re- peated plug tracers at the plug surface endpoint are removed. One of these will define the top right corner of the plug as it is pushed out; the other will define the top left corner of the target and will also be the endpoint of the slipline. The target package is no longer attached to the axis, and the first and last target tracers define the new slipline endpoint. Accordingly, two free surface tracers are added at the slipline endpoint; one will move with the plug, the other will move with the target and remain at the slipline endpoint. 6-14 6 . 2 . 4 . 2 Redefinition of Target and Plug Masses When the plug edge is extended and the plug package boundary is enlarged, the code converts some of the target material into plug material. To do this, the code computes the volume to be associated with plug material in cells containing the plug edge extension. Using the intercepts of the plug edge extension with the cell boundaries, this volume is computed, as denoted by the shaded regions in the figure below. The initial and final extensions of the plug edge are special cases. The initial extension usually crosses one or more cells that contain projectile as well as target material. In this case, any curvature of the top of the projectile is ignored and the entire top surface of the projectile is assumed to have the same y-coordinate as the beginning point of the plug surface. The partial volume associated with the plug, therefore, is measured from the approximate position of the top of the projectile rather than from the bottom cell boundary, as illustrated below. 6-15 TAH Clii I'l.l C. I'lu: ■ MAT I The partial volume associated with the target is defined as the total cell volume minus the plug and projectile partial volumes, where the projectile partial volume is based on its mass and density. In the second special case, when the plug edge is ex- tended to the free surface, the new top boundary of the plug is the top of the target rather than a grid line. Ignoring the small curvature of the target, the partial volume associ- ated with the plug in the cell containing the endpoint is based on the y-coordinate of the new endpoint of the slipline, as illustrated below. TAR(jhT 6-16 Once the volumes associated with the plug material are computed, the target material in the region of the plug ex- tension is converted into plug material. The pure cells to the left of the plug surface endpoint are converted simply by changing the value of MFLAG for those cells from 3 to 2. (Package 3 is the target; package 2 is the plug.) The XMASS, SIE, and RHO arrays are redefined for the interface cells that contain both plug and target material. Given m = MFLAG(k) - 100, the plug mass is stored in XMASS (2,m) and is defined as the product of the volume of the plug material for that cell and the average material density for that cell. XMASS ( 2, m) = V plug av (6.4) where av "‘cell^'^cell (6.5) ;'he target material is stored in XMASS(3,m) and is the difference of the total cell mass and the plug mass: XMASS(3,m) = AMX(k) - XMASS(2,m) .(6.6) The density of both the plug and target materials is defined to be p , i.e., RHO(2,m) = RH0(3,m) = p . The specific internal energy of both materials is defined to be the average specific internal energy for the cell: SIE(2,m) » SIE(3,m) = AlX(k). The procedure for an interface cell that will contain all three material packages varies from that described above. In this case, the volume of the target material is used to compute its mass, so that m liili XMASS(3,m) - • p,„g„ C6.7) where V^ ^=Vn-V, -V . *^i^. The density of target cell plug projectile * the target material is used to define the plug mass: XMASS(2,m) - ( 6 . 8 ) This follows from the fact that the plug is simply a part of the target. Therefore the density and specific internal energy of the plug material is defined to be the same as the density and specific internal energy of the target material: RH0(2,m) = RH0(3,m) SIF.(2,m) = SIE(3,m) All target material is converted to plug material in the interface cells to the left of the plug edge and at the free surface, i.e., the plug quantities are defined: XMASS(2,m) = XMASS(3,m) RH0(2,m) = RH0(3,m) SIE(2,m) = SIE(3,m) and then the target quantities are set to zero: XMASS(3,m) = 0 RH0(3,m) « 0 SIE(3,m) » 0 6-18 CHAPTER VII -o I' j, DESCRIPTION OE INPUT VARIABLES AND RESTART PROCEDURES GENERAL CAPABILITIES AND LIMITATIONS A problem generator which uses the position of material interfaces to define cell quantities has been incorporated into the HELP code. Currently the generator includes a sub- routine, TSETUP, whicli allows the user to easily generate material interfaces which are composed of straight line segments or portions of circles or ellipses. However, the method is applicable to any topological configuration, provided the user generates the tracer particles wliich define the material interfaces. The following general capabilities and limitations of the HELP code should be noted before applying it to a specific problem: 1. The .number of material packages in a calcu- lation is limited only by the dimensions of the material arrays, which can be increased up to the limit of the user's computer core storage. (See Appendix A for dimensioning the arrays.) 2. A single cell can contain any number of material interfaces and any number of distinct materials . 3. Two material packages can actually contain the same material (use the same equation of state) but have different initial conditions (velocities, energies and/or densities). 4. A package can be divided into any number of spatially disconnected subpackages. (See instructions for defining material tracer particles in Section 7.2.5.) 7-1 5 . 10 , The material for each package must be homogeneous; i.e., each pure cell within a package must have the same density, velocity and internal energy. This restriction can be removed by appropriate user-supplied statements in FILGRD. The grid boundaries cannot serve as free surfaces because the code assumes, for example, tliat a package which extends to the right edge of the grid is infinite in the x-direction. The grid must have at least three rows and/or three columns. The DATA statements in COMDIM and EQST can be redefined if the user wants to model a material which is not in the EQST list of materials modeled by HELP. An inert material must be assigned a material number from 1 to 19, inclusively. A gas or high explosive must be assigned a material number from 20 to 30, inclusively. In order for the pressure iteration to converge, any equation of state which is used in place of the Tillotson, the y-law or the JWL, must be continuous and liave a negative slope every- where in the P-V plane, or at least prevent the iteration from entering any region which violates these restrictions. (See Section 4. 6. 1.3.) A calculation can be made in plane or in cylindrical coordinates. (See definition of IGM in Section 7.2.1.) The top and right grid boundaries are trans- mittive. The left boundary is reflective (an axis of symmetry in cylindrical coordinates) 7-2 The bottom boundary is the only one which can be either transmittive or reflective (see definition of CVIS in Section 7.2.1). 7.2 INSTRUCTIONS FOR GENERATING A PROBLEM Instructions for generating a liRLP problem are given in the sections that follow. An abbreviated set of input in- structions on keypunch forms is given in Appendix B. 7.2.1 Z-Block Variables Most of the options, flags, and problem parameters (such as the number of rows and columns in the grid) are variables located within the first 150 words of blank common and are defined by the first set of input cards. These variables are called Z-block variables because they all follow Z(l), the first word of blank common, and can be referenced using an appropriate subscript, between 1 and 150, for Z. For example, CVIS and Z(27) have the same location in blank common. Not all variables in the Z-block are defined by input cards; some are defined and updated by the code as the calculation proceeds. All Z-block variables are written on the restart file, and most are parameters which must be saved in order to restart a calculation. However, the variables such as UN93 or NUN32, which correspond to Z(93) and Z(32), respectively, are not currently being used by the code and are available to the user to define additional problem parameters as the need arises. The Z-block variables are defined by cards that are read by subroutine CARDS. The format (II, 15, II, 7E9.4) is used to read the variables lEND, LOC, NUMVfPC, (CARD(I), 1=1, NUMWPC) which are defined as below: lEND IEND=0 indicates that the variable (s) defined by this card are all typed real . I I inND*2 indicates that the variable (s) defined by this card are all typed integer. inND“l indicates that the variable(s) defined by this card are all typed real, and this is the last card CARDS will read until it is called again by INPUT. (Note: a setup deck has three cards where IEND=1, a restart deck has only two such cards. Wlien generating a problem, INPUT calls CARDS three times; when restarting a problem, INPUT calls CARDS twice.) LOC LOG is the location in blank common of the first variable being defined by this card. NUMWPC NUMWPC is the number of consecutively stored Z-block variables of the same type being defined by this card. CARD(I) The value assigned to the I^^’ variable. NOTE: even if the variable is typed integer, the value must be entered witli a decimal point. (See sample input decks in Chapter XI.) The Z-block variables defined by input cards are listed and defined on the following pages. The type (integer or real) and location in blank common are indicated for each variable. The code assigns non-zero default values to some of these variables at the beginning of subroutine INPUT. The non- zero default values are given in parentlieses after the variable definition. If no value is indicated, the default value is assumed to be zero. INPUT Z- BLOCK VARIABLES Column Location Variable One in Blank Name Flag Common Definition (Default Values are in Parentheses) PK(1) PROB NFRELP NDUMP7 ICSTOP NUMREZ KUNITR Problem number. Any number from .0001 to 99.9999. It must be identical to PROB. Problem number. (Same value as PK(1).) Must be included in the setup deck. Gives frequency of "long" EDIT prints, which give velocities, energy, com- pression (or density) and stresses for all cells in active grid. A "short" EDIT gives this information only for the first column of cells, e.g., when NFRELP®2, every second EDIT print is "long," when NFRELP^S, every fifth EDIT print is "long", etc. Gives frequency of restart tape dumps relative to EDIT prints. (Program dumps only when EDIT prints.) E.g., if NDUMP7=5, a restart dump will be made every 5th time EDIT prints. Gives cycle at which execution will stop if the user wants to specify a cycle, rather than a value of T (time), on which to stop. If stopping on time omit this card. The total number of times the grid will be automatically rezoned. If not automatically rezoning the grid, omit this card. The name of the file INPUT will read from. (7) ■W ■ INPUT Z-BL( CK VARIABLES Variable Name Column One Flag Location in Blank Common Definition (Default Values are in Parentheses) IPR 2 15 Maximum number of iterations CDT will perform when attempting to equilibrate the pressures of materials in a multi- material cell. If the pressures are not within an epsilon (PRCNT) of the average pressure after IPR iterations, an error exit occurs. (35) PRCNT 16 Convergence criterion tor the pressure equilibration of material in a multi- material cell. All material pressures P^, must satisfy the following: 1 (P - P.)/P 1<PRCNT. (.001) 1 '■ av 1 av ' KUNITW 2 17 The name of the file EDIT and SETUP will write on. (7) I CM 2 21 When IGM = 0 code uses cylindrical coordinates. When IGM = 1, code uses plane coordinates. DM IN 24 The maximum relative error in the energy sum that the user wants to tolerate. The error is computed as follows : ,KMAX \ - ETH^/ETH where E. is total energy of cell k and ETH^is theoretical energy in the grid- -computed in SETUP and updated in HPHASE, TPHASE, SPHASE, ADDENG, REZONE, and RNDOFF. (lO'^) CVIS 27 A flag that describes the condition of the bottom grid boundary. If CVIS = 0, the bottom boundary will be reflective. If CVIS = - 1, the bottom boundary will be transmitt ive . 7-6 INPUT Z-BLOCK VARIABLES Column Location Variable One in Blank Name Flag Common I MAX JMAX MAPS NUMSCA PRLIM Definition (Default Values are in Parentheses) The number of columns in the grid. If planning to rezone the grid in the x-direction, IMAX must be an even number. IMAX also must be at least 3. The number of rows in the grid. If planning to rezone in the y-direction, JMAX must be an even number. JMAX also must be at least 3. When MAPS > 0, part of the EDIT print is a set of symbolic maps of the compression (or density), pressure, radial and axial velocities, and specific internal energy of the cells in the active grid. MAPS=1 gives a compression map. MAPS=2 gives a density map. The number of times the cycle or time interval between EDIT prints is in- creased. (See PRLIM) The time or cycle at which the EDIT print frequency is ciianged. -Example : you wish to print every 10 sec until T = 10"' sec; and every 10"' sec until T = 10'^ sec, and every lO"^ sec there- after. Set: - PRDELT =10'® IPCYCL = 0. PRLIM =10"^ PRFACT =10. NUMSCA = 2. 7-7 input Z'BLOCK variables Column Location Variable One in Blank '"'Same Flag Common PRDELT PRFACT II 45 46 47 12 IPCYCL TSTOP 48 49 50 IPLGRT 56 CPe Definition fault Values are in Parentheses) rncrcfsfa.b^ ?.;e same factor interval. omit’*NUMSCA, PRLIM, PRIACT. The time Csec) '’®^^p®^j^tervals , omit If printing on cycle interva , this card. The factor "j!j;?a‘j%hrprint‘limit PlM) rra’lSL^d. CSee PRUM) The number of have "°"-fi”J^l2''define the "active" !S-;ifrei?r^^^ire%re%rtfvrrria., See-rn-l/ro re^ocfti'erof itr plus 2. U nf rvcles between EDIT ™!„"ts U printing on time intervals. omit this card. 1 nf T ftime) at which the ^ The value of When stopping calculation will Pj.j^g,^,QP ^ 0 ),^set note: This^cajd has^a ,1_ SudeHS bfth the setup de J and the restart deck. The. rightmost column of the plugging To^i^e ;!giror;fe^ XP ctea -ara‘r’'oSit'^hifcirrwre^^rot 4* V> A II A. V — — rvliiooinp option. 7-8 :i 4 1 r 'vwfvaKhii INPUT Z-BLOCK VARIABLES Column Location Variable One in Blank Name Flag Common Definition (Default Values are in Parentheses) PLWMIN IPLGBT IPLGTP GAMMA NMAT CYCMX CYCPH3 NTRACR The required minimum specific plastic work J (ergs/gm) whicli the material near the i top of the vertical edge of the plug | must have before the plug edge is ex- I tended to the next row. See Section a 6,2.1. Omit this card when not using | the plugging option. | The bottom-most row of the plugging region of the target. (Usually the \ first row of the target.) Omit this card when not using the plugging option. The top-most row of the plugging region j of the target. (Usually the top row of the target.) Omit this card when not using the plugging option. . Y in P=(y-1)pE for material #20, an ■'< ideal gas, ^ i The number of material packages, ex- cluding the void package. The number of passes through subroutine INFACE each cycle to minimize transport noise near interfaces. (2) The number of passes through subroutine SPHASE (strength phase) each cycle. When the calculation is purely hydrodynamij (no strength effects), set CYCPH3*-!. (1) The number of material tracers per cell diagonal. Used by ADDTCR when adding tracers. ADDTCR adds tracers only when ’ the distance between tWo consecutive j tracers in a specified region is .1 greater than a cell diagonal divided 'M by NTRACR. MINX, MAXX, MINY, MAXY J specify the region of the grid considered^ by ADDTCR. NADD specifies the ^ frequency with which ADDTCR is called. 7-9 Variable Name NMXCLS NTPMX NTCC SIEMIN EMIN PMIN INPUT Z-BLOCK VARIABLES Column Location One in Blank Flag Common Definition (Default Values are in Parentheses) 73 78 81 82 85 86 The maximum number of interface cells to exist in the grid on any one cycle during the calculation. This maximum should correspond to the dimensions of the material arrays (XMASS, SIE, etc.). See Appendix A. The maximum number of tracer particles per material package. This maximum should correspond to the dimensions of the TX and TY arrays. When NTCC 0, NTCC passive tracer particles are initially placed in the center of every fourth nonempty cell and subsequently moved through the grid with the material. These are used only for producing particle plots and do not affect the material behavior. When this card is omitted (NTCC = 0) , cell centered tracers are not computed unless the plugging option is being used, in which case they are automati* cally generated in the plugging region of the target. (See Section 6.1.3) A cutoff on the total specific internal energy increment of a cell in TPHASE. This cutoff prevents small numerical signals from enlarging the active grid. ( 105 ) The minimum value of specific internal energy to be used in the ideal gas equation of state. (10^) A pressure cutoff. If lP(k) | < PMIN, then P(k) =0. (5 x 10°) 7-10 ■ " ■ v.-.^ .-.. ,.’«r INPUT Z- BLOCK VARIABLES I Column Location Variable One in Blank Name Flag Common INTER 2 87 NODIIMP NLINER 2 105 NVRTEX 2 109 Definition (Default Values are in Parentheses) A special editing flag for debugging: »2 for debug prints in CUT “3 for debug prints by EDIT after HPHASE and SPHASE *5 for debug prints ia TPHASE *7 for debug prints in SPHASE =11 for debug prints in UVMOD ■13 for debug prints in UVCALC (See Section 9.2 if more than one debug is desired.) The flag which initiates the rezoning of the grid. If NUMREZ > 0, this flag will be set automatically by the code. (See Section 8.4) However, the user can force a rezone on any restart cycle by setting REZ ■ 1 in the restart deck. (lEXTX and/or JEXTY must also be defined in order for the grid to be rezoned.) When NODUMP = 1, EDIT will not write any restart dumps. (This flag over- rides the NDUMP7 option, but does not prevent SETUP from writing the cycle 0 dump . ) The package number of the shaped charge liner material that forms the jet. Used only when calculating the collapse of a liner. The second index of the void tracer particle that is at the vertex of the void closing region, i.e., (TX(NVOID, NVRTEX), TY(NV0ID, NVRTEX)) are the coordinates of the vertex point. Used only when the code's automatic void closing routine is to be activated. (See Section 8.4) ' 4.V 7-11 INPUT Z-BLOCK VARIABLES Column Location Variable One in Blank Name Flag Common Definition (Default Values are in Parentheses) ROEPS 110 A round-off epsilon used to defipe cutoffs in the calculation. (10'“’) PLGOPT 111 If the user is generating a plugging calculation, PLGOPT must be set equal to 1. The special plugging mechanisms are activated only if PLGOPT = 1. (See Chapter VI before using the plugging option. ) NSLD 2 112 The maximum number of cells the user expects will contain the slipline on any one cycle. This maximum should also depend on the dimensions of the slipline arrays (MSLD,. etc. ) . FINAL 113 The final value of the stability fraction used in determining the time step. See STAB. (.4) NOSLIP 2 115 If no sliplines are to be generated, set NOSLIP = 1, which causes t.Le in- structions and routines that affect only slipline cells to be skipped. LVISC 2 116 A linear artificial viscosity term is added to cell boundary pressures if LVISC = 1. (See Section 2.2.2.. 5.) NADD 2 118 A flag for automatically adding material tracer particles in the region specified by MINX, MAXX, MIMY, MAXY. If NADD = 10 ADD'iCR is called every cycle which is a multiple of 10, (e.g., cycle 20, 30, ... 250, 260, etc.). NOTE: NTRACR (=Z(72)) must also be defined before material tracers can be added by ADDTCR. MINX MAXX 2 2 119 ^ 120 These parameters specify the left and right columns and the bottom and top , rows, respectively, of the region in MINY 2 121 which ADDTCR will add tracer particles. These variables must be defined when MAXY 2 122 , NADD f 0. 7-12 INPUT Z- BLOCK VARIABLES Column Variable One Name Flag Location in Blank Common Definition (Default Values are in Parentheses) lEXTX JEXTY STAB DTMIN CRATIO A rezone flag. The grid will be rezoned in the x-direction only when lEXTX « 1. Omit this card if the grid is not to be rezoned, or if it is to be rezoned in the y-direction only. A rezone flag. The grid will be rezoned in the y-direction only when JEXTY * 1. Omit this card if the grid is not to be rezoned, or if it is to be rezoned in the x-direction only. The initial value of the stability fraction used in determining the time step. If FINAL * 0, STAB is constant. Otherwise, its value is doubled each cycle until it reaches the value of FINAL. (10’3) The minimum value for the time step. The program gives an error exit if CDT calculates a At < DTMIN. (10"^^) The compression of the material in adjacent cells is computed. If the ratio of these cells’ compressions is greater than CRATIO, their pre'.sures are inverse compression weighted and their velocities are compression weighted to define cell boundary values in HPMASE. (See Section 2 . 2 . 2 . 4 . ) (10^) BBAR EMOB A constant used in the calculation of the local sound speed of materials other than ideal gases and high explo- sives. C = C^ BBAR • >/|P| . (.5) Dummy end card of a setup deck. This card is read by CARD^ after all other input cards have been read. “TlTis card must not be included in a restart deck. Set EMoB = 0 in the setup deck. 7-13 7.2.2 Cell Dimensions The second set of inpjit cards define the cell dimensions. The X and y dimensions of the cells can vary; however, the dimensions of contiguous cells should not vary by more than 10% except in parts of tlie grid well beyond the region of interest. Large aspect ratios (e.g., AX/AY > 2) are also discouraged in significant regions of the grid. The cards which define the DX, DY arrays are road by subroutine SETUP. A set of cards defining DX are read, then a set defining DY. Each card specifies up to four different DX or DY values. The informatioa is read into the NNT and TEMP arrays: (NNT(L),L = 1,4), (TEMP(L),L = 1,4), where NNT(L) is the number of columns (or rows) that have a DX (or DY) dimen- sion equal to TEMP(L). After the x-dimension (DX) of all columns has been defined, the next NNT is set equal to 999, which indicates that the entire DX array is defined. Then the DY array is similarly defined by the next group of cards. After the y-dimension (DY) of all rows has been defined, the next NNT is set equal to 999, indicating that all y-dimensions have been defined. An error exit will result if the user does not define exactly IMAX DX's and exactly JMAX DY's. 7.2.3 Initial Density, Velocity, and Specific Internal Hnergy of Each Material Package In addition to the material tracer particle positions, the definition of material packages requires identification of MFLAG(2), the flag of the cell in the bottom left corner of the grid, and the material code number and initial conditions of each material package. The card which follows those defining the cell dimensions defines MFLAG(2). The user should be able to predict from the placement of the interfaces in the grid whether cell k = 2 will be a pure cell, an interface cell, or a void cell at 7-14 t = 0. If it will be a pure cell of package n, then set MFLAG(2) =» n. If it will be an interface cell, then set MFLAG(2) ■ -1. (In the latter case its value will be changed to be greater than 100 after the interfaces are processed by CALFRC and VOLFND) . If it will be a void cell, then set MFLAG(2) » 0. After MFLAG(2) is defined, a set of cards which defines the material code number and initial conditions of each material package is read. The first card in this set defines material package 1, the second card defines material package 2, etc. NMAT cards must therefore be read to define all the material packages. The material code numbers for 24 materials are listed in EQST as well as in Tables 2.1 and 2.2 of this report. The material code numbers are stored in the MAT array and are used whenever the equation of state constants (Pq»A,B,y» etc.) are referenced. These equation of state constants are defined in DATA statements in EQST and in COMDIM, the "included" element. For example, if package 3 is lead (material code number 10) then MAT(3) = 10, ESCAPA(IO) is its bulk modulus. A, and RHOZ(IO) is its normal density, p^. The density (g/cc), specific internal energy (ergs/g) and the radial and axial velocities (cm/sec) of all the material in each package are specified by these cards. The knowledgeable user can vary these quantities within the package by adding to subroutine FILGRD coding which would redefine the appropriate variables. When defining the initial density, RHOIN, of a material, the user should be aware of the normal density, RHOZ, specified by the code fcr that material. (See Table 2.1 and 2.2.) The user can change the definition of for a given material by redefining the appropriate RHOZ variable in the 7-15 DATA statement in COMDIM, the "included" element. The user can also specify an initial density, RHOIN, which is dif- ferent from a material's normal density, RHOZ. See Appendix B for the exact format of the cards which define MI-’LAG(2) and tlie material code numbers and initial conditions of the material packages. 7.2.4 Strength Constants for Each Material Package The shear yield strength and the tensile strength of each material package are defined by the next set of cards. The shear yield strength constants, Y , Y, , Y^, E , G, are o o l ^ m all set to zero if the material has no yield strength. If, however, the material has no tensile strength, AMDM(n) is set to 1. The definition of the yield strength terms is given in Section 2. 3. 2. 3. The tensile failure threshold is described in Section 2.3.3. The setup deck must contain NMAT cards defining strength properties, even when none of the packages have strength and SPHASE is being hypasse<!. See Appendix B for the exact format of the cards which define the strength constants. 7.2.5 Material Tracer Particles A material package boundary is generated by dividing it into a series of straight lines, and/or arcs of circles and/or arcs of ellipses. A set of 2 or 3 cards is read by subroutine TSETUP to define each one of these segments. The first of these cards specifies: (1) LTYPE, what kind of a segment to generate - a horizontal line, a vertical line, a diagonal line, an arc of a circle, or an arc of an ellipse; (2) MPN, which material package is being generated; and (3) NPTS, how many tracer particles to put along the segment. If the segment is a 7-16 straigljt line, only one other card, which defines the end- points of the line in centimeters, is read. However, if the segment is an arc, two other cards are read which define the beginning and ending angles of the arc in radians and the location of its center, or foci, in centimeters. iVhen defining the last segment of a material package or * subpackage boundary, the user must assign a negative value to LTYPE. This signals TSETUP to generate a dummy tracer particle which marks the end of a tracer string for a sub- package or package. The coordinates of this dummy particle are (-1000, 0), and an error exit will result if this dummy particle is not the last particle in the set of tracers for a material package. After all segments of all packages, in- cluding the void package, have been generated, a flag card, setting LTYPE * 100, is read by TSP-TUP. This signals that all the material tracer particles have been defined. As the segments are generated they must form a continuous boundary interrupted only by the grid boundaries. (It is unnecessary to generate tracers along the grid boundaries.) Furthermore, the segments need to be generated in an order such that the interior of the package lies on the left as one moves between any consecutive pair of tracers. This latter requirement also dictates which endpoint must be specified first on the input cards. Since every segment of an interface is described by two identical strings of tracer particles (one for each material package lying on opposite sides of the interface) , A material package can be divided into any number of spatially disconnected subpackages. The tracer string for each subpackage must end with a dummy tracer. 7-17 tmi the number of tracers generated along a segment must be identical for both packages, the order of the two sets of tracers along that segment being exactly opposite. A good rule to use in deciding how many tracers to place along a segment is to generate two per cell edge. In regions where great distortions are likely to occur, the tracers may be generated more densely, although activating ADDTCR (see definition of NADD in Z-block variables) helps to insure that the tracers will not become too sparse. The user cannot generate more than NTPMX (a Z-block input variable) tracers for a given material package boundary; in that event TSETUP will print a message and stop. The definition of NTPMX should match the dimensions of the tracer particle arrays, TX, TY. See Appendix B for the exact format of the cards which define the material tracer particles. 7.2.6 Slipline Endpoints When sliplines are being generated, SETUP reads a set of cards which defines the siipline endpoints and specifies each material package as a master, a slave, or neither, (If there are no sliplines in the calculation and if NOSLIP (a Z-block input variable) is set equal to 1, the slipline cards will not be read by SETUP and must be omitted from the input deck.) One card is read for each package; the first card defines the first package, the second card defines the second package, etc. Each of the first two variables, MASTRD, NSLAVD, is set either to the package number or to zero, depending on whetlier the package is a master or slave. Each of the next two variables NBGMD, NBGSD is set either to the index of the first tracer of the slipline or to zero, 7-18 t again depending on whether the ^^NDSD, is ^:rri::;r^ ti: rder:i the iast tracer ei the Siipune :: to tore, depending on whether the package ts a .aster slave . . ::: rr::.:: irster'^Tt'o::' nrt^rtfe":hi:h";:c'rages are the .asters Id which are the slaves.^ cellTal centrinriire :ran^:::”is::;a::::the-^^^^^^^^^^^^^ rac'rgL'u T^Z'Z7Z7t^.. sa.e velocity co.ponents. '"’“see Appendix B for the exact for.at of the cards which define the slipline endpoints. 7 2 7 ..Cinition of Hig B Pvpiosive hetonationloints SPTUP reads two sets of cards in order to define the P7ia.ary and secondary initiation points of one or .ore g^ :-rorrr;J-::rr,, 7-19 starting with the point having the minimum detonation time. There is one card for each detonation point, plus an end flag card. The second set of cards defines the approximate region of detonation associated with each detonation point. These regions can overlap. There is one card in this set for eacli detonation point, hut no end flag card. If the calculation does not involve a high explosive, one blank card must be inserted in the deck. See Appendix B for the exact format of the cards which define the detonation points. 7.3 RESTART PROCEDURES During the calculation, HELP periodically writes on a tape or drum file tlie problem parameters and the current state of the material in each cell. By reading this file and a few input cards, the code can "restart" and continue a calculation from an intermediate point. In general, a restart deck consists of five parts: 1. the heading card 2. the "restart" card 3. cards redefining various Z-block variables (if desired) 4. the "end of data" card 5. cards defining slipline endpoints (if the slipline option is activated; i.e., if NOSLIP = 0) . The heading card can contain any alphanumeric symbols between columns 2 and 72. The "restart" card defines three variables in the PK array whicli are described below. (Note: since the "restart" card contains three words, it must have a "3" in 7-20 ■ 1 column 7. See Section 7.2.1.) Also, since it is the only card processed by CARDS before the restart file is read, it must have a "1" in column 1. The third part of the restart deck allows the user to redefine any of the Z-block variables he desires; e.g. variables which make the code add tracer particles, change edit frequency, obtain debug prints, etc. If the user is stopping the calculation on cycles, he must define a new stop cycle, ICSTOP, and this is the only case in which there must be a card in this third part of the restart deck. Tlie fourth part of the restart deck is the "end of data" card. This card, which must have a "1" in column one, redefines TSTOP, the value of time at which the current run will stop. If the calculation is stopping on cycles, rather than time, a zero value is entered. The last part of the restart deck, which defines the slipline endpoints, is necessary only if the slipline fcapability is being used in the calculation. (See Section 7.2.6 for discussion of defining the slipline endpoints.) As noted in Section 8.3, the slipline can be completely redefined or even removed during a restart. The user should be particularly aware of the possible change of endpoint in- dices if ADDTCR (see Section 8.2) and/or VDCLOS (see Section 8.4) have been called during the previous run. In general, unless the user wishes to redefine the endpoints of the slipline, he sl\ould input those indices listed in the EDIT print of the cycle from which the calculation is being restarted . The restart deck of 6 cards listed below could be used to restart at cycle 78 a calculatipn whose problem number is 26.5 and to run it until T = 3.0 ysec. A "short" EDIT print of the restart cycle is indicated by PK(3) * — 2. The third card redefines NADD, which controls the frequency 7-21 1 of the calls to subroutine ADDTCR (let's assume that MINX, MAXX, MINY, MAXY were defined previously). Until NADI) is redefined again ADDTCR will be called every cycle which is a multiple of 10. The slipline cards indicate package 1 is a slave and package 2 is a master; 23 tracers in each package define the slipline. j 8-16 1 17-25 I 26-34 | 35-43 | 44-52 j 53-61 j 62-70 | I 2 3‘)5a7»’o I2 3*)567^<»{j( . | 2 H 1 23jM54789o J ^3‘)S67 B vc i|23 MSA 789o tBl|324.! I lali to. C 1 2 n' IHf>n - aLUiiIuUn P‘PACT 79. «2. -CA n 147 Let's assume that this same problem was to be restarted on cycle 325, the slipline was to be removed and the calcula- tion was to be stopped on cycle 326. The following 5 input cards would be used; ! 8-16 I 17-25 1 26-34 j 35-43 I 44-52 1 53-61 I 62-70 1 1 2 3 H b 4 7 8 9 7 I 2 3 9^4 7 8 9j I 2 3 J 54 7 0 9 ^ | 2 3 9 b 47 8 9 p I 2 3 9 b>7 7 8 9 q 1 151324,5 2 I ) bl I • 2 7132a. I Solo. iKnn - /LUnlHUi Impact 325. -2. Card 3 sets NOSLIP = 1, which deactivates the slipline option. Thus it is not necessary to define the slipline endpoints in this restart deck. n 7-22 rt should be noted that any blank common variable that is not in the Z-block can be changed using the CARDS routine format if the user can compute its location in blank common. However, usually it is easier to code the changes into subroutine INPUT (following statement number 40) after the restart tape has been read and the pressure array has been initialized. Of course, this is the only way variables not in blank common can be redefined by the user. In the course of some calculations, problems can arise because of the position of the material tracer particles. % This is most likely to occur in calculations which involve severe velocity gradients or which employ VDCLOS, the auto- matic void closing routine (see Section 8.4). In the first case, the speed at which the tracer particles become sparse makes it easier for the tracers to cross each other, destroying the logic of FRACS. In the second case, it may happen, for example, that a tracer particle may pass to the left of the vertex point. This could create a situation in which VDCLOS forces the tracer string to cross itself when the tracers are moved to the new vertex position. If these situations arise, it is necessary for the user to add, delete, or reposition some of the interface tracer particles on a restart cycle (probably the last one whic*^ was written before the problem arose) in order for the calculation to proceed smoothly. The procedure is not difficult, but. must be carefully thought out before being implemented. Basically these are four principles which must be firmly adhered to: 1. Every interface tracer of each package, including the void, has a corresponding tracer from the other t package(s) which adjoins that interface. In the case of a corner point among more than two packages. 7-23 such as the intersection of the liner, casing, and HE package in a shaped charge calculation, there may be more than one other corresponding tracer. If any tracers from one package are added, deleted, or moved to a new location, its corresponding twin(s) must also be added, deleted or moved to the identical new location. 2. If tracers are being added to or deleted from a material or void package n, the value of NMP(n), which stores the number of particles, including the dummy endpoint, defining the package n boundary, must be increased or decreased accordingly. 3. When adding or deleting tracers, the user must in- sure that the dummy endpoint (-1000, 0) is the value of the final tracer of each package and sub- package . 4. When adding or deleting tracers from calculations which use the automatic void closing routine, and/or the slipline option of the code, care must be taken to insure that the values of NVRTEX and/or the NBGMD, NENUMD, NBGSD , NHNDSD arrays reflect the changes being made. An example of a calculation which uses the automatic void closing routine will serve to illustrate the above procedures. Consider the situation in Figure 7.1, in which the package 2 tracer particles near the vertex of the void closing region have become too sparse to satisfy the conditions considered by VDCLOS, i.e., even if the pressure is positive in the cell containing the vertex point, the angle between the two surfaces is greater than .2 radians and is becoming larger. (See Section 8.4.2.) Therefore, this void must be closed by t)ie user. V-24 In this example, the current vertex of the void closing region is defined by the 79^ package 1 particle, the 17^ package 2 particle and the 140^^ free surface particle (package 3). Let us assume the user decides the new vertex i" K point will be defined by the 18th package 2 tracer. The 75^ package 1 tracer will be moved to the new vertex point, while package 1 tracer particles 76 through 78 and free surface » tracers 140 through 144 will be removed. Figure 7.2 illus- trates the modified interface positions resulting from these changes to the tracer particles. Note that the 139^^ free surface tracer is now at the vertex, so NVRTEX must be set equal to 139 for VDCLOS to operate correctly. These changes to the tracer particles can be made in subroutine INPUT as indicated by the following lines of FORTRAN: i I I I I I c 30 MO c c r C C 351 C 352 C ••• initialise T'-STORAgE # DO MC P ( M «C *C • •• rHA luL COHtiOij I F ( NC • nT . I ‘K, » (,'J TO 355 TX( I ,7S)«TX(2,UI Ty(i,75)»TY«2,lB> LO 3bl L»7V,96 T X < I ,L-3 I »T X ( I tL I Ty(l,L-3)«lYU»Ll CONTI NUf. NMP I I ) •'»3 DO 352 L-IlSiZlV TXC3,L-5)«Tx(3tu> TY(3,l.-5»«TY(3»U continue NMP ( 3 ) *2 I M yapiaRLEs here on a lines added to INPUT to y change tracer particle positions c 355 NVR ’EX»139 CONTINUE V ke5Tap r I.. - 7-26 I (grid lines) Figure 7. 2 --Tracer particle positions after being redefined by user. In this cast, the changes arc being made after cycle In this , restart 145 has bcfin read off the restar dump will have the unmodified positions (Figure 7.1) of t ^ 0 ^ rticles; however, in the next cycle to be exocu ed, C tnfaCE will process the modified tracers as they number x4b, INFACE win pro^c appear iu Figure 7 . 2 , and the next restart dump wrll the modified tracer positions. CIIAPTLR VIII INTERACTIVE CAPABILITIES There are several ways the user can interact with the calculation as it progresses. The area encompassed by the grid can be enlarged and the resolution reduced by rezoning; the density of the tracer particles in a specified region can be increased; interfaces can become slip surfaces or slip surfaces can be removed; and a vertex of two intersecting free surfaces can be specified where the code will automat- ically "close a void." These interactive capabilities will be discussed in more detail in the sections that follow. 8.1 RE ZONING THE GRID The HELP rezone reduces the resolution of the grid by combining cells and enlarging the area encompassed by the grid, keeping the number of cells in the grid constant. The HELP rezone combines two adjacent cells into one. Therefore, the grid can be rezoned radially and/or axially only if there is an even number of cells in the direction to be rezoned. The user should note that this rezone routine is not capable of rezoning plugging problems. 8.1.1 Combining Cells When two cells, k and kn, are combined into one, the masses of the two cells are summed, and the velocity components (u',v*) of the new cell are a mass weighted average of the velocity components of the original cells, thereby conserving mass and momentum. 5 m' = m(k) + m(kn) u' = (u(k)m(k) + u(kn)m(kn))/(m(k) + m(kn)) v' = (v(k)m(k) + v(kn)m(kn))/(m(k) + m(kn)) The material masses and velocities of interface cells are averaged in the same manner. For example, given m=MFLAG (k) - 100 and mn*MFLAG(kn) -100, the mass of material n in the new cell, will be m^'^ = XMASS(n,m) + XMASS(n,mn) The internal energy of the new cell is defined so as to conserve the total energy of the two cells being combined. Therefore, the internal energy of the new cell, IE*, becomes the sum of the internal energies of the two cells plus the thermalized kinetic energy resulting from the definition of the new cell velocity components. IE* = [lE(k)*mrkj + IE(kn)*m(kn) + m(k) • [(u(k)-u')^ + (v(k)-v’)^] + m(kn) • [ (u(kn) -u* ) + (v(kn) - v* ) ^]] / (m (k) + m(kn)) The internal energies of the materials in interface cells are redefined in the same manner. The deviator stresses of the new cell are a mass weighted average of the deviator stresses of the two cells, whereas the slipline angle and the detonation time associated with the new cell are simple averages of those quantities in the two original cells. 8.1.2 Adding Material The code automatically adds material above and/or to the right of the original grid, depending upon the direction 8-2 in which the grid is being rezoned. If, for example, the grid is rezoned in the x-direction, IMAX/2 cells are added to each row, since the combination of pairs of original cells has resulted in the original area being encompassed by half as many cells. The cells that are added on the right have the same x-dimension as the cell in the last column of the original grid. These added cells are also assumed to have the same initial properties (density, velocity and internal energy) that the material had in the initial problem setup. If the last coll in a row is an interface cell, the code assumes that the interface extends throughout the row being added; therefore, all the added cells ift that row will be interface cells. The user, in this case, should have generated tracer particles on a horizontal line beyond the area of the original grid. (Note: The rezone assumes that interfaces extending beyond the original grid are horizontal in the x-direction and vertical in the y-direction. Also, the rezone will not automatically generate a free- surface and a void unless they are an extension of a free surface and a void in the original grid.) Furthermore, the rezone routine, as it exists, does not add material on the left (in plane coordinates) nor below the grid. However, the knowledgeable user can remove these limitations quite easily. 8.1.3 Activating the Rezon e Basically the code implements two methods for rezoning the grid. One of these is an automatic rezone, in which the code automatically determines, during a run, when to rezone, does so, and continues the run to the desired final time or 8-3 cycle. The other method is a rezone which is triggered by the user through the restart deck. No matter which method is being used, the user must set the Z-block variables lEXTX and/or JEXTY to 1 to indicate the direct ion (s) , radial and/or axial, respectively, in which the grid is to be rezoned. As noted before, the grid dimension(s) (IMAX and/or JMAX) in the direction(s) to be rezoned must be an even number. 8. 1.3.1 Automatic Rezone In order to activate the automatic rezone capability, the user must define the Z-block variable NUMREZ to be the maximum number of automatic rezones he wishes to allow. Subroutine TPHASE then tests to see if any mass is trans- ported through the right grid boundary (if lEXTX = 1) and/or the top grid boundary (if JEXTY =1). If any mass is lost through the boundary being tested, the rezone flag REZ is set to 1 and the code rezones on the next cycle. The code automatically edits both before and after the rezone. Also the value of NUMREZ is decreased by 1 follow- ing each rezone so that a maximum of NUMREZ automatic rezones is allowed. 8. 1.3. 2 User-activated Rezone The user has the option of triggering the rezone through the restart deck by defining the rezone flag REZ * 1. In this case the pickup cycle is rezoned in the directions indicated by the lEXTX and JEXTY flags. The code gives an edit following the rezone and continues the execution of the problem. 8-4 8.2 AUTOMATIC ADDITION OF MATERIAL TRACHR PARTICLES By defining the appropriate variables (NADI), NTRACR, MAXX, MINX, MAXY, MINY) , the user can insure that the tracer particles will not become too sparse in the regions of the grid where there are great expansions or distortions of the material interfaces. Every cycle which is a multiple of NADD, subroutine ADDTCR is called. ADDTCR searches the region bounded by the MINX, MAXX codumns and the MINY, MAXY rows. If a pair of consecutive tra^rs in that region are more than 1/NTRACR of a cell diagonal apart, ADDTCR adds enough tracers between the pair to achieve the desired density of NTRACR tracers per cell diagonal. (Section 7.4 describes how tracers can be "manually" added by inserting coding changes in INPUT when a calculation is restarted.) If, by adding a set of tracers, tlie dimension (NTPMX) of the tracer particle arrays would be exceeded, ADDTCR instead prints a warning and sets NADD to zero. The code will not call ADDTCR until NADD is redefined in the restart input deck, which, in this case, should be done only after the tracer particle arrays have been enlarged and the code recompiled. Also, ADDTCR redefines the slipline endpoint indices and the void closing vertex, NVRTEX, when necessary. 8 . 3 REDEFINING SLIPLINES The sliplines in HELP can be lengthened, shortened, created or eliminated each time a calculation is restarted, simply by redefining the endpoints of the slipline. During the first cycle following a restart the code will reflect the change in the slipline definitions when it computes the value of THETA(m) for each interface cell. An interface cell whose interface is changed from a slip to a non-slip surface, will have all material velocity component set equal to the cell centered velocity components and the resulting thermalized kinetic energy will be added to the material internal energies. This will be done automatically in TPIIASE. The material velocity components of interface cells that have changed from non-slip to slip cells, will gradually become different from the cell -centered values as the master and slave materials are treated differently in TPHASE and HPHASE. It should be remembered (when redefining the endpoints of the slipline), that a slip cell has no shear yield strength. 8.4 CLOSING A VOID The continuous velocity field used to move tracer parti- cles in HELP prevents two particles on a collision course from ever actually meeting. As the particles come closer and closer to one another their velocity fields become more and more alike until they are being moved with almost identical velocity components. Therefore, if two free surfaces are moving toward one another, they will never quite meet, and the void will never quite close. The void closing routine in HELP will automatically close the void in one specified region of the grid where two free surfaces initially have a common point; this common point is called the vertex in the discussion that follows. 8-6 8.4.1 Identifying Void Closing Region The vertex point is identified by an input variable, NVRTEX, where (TX(NVOID, NVRTDX) , TY(NV0ID, NVRTHX) ) are the coordinates of the void tracer which is at the common point. As the void is closed, NVRTEX is redefined, identifying another free surface particle as the vertex of the void closing region. 8.4.2 Criteria for Closing a Void The void is closed if at least one of the following three criteria is met: 1. The pressure, PVRTEX, of the cell which contained the vertex particle at the start of the cycle is positive, and the angle between the two surfaces is less than .2 radians. 2. The two free surface interfaces have crossed. (This occurs only if the tracer particles have become too sparse in that region, or if the plugging option, which invokes special noncontinuous rules for moving tracer particles, is being used.) 3. The particle next to the vertex on one free surface is less than 1/10 of a cell diagonal from the other free surface. 8.4.3 Method for Closing the Void VDCLOS defines four points which are uSed in testing the criteria listed above as well as in the void closing procedure. Referring to Figure 8.1, (TXl, TYl), (TX2, TY2) , (TX3, TY3) are the void tracers indexed by NVRTEX- 1, NVRTEX and NVRTEX+1, respectively (i.e., (TX2, TY2) = [TX(NV0ID, NVRTEX), TY(NV0ID, NVRTEX)]) . The fourth point (TXP, TYP) is defined by VDCLOS 8-7 I such that the line through (TX3, TY3) and (TXP, TYP) is perpendicular to the interface line through (TXl, TYl) and (TX2, TY2). If one of the three criteria listed above is met, the following five steps are taken to close the void. These steps are described with reference to the example illustrated by Figures 8.1 and 8.2. 1. The package B tracer at (TX3, TY3) is moved to (TXP, TYP). 2. The package A tracer at (TXl, TYl) is moved to (TXP, TYP) . 3. The void tracer at (TXl, TYl) is moved to (TXP, TYP) . 4. The void tracers at (TX2, TY2) and (TX3, TY3) are removed from the string of void tracers, and NMP(NVOID) is decremented by 2. 5. Since the new vertex point is at (TXP, TYP), which is the location of the NVRTEX - 1 tracer (see step 3), NVRTEX is decremented by 1. Each time the void is closed, the position of the vertex point is changed, and the void closing region is redefined. VDCLOS then applies the same three criteria to the new void closing region, closing the void and redefining another void closing region if any of the criteria are met. This process continues until a void closing region is defined * Since PVRTEX is the pressure of the cell containing the vertex point at the beginning of the cycle, the first criterion is not applicable if the new vertex point is not in that cell. >1 ft i' Figure 8.2-- Void closing region after void is closed. The old positions are ^ represented by dashed lines, and the original position of (TXl, TYl) is represented by an open circle. 8-10 which satisfies none of the applicable void closing criteria. VDCLOS is then exited and not invoked again until the follow- ing subcycle of INFACIi. S 8-11 CHAPTER IX ERROR CONDITIONS The following sections should help the user diagnose a variety of error conditions which may occur during a HELP calculation. The error messages printed by tlie code are listed alphabetically in Section 9.1, along with the sub- routines which print each message and comments which suggest what the user may do to rectify each problem. Section 9.2 describes the optional diagnostic prints the user may obtain in order to clarify an aspect of a calculation or an error condition . 9.1 ALPHABETIC LIST OF ERROR MESSAGES The table that follows is an alphabetic list of the error messages printed by the HELP code. Each error is classified as fatal (F) or nonfatal (Nf); only fatal errors cause execution to stop. The subroutine line number of the statement whicli prints the error message is indicated in parentheses under the subroutine name. The comments direct the user to the most likely cause of the error condition; however, the user is encouraged to pursue the matter further, as there are undoubtedly situations which are not described here . In many cases the user can correct the error by redefining an input parameter or by coding changes to cell quantities or tracer positions in subroutine INPUT. (See Section 7.4.) The calculation can then be restarted on the dump cycle before the error occurred unless the error occurred in one of the generator routines. However, if the error occurs at the beginning of a cycle, before SPHASE, HPHASE or TPHASE have changed any of the cell quantities, and an error dump has been made as a result of a call to ERROR, the user can correct the error and restart from the error dump. (Error #32 is an example.,) Error Message Type Source Comment 1. ALL SPECIFIC VOLUMES Nf CDT MULTIPLIED BY (147) 2. BLANK CARD Nf CARDS (26) 3. CHECK FIRST RECORD OF F INPUT THE DUMP AND FIRST (199) DATA CARD OF INPUT DECK ERROR CONDITION IN F THETAS THETAS M, MOS, (MSLD (78) (N), N»l, MOS) ERROR CONDITION - F THETAS THETAS: THETA, XTP, (115) XBT, YRT, YLF 6. ERROR EXIT - SEE F ERROR STATEMENT NUMBER (9-36) IN 7. ERROR IN ADDTCR. Nf ADDTCR NUMBER OF TRACERS (170) IN PACKAGE EXCEEDED NTTOTI 8. ERROR IN DETIME CELL F DETIME 1= J- HAS NO DETONATION (131) TIME, K,N« Normal iteration pro- cedure not sufficient to conserve volume of cell . A blank card is in the set of cards read by CARDS. Indicates user is read- ing wrong tape or has asked for wrong cycle or problem number on tape . Flags (MSLD) definedby PTSAV are in error. Pos- sibly a storage error. Slipline does not cross this cell and yet PTSAV has indicated that it does. Possibly a storage error. Refer to coding near the statement number in the indicated subroutine. To activate ADDTCR on subsequent run, change NTPMX and reset NADD. Tracer particle arrays may need to be enlarged. A detonation time has not been calculated for a cell containing explosive. Check detona- tion points. Error Message Type Source Comment ERROR IN DETIME ROUTINE, SEARCHING FOR K- KW WAS IN- CREASED TO A VALUE .GT. KMAX OR .LT. I 10. ERROR IN INPUT CARD, F IDET= INPUT CARD IS 11. ERROR IN MFLAG DETER- F MINATION 1= J= K= 12. ERROR IN PTSAV I,J, F Ml, M2, MT, THETA (MT) DETIMfi Indicates the search (425) for the K-th cell is outside the grid. Can result from machine round off error, giving an erroneous value of I or J . DETIME Indicates the type of (444) initiation point (IP) read is in error. Attempted to read primary IP data into the secondary IP array. Input cards are probably out of order. FILGRD Usually a sign that (134) MFLAG(2) was defined incorrectly in the in- put deck. May also be an error in material interface definition. PTSAV Slide line flag (MSLD) (35) for cell I,J improperly defined. Probably a storage error. 13. ERROR IN VOLUMES 1= Nf FILGRD J= MFK= VSUM= VCELL= (31) Sum of material partial volumes do not equal the cell's volume. Can prob- ably ignore if error is very small and occurs for only one or two cells. Can be a result of inter- faces being incorrectly defined by the input cards. Check especially that cards defining tracer particle positions for different materials adjoin- ing a common interface segment have the same end points and define the same number of tracers. ERROR ON PRECEED- ING DATA CARD EVACUATION OF NIATE- RIAL N I,J,N,M,ML, MB TFLUX, XMASS(N,M) , SAMMP(N.M), SAMPY(N.M) FRACS DETECTED AN IN- COMPLETE SUBPACKAGE OF MATERIAL PACKAGE CHECK REMAINING PA^ 'IA.GES Ml, M2, NN- , TXi, TYl, TX2, TY2 = IMAX= WHICH IS NOT AN EVEN NUMBER. THE GRID WAS NOT REZONED IN THE X-DIRECTION IMPROPER GRID DEFINI- TION IN SETUP. IMAX = I *= IMPROPER GRID DEFINI- TION IN SETUP. UMAX - J = INSUFFICIENT AMOUNT OF MIXED CELL STORAGE 21. UMAX * WHICH IS NOT AN EVEN NUMBER. THE GRID WAS NOT REZONED IN THE Y-DIRECTION 22. MASS EVAPORATED DUE TO ROUND-OFF IN TPHASE I,J,MFLAG,AMX,T.E. 23. MASS EVAPORATED DUE TO ROUND-OFF IN TPHASE, I,J,N,MFK,XMASS,SIE,RHO CARDS C28) DMADJ (126\ I 145/ FRACS (562) REZONE (27) SETUP (149) SETUP (151)- NEWFLG ( 6 ) REZONE (156) TPHASE (672) TPHASE (700) Probably numbers are punched in the wrong column on a card read by CARDS. Not an error condition unless XMASS(N,M) is relatively large for evacuation into a neigh- bor cell. Can be caused by interfaces crossing; look for error 31 . NMP(N) not properly de- fined or dummy endpoint not added when tracers were generated or manually changed. Grid must have an even number of columns to be rezoned in the X- direction. User did not define IMAX values of DX. User did not define JMAX values of DY. Code wants to generate more than NMXCLS inter- face cells. Redimension material arrays and re- define NMXCLS. Grid must have an even number of rows to be rezoned in the Y- direction. Mass of a pure cell is too small. This should not occur frequently. Mass of material N in an interface cell is too small or could not be evacuated. Not an error condition unless XMASS is relatively large. 9-4 Error Message T Type Source Comment Nf 26. MIXED CELL MASS TRANSPORT ADJUSTED TO PREVENT OVER- EMPTYING I,J,M,N, SAMPY , SAMMP , SAMMY , SGAMC NEW SLIP ENDPOINT NOT DEFINED L, SLOPS, SLOPT,TXCL,TYCL,TXl , TYl, TX2, TY2= ERROR IN PLGTCR NEW SLIP ENDPOINT F OUTSIDE GRID SLPNDX , SLPNDY, ALPHA, TXCL, TYCL. ERROR IN PLGTCR DMADJ (59) PLGTCR (118) PLGTCR (23) 27. NUMBER OF PARTICLES F TSETUP i HAS BEEN EXCEEDED (156) f i 28. NUMBER OF SLIDE CELLS F PTSAV i EXCEEDS INPUT VALUE OF NSLD (69) ' 2Q 1 . REZONE ROUTINE CAN NOT HANDLE PLUGGING PROBLEMS F REZONE (9) ' 30. SUM OF FRACS NOT Nf FLGSET [ EQUAL TO TOTAL AREA I,J,TAU,XDY2PI = (191) 9-5 Not an error condition unless it recurs for the same cell on many successive cycles. Plug surface extended into free surface cell but can't find an in- tercept with back of target. Check tracer particles of both packages . Plug surface (slipline) can not be extended to a point outside the grid. Check the angle ALPHA. Check input defining tracer particles as well as value assigned to NTPMX. Redimension slide arrays and redefine NSLD. Plugging calculations do not, in general, require rezoning . Interfaces may have crossed due to tracers becoming too sparse. If print occurs only once for a given cell, the error is caused by roundoff and can be ignored. Otherwise back up and move tracers i manually in INPUT or have ADDTCR add tracers in that region. (See Section 7.4.) * (This error can also occur if the interfaces are in- itially defined incorrectly.) ...„«i^,.^. ,i> .. i i i i M.ww ^ *f#.W«S»<SOTffi5 Error Message Type Source Comment 31. TROUBLE DEFINING MATERIAL OF CELL THAT HAS BECOME PURE IN FLGSET, I,J » 32. TROUBLE WITH PRES- SURE ITERATION I = J = ITC = PAV = FLGSET This can be a result of (145) interfaces crossing. Look for error 31. It also can result from an isolated mass. The user can "evaporate" this mass manually in INPUT on a restart. CDT This can be a result of /142 \ too small a value of \229 / PMIN, or it can be caused by a discontinuous equation of state or one which has a negative slope in the P-V plane. I 9.2 INTER DIAGNOSTIC PRINTS Additional diagnostic prints can be generated by the user by appropriately defining the input variable, INTER. The table below indicates the nature of the diagnostic prints for each value Ca prime number) of INTER. If the user wants more than one of the prints, set INTER equal to the product of the corresponding prime numbers. For example, if total energy is not being conserved, the user might set INTER = 5 x 7 to see which cellCs) are responsible for the error. In order to * also see what the cell, quantities are before and after each of the phases the user would set INTER =3x5x7. VALUE OF INTER PRINTS FROM SUBROUTINE NATURE OF INFORMATION PRINTED 2 CDT . All intermediate steps in pressure iteration for each multimaterial cell in active grid. 3 EDIT An EDIT print of all cell quanti- ties after both SPHASE and HPHASE, as well as the normal EDIT follow- ing CDT (on EDIT print cycles only). 5 TPHASE Mass transport terms at each boundary of every cell, plus total energy in grid after each cell in active grid is updated in TPHASE. 7 SPHASE . Total energy in grid after each cell in active grid is updated in SPHASE. 11 UVMOD Post-TPHASE material velocities in slip cells before and after the slipline conditions are invoked. 13 UVCALC Post-HPHASE material velocities in slip cells before and after the slipline conditions are invoked. 9-7 •;nw'.W'Vn''.'*'*i y i ^J» m V "*w>«rrF««'Hf>rirt(0«a>*M« -;<»-r*^*“^V*l'*-«' NOTH: The user should be aware that invoking several of the above options and running a calculation for many cycles can generate a prodigious amount of output, particularly if there are many mixed cells and/or a large grid. Frequently, especially in the case of pressure iteration problems or an energy error, it is necessary to invoke the diagnostic prints for only one cycle in order to get the necessary information. Vv. (The reverse 9-8 of this page is blank)