SALES 2 - Applied Modelling & Computation Group
Contents
1. t 1 00000E 05 cycle 1 SALES 2 v2 1 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face pressure min 5 88418E 06 max 5 30000E 01 L 5 29947E 02 H 4 76999E 01 dq 5 30006 02 contours for density density SALES 2 v2 1 min 2 90000 t 1 00000E 05 cycle Fri Aug 2 13 07 50 2002 E 03 max 3 00000 E 03 L 2 91000 Example Problem 1 E 03 H 2 99000 Collapsing Cliff Face E 03 dq 1 00000 E 01 zone plot t 1 00021E 01 cycle 3431 SALES 2 v2 1 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face velplot A a E e Sa oe e lai 024 S K N e Reg e ET Ro GL BG Re ai GE eg Pa TE A KV AN GE so EL N n i OE Ri JA d J Sf d of of d vo e Kne Tia PR NE caro EN A AAA 2 01 vmax 2 43666E Collapsing Cliff Face 3431 01 cycle t 1 00021 Fri Aug 2 13 07 50 2002 Example Problem SALES_2 v2 1 contours for pressure 4 pressure SALES 2 v2 1 min 12. 14 14 14 15 16 18 18 19 19 21 21 23 23 24 5 Description of parameters 29 5 1 Input quantities one De eean DE ee ee 29 5 2 Derived and scalar quantities in the common block 32 5 3 Cell and vertex variables a 34 6 Compilation running and a case example 36 6 1 Compilation issues and portability 36 6 2 Example calculation o oo 36 Chapter 1 Introduction Hydrodynamic computer codes hydrocodes are powerful numerical tools for studying the dynamics of continuous media Numerous hydrocodes exist and have been applied to a vast range of problems from simple fluid flow to complex impact cratering simulations involving phase changes and mul tiple materials However there is no such thing as a universal code each hydrocode is developed to study only a small subset of the infinite number of possible problems Thus careful attention must be paid when choosing the appropriate hydrocode for a given problem SALES 2 is an extensively modified version of the Simplified Arbitrary Lagrangian Eulerian computer program SALE Amsden et al 1980 The original SALE code was developed to study only viscous fluid flow for ex ample fluid flow through a channel or pipe SALES_2 has retained all the functionality of the original SALE code however it can also model elas tic and plastic deformation and tensile failure Problems that have been addressed by SALES_2 or a close facs
3. f f 7 PR o JA M M P o p i j p A i l f f p A J M gt P o ENN SA e Es ls E of f P P a P f i i P o I p o f o A ve NN ES 2 v2 1 LU Le DAR E Li RR RER LENS EE RE pa n BD e bag Li 1 EE A cy Ne ON g ES RT ee ag al E 8 FS i ei ad 1 5 CERA T ra Ea gt 0 EE RR n ea aa ca ball do i e t kle A TT a TE Ry eee OER N OON N s de e A _ A A a aaa AI GE N EEN EN 2 PRR ES EE N ON a oe 2 at t 2 00003E 01 cycle Fri Aug 2 13 07 50 2002 6819 vmax 2 28753E 01 Example Problem Collapsing Cliff Face contours for pressure SALES 2 v2 1 min 1 04018E t 2 00003 Fri Aug 2 13 07 50 2002 08 max 8 64073 pressure E 07 L Example Problem E 01 cycle 6819 8 49752E 07 H 6 73648 Collapsing Cliff Face E 07 dq 1 90425 E 07 contours for density density SALES 2 v2 1 min 2 88338 t 2 00005E 01 cycle Fri Aug 2 13 07 50 2002 E 03 max 3 00167 E 03 L 2 89521 Example Problem 6819 E 03 H 2 98985 Collapsing Cliff Face E 03 dg 1 18293 E 01 contours for int energy NA t 2 00003E 01 cycle 6819 SALES 2 v2 1 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face int energy min 1 09566 E 05 max 1 91378
4. SAL contours for pressure pressure SALES 2 v2 1 min 9 50320 t 3 00012E 01 cycle Fri Aug 2 13 07 50 2002 E 07 max 8 21996 10738 Example Problem E 0 EEA 3089 E 07 H 6 44764E Collapsing Cliff Face 07 dq 1 77232 E 07 contours for density density SALES 2 v2 1 min 2 85950 t 3 00012E 01 cycle Fri Aug 2 13 07 50 2002 E 03 max 3 00169 E 03 L 2 87371 10738 Example Problem E 03 H 2 98747 Collapsing Cliff Face E 03 da 1 42193 E 01 contours for int energy SALES_2 v2 1 int energy min 4 62173 t 3 00012E 01 cycle Fri Aug 2 13 07 50 2002 E 05 max 5 17730 E 05 L 3 64182 10738 Example Problem E 05 H 4 19740 Collapsing Cliff Face E 05 dq 9 79903 E 04 zone plot t 4 00020E 01 cycle 15370 SALES 2 v2 1 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face velplot a 2 2 a 2 Le 01 Collapsing Cliff Face 1 82719 xample Problem vmax 15370 01 cycle Fri Aug 2 13 07 50 2002 4 00020 t S 2 v2 1 SAL contours for pressure t 4 00020E 01 cycle 15370 SALES 2 v2 1 Fri Aug 2 13 07
5. E 05 L 7 94715 E 04 H 1 61283 E 05 dq 3 00944 E 04 zone plot fe eo one Pees Ts PaE ui reni aa pe E i ME Aa DR EE RE EE A ee FA MENG EE AA TI t 3 00012E 01 cycle 10738 SALES 2 v2 1 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face A i nio gt R pva EE NE JER o o e o gt deu 7 ni J i in 2 LS Y ES a E e wa 4 J J s N N ea e y of J J 2 SE si velplot f N o AN oe CE e os dl of FR a e dd e at o J of e oe A o SS s of 5 of ES er o JA SS a eo DI of n SL oy 9 o o SL e e o o xd ae SE o os e o Y d o o oy Y o o EVA CI o e e i Vg ES e n e o e e ear o KA AK o e o rn Y e fe a e ut o I Pa o x gt e ea o Ra o o p tap o KA o 7 o od p o ee kad id e o m o o o 2 e o o Ex KA s o ao oe amp amp o o e 5 o _ de o o Sag o e kg o Di a o o o a i R E e e o o o Ki 7 e 3 a pS o o R KN LED o d KJE ee KA LA o o o A es 2 E o gt o o o o o o o o o o o o o o KS o o o o o o o o o o o e o o o o a o o 2 2 2 2 2 01 Collapsing Cliff Face 10738 vmax 2 07053 t 3 00012 xample Problem 01 cycle Fri Aug 2 13 07 50 2002 S_2 v2 1
6. Ou u 4 15 Tey 32 I cyl 2n amp AD in which A is the bulk viscosity coefficient and D is the velocity divergence The coefficient cyl is zero for Cartesian coordinates and unity for cylindrical coordinates By setting the yield strength at zero temperature y0 to zero or by setting the shear modulus mu to zero the effects of elasticity and strength which will be discussed next may be ignored In this case SALES_2 behaves in the same way as the original SALE modelling a Newtonian fluid 4 3 2 Elasticity and the Kelvin Substance The deviatoric strain may be obtained simply by multiplying the deviatoric strain rate by the time step dt thus Era Ge edit Gy 3 eat 4 16 dy E bea The elastic deviatoric stress tensor may then be calculated from Hooke s law namely 23 IL 21 4 17 elastic where u is the shear modulus Note that in contrast to the original SALE manual 7 represents shear viscosity and u represents the shear modulus In the solution algorithm of SALES_2 the total deviatoric stress tensor is the vector sum of the elastic and viscous contributions Il 2ue HE 4 18 Thus the SALES_2 hydrocode models a solid material as a spring and a dash pot in parallel a rheological model known as a Kelvin substance If a constant deviatoric stress S is applied to a Kelvin substance unstrained at time t 0 then the deviatoric strain at time t is given by e
7. Fig ure 4 1 illustrates two such patterns the bow tie pattern and the herringbone pattern Typically vertex coasting is associated with the shortest resolvable wavelengths 2dx in the mesh a b Figure 4 1 Examples of degenerate mesh deformations a the bow tie pattern b the herringbone pattern To prevent such deformations from slowly degrading the solution al ternate mesh nodes can be coupled with a small artificial restoring force Amsden et al 1980 Ideally this force should affect flows only at the 2dx wavelength level to guard against degenerate mesh instabilities but have no influence on the larger flow variations being studied In SALE small accelerations are introduced at each vertex which are based on the sur rounding velocity field and tend to keep the vertex velocities from deviating too strongly from adjacent node velocities As discussed by Amsden et al 1980 a fourth order coupling scheme is effective for the bow tie mode but a more diffusive second order scheme must be used for the herringbone pattern The fourth order form used in SALE is given by Anc i _ gd UW Ut j PL j 1 tt i 1 Pd ul Uuj Uj Upi T Up TW UH 4 16 du in which ane is a coefficient that governs the amount of coupling The second order form is given by wi ul ala af 43 These equations are implemented in SALE in such a way as to enable both order terms to be applied The coeff
8. and y0 As explained in section 4 3 this parameter is the Hugoniot elastic limit for the material at zero temperature By setting this parameter to a high value the crushing strength of competent rock for the fourth layer and zero for the other three layers the upper layer is defined as rigid and the lower layers fluid as though they have failed under the weight of the upper layer For the example problem the fragmentation algorithm is disabled For information regarding the use of this algorithm the user is referred to Melosh et al 1992 The accompanying pdf document example pdf and output file exam ple out are for the example calculation described above and are included for aiding in code verification at other installations 37 Bibliography Amsden A A Ruppel H M and Hirt C W SALE Simplified ALE Computer Program for Fluid Flow at all Speeds Technical Report LA 8095 Los Alamos National Laboratory Los Alamos NM 1980 Anderson C E 1987 An overview of the theory of Hydrocodes Interna tional Journal of Impact Engineering 5 33 59 Ivanov B A Deniem D and Neukum G 1997 Implementation of dy namic strength models into 2D hydrocodes Applications for atmospheric breakup and impact cratering Int J Impact Engineering 20 411 430 Melosh H J Ryan E V and Asphaug E September 1992 Dynamic Fragmentation in Impacts Hydrocode Simulation of Laboratory Impacts JGR 97 No E9
9. Ko and 1 B a b S 1 2 A 2 In the expanded state n lt 1 when the internal energy exceeds the energy of complete vaporization E gt Fey the pressure is given by bpE p ape 2 Bloleo 1 a 0 p0 1 i DEERE E 8 where a and are constants that control the rate of convergence of this equation to the perfect gas law In order to make the transition between the two regimes smooth it has been found that in the partial vaporization regime when p po lt 1 and Ey lt E lt Ew the pressure is best computed from a hybrid formula E Eiw Ph Eev E pe 4 9 5 Ew EV E where pp is computed from the expanded or hot Tillotson equation and pe is computed from the compressed Tillotson equation When high pressure phase transformations occur in a shocked solid the internal energy E computed from the Hugoniot equation includes the transi tion energy E from the low to high pressure phases The internal energy 20 Enpp in the Tillotson equations for the high pressure phase must be cor rected for this transition energy by using the following relation to compute the internal energy p Voo V 2 where Voo is the specific volume of the low pressure phase Vog is taken to be zero at standard temperature and pressure Ehpp re Etr 4 2 5 SALES equation of state algorithm In SALES_2 the equation of state is defined in the subroutine eos Here the subroutine has been modified to repre
10. SALE however if the time step is less than dtane then the amount of coupling is reduced such that the absolute amount of coupling is the same To ensure equivalence in the absolute coupling effect for two calculations using different mesh resolutions an identical time step must be used in each run If this is the case the factor dt dtanc which is inversely proportional to the cell dimension will be smaller for the calculation with the lower mesh resolution Thus the amount of coupling in the lower resolution calculation is reduced compared to that used in simulations with the higher resolution by the ratio of the cell dimensions It is important however that the cho sen time step is that which is appropriate for the highest mesh resolution smallest cell size used This is because the time step restrictions employed by SALES_2 will not allow a large time step which might be appropriate for a low resolution simulation to be used in a higher resolution simulation that is dt dtane cannot be greater than 1 Furthermore the user should be aware that the value of an defined in the input file will only be applied when the time step is equal to the time step as defined by the Courant condition at the beginning of the calculation dt dtane 1 Consequently finding the appropriate value for ane should be carried out using the mesh resolution likely to be used for most of the simulations 17 4 2 Equation of state Before discussing th
11. be sure that the code is reading what the user has input Many errors have been made in the past because this simple step has been overlooked 3 1 2 Start information jname This variable is an eighty character string for the problem title This string will appear at the top of any output file and in the title of any graphics plots The string should be enclosed in single quotes ibegcont This integer variable defines whether the problem is a new prob lem ibegcont 0 a continuation of an old problem 1 or a restart of a previously terminated problem 2 If the problem is a new problem the entire input file will be read If however the problem is a continuation or a restart the input file will only be read until after the time information labelled in the input file The rest of the input variables as well as de rived variables and cell and vertex quantities are read from the dump file by the subroutine tape read The only difference between a continuation and a restart is that for a restart the problem time and cycle number are reset to zero instead of resuming at the problem time that the dump file was made ndump This integer variable defines the number of the dump file from which the input parameters and cell and vertex quantities will be restored 3 2 Improvements to output routines 3 2 1 Plotting routines SALES 2 uses the FORTRAN77 graphics libraries of PGPLOT to produce some basic graphical representation of various ve
12. comed Please contact the author garethQlpl arizona edu Contents Introduction Overview 2 1 Problem set up e 2 2 Cycle hnitializgation i meca eske vo sleep ket 2 3 Lagrangian computation rev vr rv rn rann 2 4 Cycle completion rn nen Improvements to the I O routines 3 1 Problem set up improvements 3 11 New input Ale is sani ie A ane ee 3 1 2 Start information 2 aar vr vr rn nrk 3 2 Improvements to output routines 3 2 1 Plotting routines 204 3 2 2 Monitoring tracers e ss ne 3 2 3 Output data files and dump files Modifications to the SALE solution algorithm 4 1 Problem Set up colo du DL a Fe gt 4 1 1 Mesh generation and multi material capability 4 1 2 Resolution ss atea ska hb eed ea a ad 4 1 3 Alternate Node Coupler 4 2 Equation of state rv vr rv vr rn rn 4 2 1 The liquid equation ofstate 4 2 2 The perfect gas equation of state 4 2 3 The Mie Grunsien and Murnaghan equation of state 4 2 4 The Tillotson equation of state 4 2 5 SALES 2 equation of state algorithm 4 3 Constitutive model 0 2 2220 4 4 3 1 Viscosity and Newtonian fluid flow 4 3 2 Elasticity and the Kelvin Substance 4 3 3 Plasticity Yield Strength and the Bingham Substance 4 3 4 Grady Kipp fragmentation algorithm
13. kg m double node 5 O double node specific internal energy J kg double node n O imat material type identifier integer ey double node ey doublefnode cet double node pivy donble node th double node Tri Trag tac alr simin Pimax signal double node sigma double node B 3 kg double node Pa 3 i i sgmad double node umom VEO double node double node EE ii nom vmonp double node double node node i Hi O artificial viscous pressure Pa double node double node ET i node pixy double node i node 34 TTT TTG fs radial coordinate GELD for cartesian r x Tor cylindrical donbletnode fi direction vertex vosy donbletnode double node double uode m double uode side a a mv mv KE ETT double node i ue ve 35 Chapter 6 Compilation running and a case example 6 1 Compilation issues and portability As mentioned in the introduction SALES_2 has been written almost exclu sively on a Sun Unix platform and compiled using SUN F77 This compiler has the disadvantage of hiding a multitude of sins that more rigorous compil ers may complain about Such issues have lead to changes to version 2 0 of SALES_2 please ensure that you have downloaded version 2 1 before attempting to compile on a different platform If you encounter any compiler errors please report them to the author email garethQlpl a
14. one should remember that while there are nx x ny cells in the mesh there are nx 1 x ny 1 vertices The number of materials required is set by the parameter numat in the input file A maximum of four materials is permitted in the current configuration The input file also requires a list of material descriptions a string of 32 characters to label the materials to be used The remaining material input parameters require a value for each material selected These parameters are discussed in the relevant sections below or in the original SALE manual 4 1 2 Resolution The choice of resolution in space and time is an important issue when using hydrocodes simulation should be conducted with a high enough resolu tion to resolve all the important flow variations in space and time However these needs must be balanced by the specifications of available hardware and the time available to run the simulation The amount of memory required for a given simulation is roughly proportional to the number of cells in the mesh and the duration of the calculation is dependent on the processor speed and more importantly the size of the time step The choice of cell size is both hardware and problem specific and is left to the user to decide The time step size however is often restricted by certain criteria some of which are discussed in the original SALE manual good rule of thumb however is to use a time step that is as large as possible giv
15. phasel If the calculation is in the very low speed regime or incompressibilty is assumed then the pressure forces may be determined by a Newton Raphson iteration This method offers a numerically stable means by which pressure signals can traverse more than one cell in a time step Thus time advanced pressures and velocities can be calculated with greater efficiency than using the explicit pressure calculations with a small time step The last calculation in this section of the code is the energy calculation which is done in subroutine energy Here the contributions of the pressure work and viscous dissipation are combined to compute the change in the internal energy of the cell 2 4 Cycle completion Depending on the type of calculation i e Lagrangian Eulerian or other this section of the code modifies the shape of the mesh or fluxes material through the mesh In the purely Lagrangian case the new grid coordinates are computed in subroutine regrid using the combined velocities computed in the previous section The cell volumes are then modified subroutine volume and the vertex velocities are stored for use next cycle subroutine ultou In the Eulerian case the velocity of the material in a cell relative to the fixed mesh is computed in subroutine regrid The cell centered quantities of internal energy density and momentum are then moved through the mesh to adjacent cells in subroutine advect before the momenta are converted back to v
16. 14735 14759 Melosh H J 1989 Impact cratering geological process number 11 in Oxford monographs on geology and geophysics Oxford University Press 38 zone plot SAL ES 2 v2 1 t 0 00000E 00 cycle O Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face contours for density density SALES 2 v2 1 min 2 90000 t 0 00000E 00 cycle Fri Aug 2 13 07 50 2002 E 03 max 3 00000 E 03 L 2 91000 Example Problem O E 03 H 2 99000 Collapsing Cliff Face E 03 dq 1 00000 E 01 zone plot SAL ESJ2V2 Fri Aug 2 3 07 50 2002 t 1 00000E 05 cycle Example Pro 1 blem Collapsing Cliff Face velplot Fa p DI P p 6 3 6 p 6 6 3 3 4 6 6 3 9 6 3 4 3 P 3 P 4 6 3 6 6 6 3 6 P P 4 6 3 3 P p 6 P 6 3 6 6 3 P 6 6 3 3 P 6 6 3 3 3 o gt o o o gt gt o 9 gt gt o 9 o t 1 00000E 05 cycle 1 vmax 1 00000E 04 SALES 2 v2 1 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face c ntours Tor pressure
17. 15542 t 1 00021E 01 cycle 3451 Fri Aug 2 13 07 50 2002 E 08 max 8 77972 E 07 L 9 52081 Example Problem E 07 H 6 74632 Collapsing Cliff Face E 07 dq 2 03339 E 07 contours for density density SALES_2 v2 1 min 2 89334 t 1 00021E 01 cycle Fri Aug 2 13 07 50 2002 E 03 max 3 00172 E 03 L 2 90418 Example Problem 5431 E 03 H 2 99088 Collapsing Cliff Face E 03 dq 1 08382 E 01 contours for int energy SALES_2 v2 1 int energy min 1 64633 t 1 00021E 01 cycle Fri Aug 2 13 07 50 2002 E 04 max 3 69400 E 04 L 1 11230 Example Problem 3431 E 04 H 3 15997 Collapsing Cliff Face E 04 dq 5 34034 zone plot Map de cessa A I TT ER Ee EE EE mase L ER AE ee IE er VG ae es ae TETEE a a Example Problem Collapsing Cliff Face velplot i er GE ae e N A o i i e M e f f f o Cp us AAA E e ao 2 or f x f N M o ig o I f i
18. 50 2002 Example Problem Collapsing Cliff Face pressure min 9 36437E 07 max 7 52215E 07 L 7 67572E 07 H 5 83350E 07 dq 1 68865E 07 contours for density density SALES 2 v2 1 min 2 83825 t 4 00020E 01 cycle 15370 Fri Aug 2 13 07 50 2002 Example Problem Collapsing Cliff Face E 03 max 3 00170 E 03 L 2 85460E 03 H 2 98536 E 03 dg 1 63455 E 01 contours for int energy SALES_2 v2 1 int energy min 1 08119 t 4 00020E 01 cycle 15370 Fri Aug 2 13 07 50 2002 E 06 max 7 86435 E 05 L 8 94424 Example Problem E 05 H 5 99673 Collapsing Cliff Face E 05 dq 1 86762 E 05 zone plot AL ES 2 v2 1 Fri Aug 2 t 5 00008E 01 cycle 20739 3 07 50 2002 Example Problem Collapsing Cliff Face velplot st 9 Bay SE ye J ee dd A S id e J J one oe A J J e o o 27 Tg A o e DA of d d er 2 jd 9 6 DA mi JER e of or af o siv FA ee add Es e or o a e se E o 32 of ae o o of JA d A VE y Y gt o aa p oe o 75 y A A 4 as Ac se o or v DE a gt o y y ya o ole 2 e e 5 gt N gt o o Y a o e gt GN E og A TE A IS o oe ee ee oe TO o 9 e 2 a o 01 Collapsing Cliff Face 1 62642 xam
19. Also in celset is the algorithm for assigning a material number for each cell The multi material facility does not allow for more than one material to occupy a single cell and hence is limited to purely Lagrangian calculations only The multi material facility included in SALES_2 is very basic as it is envisioned that each user will wish to customize the algorithm extensively The algorithm included in this version divides the mesh into numat equal layers where numat is the number of materials desired maximum 4 layers Any change in the geometry or sizes of these layers should be made by modifying the multi material part of the subroutine celset The algorithm is straightforward The mesh is divided into layers by three boundary lines bndry1 bndry2 and bndry3 depending on how many materials are required any number of these may be superposed and placed at the top of the mesh such that they have no effect on dividing the mesh For example if a two material mesh is desired then bndry1 is set at 14 half way up the mesh while bndry2 and bndry3 are both set at the top of the mesh The program then sweeps over every cell defining a material number imat for every cell depending on whether the cell is below bndry1 or above If all three boundaries are set within the mesh then a maximum number of 4 materials is defined Modifying this algorithm to include more than four material layers or unequal layer heights should be intuitive however
20. LE The i and j corresponding to the cell quantity is that which corresponds to the bottom left vertex of the cell In addition to generating the actual output data files the dlout routine also creates a summary file containing information regarding how many and which cells were traced what properties were traced how many data points were created and so on This information is required by the plotting routines in order to generate the plots of each property for every cell 3 2 3 Output data files and dump files SALES 2 also provides the facility to output all important cell centered and vertex quantities to an output file longprint dat Longprints are activated by setting the ilprt flag to 1 in the input file Longprints also create a file global dat which provides a printout of the total internal energy in the mesh the total energy and the total momentum The problem time interval separating long printouts is defined by the input parameter dtlprt The subroutines which generate these output files are Ingprt and global The last output type is the tape dump These dumps are forced by setting the idump flag to 1 and occur at problem time intervals of dtdump The dump file created contains all the cell centered and vertex quantities for the entire mesh input parameters and derived parameters necessary for restarting the simulation following termination This is a useful feature if the user wishes to continue with an unstable calcul
21. S 24 1 ein 4 19 Thus if the viscosity n is much less than the stiffness u the substance behaves essentially elastically set vise lt lt mu For non zero viscosity however the strain will eventually tend to S 2u This model is the simplest which corresponds to a solid in which vibrations are damped by internal friction Consequently it has frequently been used to model the behaviour of crustal rocks over short time scales The simulations of complex crater collapse discussed in subsequent chapters use this rheological model to rep resent the target under normal pre impact conditions Yt Figure 4 3 Mechanical model of a Kelvin substance 4 3 3 Plasticity Yield Strength and the Bingham Substance The Bingham substance Fig 4 4 may be represented by a spring connected in series with a frictional block and a dash pot in parallel Thus it has much in common with a plastic material The Bingham substance behaves as an elastic body for stresses less than the yield stress above the yield stress the Bingham substance starts to deform permanently However the nature of the flow of a Bingham substance differs to that of a plastic in that flow occurs throughout the material a Bingham substance flows with a viscosity The behaviour of a Bingham substance is described by the equations 24 Oo Oy Figure 4 4 Mechanical model of a Bingham substance oy 2keij if o lt oy 4 20 Cij Oy 2n ij if o gt oy 4 21 To adequate
22. SALES 2 A multi material extension to the SALE hydrocode with improved equation of state and constitutive model Gareth S Collins and H Jay Melosh September 13 2002 Foreword Thank you for downloading SALES_2 This document is a not a complete user manual it merely describes modifications made to the original SALE hydrocode and has been designed to complement the original SALE manual A copy of the SALE manual is available at http www lpl arizona edu tekton sales_2 html SALES 2 is distributed as a number of related files sales 2 f The source code Just over 4000 lines including comments sales_2 comm The common block included by most subroutines sales_2 input The input file set up for an example calculation sales_2 pdf This manual example pdf Graphics output file from an example calculation example out Output file from an example calculation Old users please note that this document is distributed with version 2 1 of SALES 2 SALES2 also requires the PGPLOT graphics libraries These can be downloaded from http www astro caltech edu tjp pgplot The source code is heavily annotated as it is expected that individual users will want to modify portions of it for their own particular problems H J Melosh and his colleagues will attempt to answer technical questions but cannot undertake the task of educating users in how to do hydrocode simulations Enjoy working with SALES_2 Any comments or suggestions are wel
23. ation using a smaller time step continue a previously terminated calculation or use a previous calculation as starting conditions for a new calculation The dump files are unformatted meaning that a human user cannot in terpret them this is to minimize the file size Consequently care must be taken to ensure that the tape dump routines tape_read and tape write are 12 coded such that the order in which each variable is written to the dump file exactly matches the order in which they are read 13 Chapter 4 Modifications to the SALE solution algorithm This section discusses the extensions made to SALE that allow SALES_2 to model a variety of different materials and rheologies First the multi material capability is discussed Then the equation of state and the constitutive model that complete the set of equations solved by SALES_2 are described in depth Finally the fragmentation algorithm is briefly introduced 4 1 Problem Set up 4 1 1 Mesh generation and multi material capability The mesh generating algorithm in SALES_2 generates a rectangular mesh with the origin at the bottom left corner Two optional subroutines have been included sphgrd which converts a rectangular grid into a spherical grid and project which deforms the top of the mesh to simulate a projectile having just made contact with the mesh For more complicated mesh configurations the user will have to modify the algorithm directly in the subroutine celset
24. cold or compressed h for hot or expanded and t for transition 4 3 Constitutive model The original SALE code used a very basic constitutive model which only permitted the study of one rheological model Newtonian fluid flow The constituive model implemented in SALES 2 allows a number of different ma terial rheologies to be modelled Rheological models combine the concepts of elasticity plasticity and viscosity In the aid of understanding each type of behaviour is conventionally represented by a simple mechanical model Elasticity is represented by a spring of stiffness k viscosity 1 is represented 21 by a dash pot and plasticity is represented by a frictional contact of resis tance oy see figure 4 2 Thus a rheological model becomes a combination of these simple elements a b Figure 4 2 Mechanical model analogues for a elasticity b viscosity and c plasticity By defining the input parameters correctly SALES_2 may combine these elements to model three idealised rheologies a Newtonian fluid a Kelvin substance and a Bingham substance The constitutive model is applied in stresd which stands for stress devia tor The stress tensor 0 may be separated into the hydrostatic component p and deviatoric components IL The definition of the deviatoric stress tensor is as follows IL Oij P fori j IL 055 for i j 4 10 where the indices refer to the standard tensor nota
25. e equation of state currently used in SALES_2 it is worth discussing equations of state more generally An equation of state relates three thermodynamic variables pressure p density p or specific volume V 1 p and specific internal energy E or temperature T The equation of state is different for different materials and is a complex function of the molecular and atomic structure of a given substance The response of a material to the shock of an impact is primarily governed by its equation of state Geological materials are too complex for any analytical derivation of their equation of state Therefore the few methods that have been devel oped for describing the equation of state of these materials are based on empirical observations s a result there is still considerable uncertainty in the equation of state of geological materials in the velocity range at which most impacts occur in the solar system Numerical hydrocodes typically compute the material density p and the specific internal energy E every cycle from which the pressure p must be calculated in the next cycle The most convenient form of the equation of state for hydrocode calculations is thus p plp E One should note that most thermodynamic studies derive equations of state in terms of temperature T not E Real physical analysis must use the temperature form and work in this area has lead to ANEOS a semi analytical equation of state However although one can in
26. een added to SALE The routine outputs information about certain cells or vertices in the mesh to separate files for each cell or vertex 10 NB We use IDL Interactive Data Language to interpret the output files produced by idlout Due to the wide use of IDL we have decided to leave this routine in the released version of SALES_2 However we have not released the IDL routines as they are quite problem specific It is envisioned that users of SALES_2 will either disable the idlout routine by setting iidlp 0 or modify this routine to make their own output files for analysis using IDL or some other plotting package iidlp This integer flag defines whether or not the idlout routine is called Set iidlp 1 to enable the dlout routine and iidlp 0 to disable the routine Note however that when disabling idlout all the other information must be consistent That is if ncells is defined as 5 for example there must be 5 sets of i j coordinates defining those cells dtidlp This double precision variable defines the problem time interval between writing to the IDL output files ncells This integer variable defines the number of cells for which data will be written to the output files SALES_2 is currently limited to output data for 10 cells only This may be increased by modifying the idlout routine as directed in the comments nplots This integer value defines the number of variables to be recorded for each of the cells maxim
27. elocities Another type of rezoning may be implemented in this section of the code if required This should be done in the subroutine rezone After this section of the code is completed the code returns to the beginning of the Cycle Initialization stage to begin another time step Chapter 3 Improvements to the I O routines This sections discusses the improvements made to the first two sections of the SALE hydrocode It describes the improved input file and restart procedure plus the modifications made to the output routines both graphical and data When quoted the names of subroutines appear in italic and variable names appear in bold face 3 1 Problem set up improvements 3 1 1 New input file The input file for SALES_2 sales_2 input is heavily annotated each pa rameter is briefly described and options are listed for various integer flags For a description of all the variables used in SALES_2 including the input parameters see section 5 The input file is read by the subroutine rinput Comments in the input file begin with the string x When reading the input file these comments are skipped using the function skip After computing some variables derived from the input parameters the subroutine echo_input writes the input parameters and some derived param eters to the output file to ensure the input file has been read correctly NB It is very important during code development for the user to eramine the output and
28. en the appropriate restrictions In many cases the largest appropriate time step size is limited by the Courant condition on sound signal propagation This time step dteou is the minimum time taken for sound to cross a cell dteou lt min 2 4 1 c C where dx and dy are the cell dimensions and c is the speed of sound in the material modelled The minimum implies that every cell in the mesh must be considered to ensure that dt o satisfies the most restrictive case SALES 2 stipulates that the time step is never larger than one fifth of the minimum sound propagation time dt lt 0 2dteou It is also worth pointing out that when conducting numerical analysis of physical phenomena it is good practice to test the effect of resolution in space and time on the calculation results If an observed result is indepen dent of the density of mesh divisions and the frequency of the sampling time 15 then it is unlikely to be a numerical artifice Such resolution studies have identified a weakness in the alternate node coupling scheme implemented in SALE a modification to correct for this flaw is discussed next 4 1 3 Alternate Node Coupler In a Lagrangian calculation of fluid flow using quadrilateral mesh cells there are certain mesh deformations that do not result in net pressure changes that is there is change in cell shape but no change in cell volume Such nonphysical motion of vertices is sometimes termed vertex coasting
29. erties numat number of materials roi initial density Kg m double numat gt siei initial internal energy J Kg double numat csound speed of sound m s double numat on density of the reference state po Kg m3 double numat tilla Tillotson a parameter tillb Tillotson b parameter a b y atill A linear bulk modulus Pa double numat btill B coefficient of nonlinear elastic term Pa double numat ezero Eo parameter categorizing energy dependence J Kg double numat alfa a constant controlling convergence to perfect gas eos double numat beta B constant controlling convergence to perfect gas eos double numat eiv Ej Specific energy of incipient vaporisation J Kg double numat ecv Ew Specific energy of complete vaporisation J Kg double numat Constitutive model parameters visc linear shear viscosity Pa s double numat lambda linear bulk viscosity Pa s double numat u 0 siem specific internal energy at melting J Kg double numat pelo vb elastic shear modulus Pa double numat hugoniot elastic limit at zero temperature Pa double numat lt Specified inflow conditions density of inflow material kg m double 4 internal energy of inflow material J kg double 4 pap pressure at boundary Pa double 4 velocity of inflow material m s double 4 31 5 2 Derived and scalar quantities in the common block ET double uma dpeo oe s
30. icient xi controls the relative order of the coupling ranging from 0 0 for second order coupling alone to 1 0 for purely fourth order coupling Note that the indices in the two equations above correspond to the coordinates of the vertices within the mesh Although representative values for an and xi are provided with the SALE manual Anc 0 05 xi 0 3 it is essential that the values chosen are appropriate for the particular problem in hand This may be achieved by simulating a known phenomenon under conditions similar to the one being studied and comparing the numerical results with known ones As alluded to above a weakness of the alternate node coupling scheme just described is that it is independent of cell size and time step Thus the velocity diffusion occurs more frequently when the time step is smaller and over greater distances when the cell size is larger This leads to the very undesirable trait of different calculation results being produced for the same problem when modelled using different cell dimensions or time intervals The dependence of the alternate node coupler on time step and cell size has been eliminated in SALES_2 by multiplying the alternate node coupler coefficient anc by the ratio dt dtane where dtane is the maximum time step allowed by the Courant condition on sound signal propagation at the start of the calculation dtane dtcoult 0 When dt dtanc the alternate node coupler is applied as in the original
31. imile of include crater collapse dy namic fragmentation during an impact spallation and atmospheric break up of meteoroids during atmospheric entry SALE is extremely tractable comprising just over 4000 lines of code This fact together with its reliability and stability have resulted in its continued use within the academic community for a wide range of problems Despite its inherent suitability as an introduction to hydrocodes however newcomers to the field are recommended to investigate the subject area before attempting to use the code For an excellent overview of the theory of hydrocodes please refer to Anderson 1987 The modifications and extensions to SALE that are documented in this report make SALES 2 a versatile computer code SALES_2 may be applied to problems involving materials of various types and rheologies and study deformations over a range of speeds and external conditions However like any hydocode SALES_2 cannot be applied to any problem In particular the Eulerian formulation within SALES_2 is restricted to one material only which limits the use of SALES_2 for problems involving extreme deforma tion such as the excavation stage of an impact event It is strongly rec ommended therefore that budding hydrocode users ensure that SALES_2 is the appropriate tool for their problem before embarking on installing or modifying the software This documentation provides only a brief overview of the SALES 2 solu t
32. ion algorithm which is presented in chapter 2 The aim of the manual is to describe in detail the modifications made by Jay Melosh and his co workers to the original version of SALE Users of SALES_2 are thus referred to the original SALE manual to complement this report The SALE manual is available online at http www lpl arizona edu tekton sales_2 html The primary modifications to SALE included in SALES_2 are two fold Firstly improvements to the input and output routines make SALES_2 more user friendly than its predecessor and it has been modified to use the for tran plotting routines of PGPLOT These modifications are discussed in chapter 3 Secondly and more importantly SALES_2 includes a far more sophisticated equation of state algorithm and a much more versatile con stitutive model which includes the Grady Kipp fragmentation algorithm Chapter 4 describes these modifications and where appropriate refers the user to the relevant literature for a more detailed explanation SALE was written in FORTRAN77 a complete listing of which accom panies the SALE manual SALES_2 contains many modified subroutines and some additions also written in FORTRAN77 The complete FORTRAN77 listing for SALES 2 can be downloaded from the SALES 2 website http www lpl arizona edu tekton sales 2 html Included in the download should be a copy of this document the code sales_2 f the common block include file sales_2 comm and a sample input fi
33. ite logical auma 0 ny x cast ogieal mumai one double anco double col double gie double pl double vinn double iil double dus double nx 1 number of x vertices integer ny 1 number of y vertices integer ifirst 1 if wl 1 2 3 or 4 2 if wl 5 or 6 integer jfirst 1 if wb 1 2 3 or 4 2 if wb 5 or 6 integer ilast nx if wr 1 2 3 or 4 nx 1 if wr 5 or 6 integer Sel jlast ny if wt 1 2 8 or 4 ny 1 if wt 5 or 6 integer ifphl 1 if wl 1 2 3 4 or 6 2 if wl 5 integer jfph1 1 if wb 1 2 3 4 or 6 2 if wb 5 integer ilph1 nx if wr 1 2 3 4 or 6 nx 1 if wr 5 integer jlph1 ny if wt 1 2 3 4 or 6 ny 1 if wt 5 integer acoef dpcof fcoef gktrac elaste omeyl anco xicof erfac twplot_ twlprt twidlp twdump ifirst first ilast ilast ifphl jph iphi phil ilastv jlastv ilastm last 32 nm character i character g growing c convection v viscosity m maximum double a double ami double neve integer du integer pmax double vier double dron double ET double dtfrac time step factor for Grady Kipp fragmentation algorithm double nota oge ind double rime double semester Di double double double double crd double with double loops no of dt cuts for non convergence for a given cycle integer dl xl xr yb yt 33 5 3 Cell and vertex variables Cell variables ETER density
34. l No Oo idlout gt tape_write Increment time and cycle number autpsy frstat Compute artificial viscous stresses avisc bc Compute deviatoric viscous stresses stresd bc and apply constitutive model Compute body forces and hydrostatic forces phaset bc If implicit calculate derivatives amp di eos iterate pressures until converged presl bc restep Compute energy changes energy What type of calculation Compute grid velocities Compute material velocity rezone Compute new grid coordinates relative to grid regrid Compute new cell volumes volume Store lagrangian velocities ultou Advect cell centered quantities advect bc Figure 2 1 A general flow chart for the SALES 2 program 2 3 Lagrangian computation It is this section of the code where the effect of all forces on the velocity field is determined In addition the changes in internal energy are com puted First the code computes the effect of any artificial forces applied in the subroutine avisc These artificial forces allow the code to treat the pas sage of shock waves through the mesh and guard against degenerate mesh instabilities See the original SALE manual for further details Secondly the code computes the effect of viscous and elastic stresses in the stresd rou tine It is here that the constitutive model capabilities of SALE have been extended This is discussed in section 4 3 Finally the effect of any external forces and the pressure forces are computed in subroutine
35. le sales_2 input with all the input parameters for a specific test problem We attempt to keep our code as well annotated and clean as possible however our use of the code has been almost exclusively on a Sun Unix plat form so we offer no guarantee on the portability of SALES_2 Compilation issues and an example calculation are discussed in the final chapter of this report Any comments or technical questions are welcome please contact the author Email garethQlpl arizona edu Chapter 2 Overview Figure 2 1 presents the SALES_2 solution algorithm for comparison with figure 4 on page 20 of the original SALE manual Highlighted are the sub routines which have been modified in SALES_2 The algorithm is separated into four sections Problem set up cycle initialization Lagrangian compu tation and cycle completion These sections are discussed below Names of subroutines are in italic 2 1 Problem set up In this part of the code the problem is initialized be it an entirely new calculation or the resumption or restarting of an old calculation First the code processes the input file If necessary 1t then restores previously stored information from a dump file before computing some derived variables used throughout the calculation If the problem is a new one the mesh is then generated It is this part of the code celset where users will need to implement their own mesh generating algorithm for a specific problem Mesh generation and
36. ly model a Bingham substance therefore the constitutive re lation in SALES_2 must also incorporate The approach adopted in SALES_2 to account for plasticity is to reduce the elastic contribution to the stress ten sor when the yield strength of the material is exceeded The yield strength YHEL at zero temperature is defined in the input file by the parameter yO This is then modified to account for material weakening at increased tem perature by the equation E 2 m where F is the specific internal energy and E is the internal energy required for melting Whether or not the yield strength is exceeded in a particular cell is determined by a comparison between the yield strength and an invariant measure of the stress in a cell That is a function of the stress tensor that is independent of the choice of coordinate axes The invariant used is the second invariant of the deviatoric stress tensor J defined as 1 1 Ja GM 3 01 02 02 09 09 ay 4 23 where summation over the indices is implied and 01 09 and og are the principal stresses If J2 exceeds the square of the yield strength then the elastic stress component is reduced in relation to how much the yield strength has been exceeded by The exact amount of stress reduction is given by s fac where ly2 Sfac To 4 24 25 If Jo exceeds the square of the yield strength a EL then the deviatoric stresses are reduced in relation to how much the yield stre
37. ngth has been exceeded by The strains are increased accordingly The exact amount of stress reduction and strain increase is given by fac where y2 TT 4 25 The rheological model available to SALES 2 that most closely analogues a Bingham substance is shown in Figure 4 5 In this rheological model the material behaves as a Kelvin substance at stresses below the yield strength At stresses above the yield stress the material behaves as a Newtonian fluid Thus the deviatoric stress for this pseudo Bingham rheology is given as Ili pei 2m jj if Ja lt Y 4 26 Ii peys fac mp if Ja gt Y 4 27 where m is the effective viscosity of the Bingham substance at stresses below and above the yield stress This rheological model may be invoked by setting the yield strength at zero temperature y0 to the Bingham yield strength of the material to be studied n Figure 4 5 Mechanical model of the pseudo Bingham substance modelled by SALES_2 4 3 4 Grady Kipp fragmentation algorithm Tensile strength is modelled in SALES_2 using the Grady Kipp model of dynamic fragmentation What follows is a very brief description of the algorithm For a far more complete exposition of the physical basis for this model the implementation into the hydrocode and a comparison between the model and laboratory impacts the interested reader is referred to Melosh et al 1992 26 In essence the approach of the fragmentation algori
38. olid 4 2 2 The perfect gas equation of state The pressure in a perfect gas is a simple function of density and internal energy p y 1 pE 4 5 where y is the ratio of the specific heat at a constant pressure to that at constant volume y C Cy 4 2 3 The Mie Grunsien and Murnaghan equation of state The liquid equation of state 4 4 and the perfect gas equation of state 4 5 are sometimes combined to yield a more realistic equation of state for solids called the stiffened gas equation of state P cg p po y 1 pE 4 6 The Mie Grunsien and Murnaghan equations of state are a combination of the perfect gas and liquid equations of state and an extension of the liquid equation of state respectively Although there application in geophysics and planetary science is widespread a more specific equation of state is required for the problem of impact cratering For further reading on these two equations of state refer to Appendix II of Melosh 1989 4 2 4 The Tillotson equation of state This is the equation of state currently used in SALES 2 It was specifically designed for high velocity impact computations in 1962 by J H Tillotson The equation was designed to duplicate the linear shock particle velocity relation at low pressures and to extrapolate to the Thomas Fermi limit at high pressures This equation of state also has parameters that allow it to describe the unloading of shocked material into the vapour
39. ormation the finite velocity of crack growth allows large stresses to accumulate activating smaller flaws Thus the tensile strength in this regime known as the dynamic fragmentation regime is apparently greater than for slow deformation The transition strain rate x at which the dynamic tensile strength exceeds the static strength is given by 4 32 ei m 4 07 4 3 ES m 3 3 m 338 UK This transition strain rate is a sensitive function of the static stress which is in turn a function of the minimum flaw activation strain Emin Ostatic Kemin 4 33 where K is the bulk modulus For a finite body the minimum flaw activation strain depends on the volume of the body V as kemin 1 V which when applied to a spherical body of radius R yields Gran gers ni i 4 34 For most materials m 10 implying a strong inverse dependence of static stress on the size of the body This factor is accounted for in the subroutine celset The fragmentation algorithm in SALES_2 also records the time that fracturing is initiated tfrac the size of fragments fragsz and the ori entation of the principal axes at failure failur in cylindrical symmetry At the time of termination this information is written to the file dam age ouput from the subroutine autpsy The subroutine autpsy also calls a routine frstat to analyse the size frequency distribution of any fragments produced during the calculation The res
40. phase The Tillotson equation has two different forms depending on whether the material is compressed to higher density than its zero pressure form or expanded to lower density The form used in the compressed region p po gt 1 and for cold expanded states where the energy density is less than the energy of incipient vaporization E lt Ey is pE Au Bu 4 7 b PE ENE 1 19 where n p po 4 n 1 and a b A B and Eo are the Tillotson parameters In spite of the notation Eo is not the initial energy density of the substance it is merely a parameter that is often close to the vaporization energy The initial energy density E must actually be zero to ensure that p 0 in the initial state Note that a low density pressure cutoff in cold expanded states for y lt 0 8 to 0 95 must be applied as discussed earlier The parameter a is usually chosen to equal 0 5 to give the correct high pressure limit Note that at high pressure and temperature the pressure tends towards apE and a 0 5 is not what would be expected from the perfect gas equation of state However a 0 5 has been found to fit obser vational data better The linear shock particle velocity model derived from observational re sults from shock wave experiments relates the shock wave velocity U to the particle velocity u by U C Sup The two constants C and S can be expressed in terms of the Tillotson parameters by c VA A is thus equivalent to
41. ple Problem vmax 20739 Fri Aug 2 13 07 50 2002 t 5 00008 01 cycle S V SAL contours for pressure pressure SALES 2 v2 1 min 9 08060 t 5 00008 Fri Aug 2 13 07 50 2002 E 07 max 7 05430 E 01 cycle 20739 Example Problem E 0 46711 E 07 H 5 44081 Collapsing Cliff Face E 07 dq 1 61349E 07 contours for density density SALES 2 v2 1 min 2 82402 t 5 00008E 01 cycle 20 59 Fri Aug 2 13 07 50 2002 E 03 max 3 04481 E 03 L 2 84610 Example Problem E 03 H 3 02273 Collapsing Cliff Face E 03 dq 2 20791E 01 contours for int energy SALES_2 v2 1 int energy min 1 87595 t 5 00008E 01 cycle 20 59 Fri Aug 2 13 07 50 2002 E 06 max 9 55879 E 05 L 1 59277 Example Problem E 06 H 6 72696 Collapsing Cliff Face E 05 dq 2 83183 E 05
42. rizona edu SALES_2 requires the PGPLOT graphics libraries which you will need to install and then link with the SALES_2 object file in order to make the executable For more information about PGPLOT and the libraries themselves the reader is referred to http www astro caltech edu tjp pgplot PGPLOT is available for multiple different compiler platform combina tions however it is worth making sure that PGPLOT is successfully in stalled and operational before attempting to compile SALES_2 6 2 Example calculation The input file supplied with the code is set up to illustrate some of the capabilities of the SALES_2 hydrocode The problem considers a collapsing cliff face The cliff consists of a rigid layer supported by three fluid layers In the example problem the mesh is a 20 x 20 grid made up of 100m x 100m cells Gravity is in the downward direction negative Cylindrical symmetry option is turned off cyl 0 0 To mimic the situation near a cliff face the boundary conditions in this problem are set such that the left wall has no slip allowed the bottom boundary allows free slip and the right and top are free lagrangian surfaces For more information on boundary conditions in SALES_2 please refer to the SALE manual In terms of material properties the mesh in the example calculation is 36 divided into 4 layers The bottom three layers have identical material prop erties The top layer differs only in its value of density roi
43. rtex and cell centered quan tities during run time The vertex positions are plotted using zonplt in lagrangian mode only the vertex velocities are plotted using velplt and various cell centered quantities are contoured using contur NB SALES 2 will not compile or run unless PGPLOT has been installed and linked correctly For more information on PGPLOT see Chapter 6 Before prompting the user for the input and output file names the code will request an output type for the plotting routines The types of out put available will depend on the drivers chosen when installing PGPLOT The standard output types for a Unix system are Xwin Output to an X window on the screen ps Output to a postscript file cps Colour postscript null suppress graphical output The plotting routines are called from fulout The plotting routines are called at regular intervals of dtplot in problem time and at the zeroth and last time step 3 2 2 Monitoring tracers The plotting routines included in SALE which in SALES_2 utilize PG PLOT display properties of the mesh at regular time intervals dtplot during the run It is frequently the case however that the user s interest is not in one particular instant in time but in the behaviour of a particular cell or particle during the entire calculation Such particles or cells are often called Lagrangian tracer particles In order to facilitate the study of tracer particles a routine called idlout has b
44. sent the Tillotson equation of state although the algorithm may be simplified to convert the equation of state to that for a perfect gas stiffened gas or a linear elastic solid By setting a b B Q 16 Eo 0 equations 4 7 and 4 8 reduce to the liquid equation of state 4 4 with A Ko the bulk modulus If instead A b B a 2 Ep 0 then equations 4 7 and 4 8 reduce to the perfect gas equation of state 4 5 with a y 1 To get the stiffened gas equation of state 4 6 therefore set b B a 8 Eo 0 A Ko and a y 1 The call to the subroutine eos passes the density ro ij internal energy sie ij the damage coefficient damag ij and the material type imat ij eos then decides which state the cell is in expansion or compression and calculates the component of the pressure pe or pp using the appropriate equation 4 7 or 4 8 If the material is fractured the compression term of the pressure is modified accordingly see section 4 3 4 If the density of the cell becomes too low p po lt density_cutoff then the pressure is set to zero and a flag is raised density_flag The appropriate pressure component pe and pp is then returned unless the cell is in the transition regime In this case the hybrid formula 4 9 is used The new value for the pressure in the cell is then returned with the density flag low_density ij and a flag identifying the equation of state regime used to calculate the pressure eos iddt ij c for
45. simple expression D nV 4 29 where n is the number of cracks per unit volume and V 4 3r1 is the volume of a spherical stress relieved region surrounding the crack of half length l This accounts for the reduction in stress in the vicinity of a crack The number of flaws activated at any time is given by a two parameter Weibull distribution N ke 4 30 where N is the number of flaws per unit volume activated at or below tensile strain and k and m are material constants defined in the input file by pweib and cweib Once activated the cracks are assumed to grow at their maximum speed cg another input parameter cg and that the number of flaws that will subseguently activate is reduced by the factor 1 D The rate of damage increase is approximated by the expression dt 3 where a 87c k m 1 m 2 m 3 1 3 a N al 3e7 3 4 31 27 Laboratory experiments of fracturing have shown that for slow defor mation rates the largest flaws dominate larger flaws begin to grow under smaller stresses When a flaw begins to grow it releases stresses in the vicinity of the crack tending to prevent the growth of other flaws Conse quently large stresses do not accumulate in the vicinity of the crack and smaller flaws are prevented from growing The result is that under slow deformation complete fracture occurs at a definite strain or stress known as the static stress Ostatic For high rates of def
46. some cases de rive one from the other by fast Newtonian iteration in general there is only an indirect relationship between temperature and internal energy There fore the most widely used equations of state have been developed with the above form 4 2 1 The liquid equation of state The pressure of a liquid subjected to a small compression or expansion is a linear function of the density fluctuation about the uncompressed density Po p chlp po 4 4 where cg y Ko po is the bulk sound speed Although it appears that this equation has no dependence on E this is a false impression as Ko is a weak function of temperature Ko is the isentropic bulk modulus for sound waves and the isothermal bulk modulus for slow static compression These two modulii differ by less than one percent for most substances In hydrocode computations densities much less than po will occur during rapid expansion Physically such densities indicate that the material in the 18 computational cell is no longer a continuum but may be in the form of droplets or solid fragments separated by gas or vacuum instead Therefore the pressure computed from this equation of state must be limited to avoid artificially large tensions developing This is achieved by testing the ratio p po if it is less than some amount often chosen between 0 8 and 0 95 the pressure is set to zero or in sophisticated equations of state to the vapour pressure over the liquid or s
47. the multiple material facility are discussed in section 4 1 1 Finally all the vertex and cell centered quantities are initialized in celset The initial pressures are computed using the equation of state via the eos routine except for the incompressible case for which the initial pressure is zero The subroutine bcset adjusts the boundary values if inflow or pressure boundaries are specified 2 2 Cycle initialization In this section the code advances to the next time step First the size of the time step is chosen according to the same rules as the original SALE The cell pressures are then calculated for the compressible case using eos discussed in section 4 2 Finally outputs are performed when appropriate and the time and cycle number are incremented At the end of the cycle the code returns to the beginning of this section of the code Subroutines Primary Secondary What type of problem rinput New Continuation o Read input parameters restart b rinput Compute derived quantities o Restore variables from dump file E restore tape read Compute derived guantities 2 Echo input and derived quantites to output file A echo input Generate mesh and setup cell variables celset volume sphgrd beset projet bc Calculate time step C timstp 0 Initialise lagrangian guantities g vinit eos acontur velplt Begin cycle Output at appropriate times us Newcyc fuloutg zonplt E Ingprt Time to stop o globa
48. thm is to assess the amount of tensile fracture in a cell and reduce the cohesive strength accordingly The model has been used to study dynamic fragmentation during impacts spallation and break up of meteoroids during atmospheric entry However for large scale impact events the algorithm should be used with caution as it does not describe the system of faults and joints typical for natural rocks In this case any pre event damage of the target should be treated specifically Ivanov et al 1997 The amount of fracturing is calculated on a cell by cell basis Fracturing in the cell is parameterised by the damage quantity D 0 lt D lt 1 D is 1 0 for wholesale fracturing in the cell the material is a collection of unconnected fragments and offers no resistance to further extensional strain D 0 0 corresponds to totally intact material in the cell The fragmentation model is implemented in such a way that the amount of damage affects only the cohesive strength of the material Thus only the shear modulus is reduced for damaged material the elastic modulus of completely damaged material is reduced to the bulk modulus The equation implemented in SALES_ to account for the reduction in the shear modulus is Ilj 2u 1 De 4 28 The fragmentation model assumes that within the continuum material there are points of weakness cracks or flaws which grow under tensile stress The damage is related to the number and size of these flaws by the
49. tion not cell number The hydrostatic component pressure is used in the equation of state to calculate volumetric changes discussed above The deviatoric component is used in the constitutive model to predict changes in distortion The strain tensor may be equivalently separated into a volumetric part e and a deviatoric part j efori j Ej eg for i j eis defined as the volumetric strain or trace and is given by the expression in cylindrical coordinates 1 Similar expressions exist for the deviatoric strain rate The deviatoric stresses may be a function of the deviatoric strain elasticity or the devia toric strain rate viscosity or both In SALES_2 the viscous deviatoric stress tensor and the elastic deviatoric stress tensor are calculated separately 22 4 3 1 Viscosity and Newtonian fluid flow Writing the equation for a Newtonian fluid in terms of deviatoric stress and strain rate gives IL N 4 12 Viscous where y is the shear viscosity In two dimensions the deviatoric strain rate in a cell may be calculated from the velocity gradient across it less the volumetric strain or trace i a Exa Fn dx i dv Eyy Dy e 1 v du If considering three dimensions with axial symmetry then the hoop strain rate p is defined as y u i 4 14 Es Therefore the viscous stresses in cylindrical coordinates may be defined as Mar P HAD My 2 AD dv
50. ults of this analysis are output to the file fragsz output 28 Chapter 5 Description of parameters 5 1 Input quantities problem title problem type flag 0 new 1 cont 2 restart dump file number for restart IDL tracer output flag 1 on 0 off IDL tracer output time interval number of tracer cells vertices number of parameters to be recorded cells to be monitored parameters to be recorded maximum time step see SALE manual problem finish time graphics plot interval complete printout flag 0 off 1 on complete printout time interval file dump flag 0 off 1 on file dump time interval 29 Mesh properties gravitational acceleration components in m s double double double boundary condition flags see SALE manual integer Computational information see SALE manual implicit pressure calculation flag 0 exp 1 imp i incompressible flag 0 compressible 1 inc integer integer rezone flag 0 eulerian 1 lagrangian 2 other relaxation coefficient used in pressure iteration pressure fraction scaling the relaxation factor rdsdp integer double double double double double double relaxation factor for the continuous grid rezoning artificial bulk viscosity coefficient alternate node coupler coefficient 1 4th order node coupling 0 2nd order node coupling double allowed relative error in the pressure iteration convective fluxing parameters double 30 Material prop
51. um of 12 quantities may be recorded 1 Pressure cell 2 Velocity in x direction particle 3 Velocity in y direction particle 4 Maximum principal stress cell 5 Intermediate principal stress cell 6 Minimum principal stress cell 7 Density cell 8 Mass cell 9 Volume cell 10 x coordinate particle 11 y coordinate particle 12 Damage cell 11 i cell j_cell These integer arrays record the i and j locations of the cells to be traced The i and j correspond to the cell quantity and the bottom left vertex of the cell There should be ncells sets of and 5 pairs plot_type This integer array defines the quantities to be recorded for each cell taken from the list above There must be nplots number of elements in plot_type In summary idlout records a number of selected cell centered or particle properties for each cell particle of interest during run time Every interval of dtidlp seconds problem time these properties are written to an output file corresponding to each tracer particle The output files may then be interpreted using various IDL plotting routines The user is required to enter the number of tracer cells particles ncells and the number of properties to be recorded nplots For each cell particle desired the user must enter an 2 j location The user must also specify the properties required chosen from a list of possibilities NB The labelling convention is the same as the original SA
Download Pdf Manuals
Related Search