Skip to main content

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)