CC BY-NC-ND
0Available online at www.sciencedirect.com
ScienceDirect Procedia
Engineering
Procedia Engineering 103 (2015) 173 - 180 ;
www.elsevier.com/locate/proeedia
The 13 th Hypervelocity Impact Symposium
Simulations of Defense Strategies for Bennu: Material Characterization
and Impulse Delivery
E.B. Herbolda*, J.M. Owena, D.C. Swifta, P.L. Millera
aLawrence Livermore National Laboratory, PO Box 808, Livermore, CA, 94551, USA
Abstract
Assessments of asteroid deflection strategies depend on material characterization to reduce the uncertainty in predictions of the deflection velocity resulting from impulsive loading. In addition to strength, equation of state, the initial state of the material including its competency (i.e. fractured or monolithic) and the amount of micro- or macroscopic porosity are important considerations to predict the thermomechanical response. There is recent interest in observing near-Earth asteroid (101955) Bennu due to its classification of being potentially hazardous with close approaches occurring every 6 years. Bennu is relatively large with a nominal diameter of 492 m, density estimates ranging from 0.9-1.26 g/cm3 and is composed mainly of carbonaceous chondrite. There is a lack of data for highly porous carbonaceous chondrite at very large pressures and temperatures. In the absence of the specific material composition and state (e.g. layering, porosity as a function of depth) on Bennu we introduce a continuum constitutive model based on the response of granular materials and provide impact and standoff explosion simulations to investigate the response of highly porous materials to these types of impulsive loading scenarios. Simulations with impact speeds of 5 km/s show that the shock wave emanating from the impact site is highly dispersive and that a 10% porous material has a larger compacted volume compared with a 40% porous material with the same bulk density due to differences in compaction response. A three-dimensional simulation of a 190 kT standoff explosion 160 m off the surface of a shape model of Bennu estimated a deflection velocity of 10 cm/s.
© 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Peer-review under responsibility of the Curators of the University of Missouri On behalf of the Missouri University of Science and Technology Keywords: Bennu, constitutive modeling, material characterization, asteroid deflection, Eulerian hydrocode
CrossMar]
1. Introduction
Impulsive Strategies for disrupting or deflecting asteroids with low collision lead-times include two main approaches: standoff explosion and hypervelocity impact. Assessment of these strategies depends on material characterization including its strength and equation of state to attain any certainty in predictions of the deflection velocity. In addition to bulk material properties, the initial state of the material including its competency (i.e. fractured or monolithic) and the amount of micro-or macroscopic porosity are important considerations [1-6]. There is a lack of data for highly porous carbonaceous chondrite at very large pressures and temperatures considering the compositional differences between collected non-porous samples [7-18]. Specifically, very little is known about the material strength at high porosity levels that have been estimated for near-Earth asteroid (101955) Bennu [19]. There is recent interest in observing Bennu due to its classification of being potentially hazardous (i.e. diameter > 150m and MOID < 0.05 AU) with close approaches occurring every 6 years [17]. Bennu is relatively large with a nominal diameter of 492 m, density of 1.26 g/cm3 [20] and is thought to be composed mainly of CI or CM carbonaceous chondrite based on spectral data [21].
The response of an asteroid to impulsive loading for deflection scenarios largely depends on the amount of material ejected from the deposition layer or impact site, respectively. The momentum of ejecta depends on the thermomechanical response, which must also account for the chemical composition of the parent material. CI and CM carbonaceous chondrites are known to have a hydrated mineralogy and there is evidence of shock-induced devolatilization of carbonaceous chondrites between 10-30 GPa [8,14,16,22] indicating a complex thermo-chemical response as well.
1877-7058 © 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Peer-review under responsibility of the Curators of the University of Missouri On behalf of the Missouri University of Science and Technology doi: 10. 1016/j .proeng.2015.04.024
A distinct element modeling (DEM) approach was used to investigate the response of unconsolidated boulders to a standoff explosion in a previous investigation [23] noting that the deflection velocity did not account for ablation of material, but rather focused on the interaction between the large boulders after the impulse was applied. Using Lagrangian finite elements, a similar problem was investigated in [24] where it was found that the morphology of contacting boulders in unconsolidated asteroids affects the pulse propagation speed due to the efficiency of the contacting surfaces to transfer the impulse. The morphological properties of Bennu to date do not indicate its internal structure. Its surface is thought to be mostly covered in a fine regolith powder with sizes ranging below 125 fim [21]. Accounting for the most recent estimations of the composition and structure of Bennu it is probable that it is a loosely consolidated body without strength.
Granular materials, such as the surface regolith typically exhibit a pressure dependent resistance to distortional loading. In the absence of the specific material composition and state (e.g. layering, porosity as a function of depth) on Bennu we introduce a continuum constitutive model based on the response of granular materials here and provide impact and standoff explosion simulations to investigate the response of highly porous materials to these types of impulsive loading scenarios.
The two different approaches differ in the amount of energy imparted to the asteroid as well as the resultant impulse propagation through the body. For relatively large standoff distances, the resulting shock wave will converge in contrast to a divergent wave emanating from the impact point, which means the material will be subjected to different dynamic loading conditions. Both strategies presented here focus on the strengths and weaknesses of each delivery mechanism given the characterization of the intrinsic properties of Bennu and Tillotson equation of state fit to carbonaceous chondrite data at high pressures and temperatures.
2. Numerical Modeling
2.1. Use of Hydrocode for deflection simulations
An Eulerian hydrocode (GEODYN) with an interface reconstruction algorithm, wide-range equation of states and a flexible constitutive model library [25,26] was used for the numerical investigations of both defense strategies (standoff bursts and hypervelocity impacts). GEODYN uses a high-order material interface reconstruction algorithm [27] and advanced constitutive models that incorporate salient features of the dynamic response of geologic media [28,29]. The details of adaptive mesh refinement (AMR) in Eulerian hydrocodes is a relatively mature technique [30] that allows the computational expense to be focused on the parts of the problem of interest. A high order Eulerian Godunov scheme has been shown to provide accurate solutions to problems resolving shock waves. The methodology implemented in GEODYN is based on adaptations of the single-phase high-order Godunov scheme. Further details of this implementation may be found in Appendix A.
2.2. A thermomechanical constitutive model for asteroid regolith
Regolith samples consist of solid material and pore space. The model described here is based on previous mechanical constitutive models developed for crushable soils [31-33]. The constitutive model development evolves from tracking the development of the stresses and energies of the solid and pore space separately. A full description of the model will be detailed separately in a forthcoming manuscript. The pressure and temperature dependent strength model begins with the concept of a porous granular medium where pore space exists between individual grains in a loosely consolidated volume. The change in volume of the material can be expressed as
dv = dv + dv (1)
s p v '
where dvs is the volume change of the solid and dvp=dvpe+dvpf is the volume change of the void space, which is composed of the elastic portion (i.e. poroelastic compaction) dvpe and the void space changes due to fragmentation and comminution of the grains dvpf. The values of each term in Eq. (1) in the reference configuration are denoted by capital letters
dV = dV + dV = dV + dV + dV. (2)
s p s pe pf '
The porosity and its reference value are defined by
dv dv dv, ,
A-_£_ =_EL +_(3)
, dv dv dv
dV dV dl' (4)
dV dV dV
P [s/cm3]
Figure 1. Fit of the Tillotson EOS to shock compression data for carbonaceous chondrite. (1) Fit of the solid EOS to the lower density chondrite classes. (2) Fit of the pressure dependent strength model using the solid EOS parameters. (3) Murchison data from [15]. (4) Bruderheim data from [15]. (5) Serpentine data from [34]. (6) Lizardite data from [22]. (7) H6 chondrite Kernouve data from [12]. (8) Antigorite data from [22]. Note that data from the hydrated chondrites (3), (5), (6) and (8) have a lower bulk density than (4) and (7).
A useful measure of the deformation comes from the definition of the Finger tensor B=FF (i.e. the left Cauchy-Green deformation tensor). The unimodular part of this tensor is defined from [35] b = JFFr, where J = det(F) = dv/dl'is the Jacobian and is a measure of the measure of volumetric compression. The initial and current density of the solid/void mixture is
A =(!-*)/>.. (5)
P=(1~*)P„ (6)
respectively. For large deformation models the symmetric unimodular tensor is a measure of pure elastic distortion through the evolution equation
B;=lb;+B;Lt - : i)Be - TA ,
which derives from [36-38]. The product characterizes the inelastic distortional deformation and Ad is defined:
Aj = B' -
For hyperelastic model evaluation, the full stress tensor derives from the Helmholtz energy function. From [39] the Helmholtz free energy (HFE) is a function of temperature and damage. In [33] the Helmholtz free energy is written for comminuting sand grains and is a function of a "breakage" parameter and porosity. A HFE capable of handling failure, damage/comminution is,
pJ:A6)
(e-e0)-e\n
The time ueiivauve oi me nrn proposeu m nquauon (?) proviues me energy uaiance equanon wnereJi{i>v), J(.&,>;, and w(0) represent functions that account for the nonlinear elastic loading of a granular material (i.e. poroelasticity) as well as the onset of comminution of the grains. Here we have coupled the mechanical breakage model derived from [33] to the Tillotson equation of state through careful selection of f\(ev) and f(ev) such that for large values of solid compression, it controls the response since the material may be fully compacted and, at low compression, the nonlinear poroelastic response controls the solid material response. The Tillotson equation of state was used here to describe the volume response of the solid material whose parameters are given in Table 1. The parameters in Table 1 are fit to CI and CM data for carbonaceous chondrites and related minerals shown in data sets (3) and (5) in Figure 1. The fit of the Tillotson equation of state was performed on three different data sets for chondrites with grain densities around 2.5 g/cm3. Curve (1) in Figure 1 is the solid equation of state fit to the hydrated chondrites. Curve (2) shows the response of a 50% porous sample using the pressure dependent strength model described in this section. The remaining data is listed from various sources for chondrites and similar minerals. The data points in (3) and (4) is for Murchison and Bruderheim carbonaceous chondrites from [15]. Data for Serpentine (5) are from [34], Lizardite (6) [22], (7) H6 chondrite Kernouve [12] and Antigorite (8) [22]
driver (solid) (high temp solid) solidus liquidus
Figure 2. Solidus and Liquidus delineations for H6 chondrite from [7] compared with the model.
are also shown. Note that data from hydrated carbonaceous chondrites (3), (5), (6) and (8) have a lower bulk density than data from (4) and (7).
Table 1. Temperature dependent heat capacities (in imitt of J/g't);
Material EOS P* A B a b or £
CI/CMOiondiite Tillot&on [gem3] [GPa] [GPa]
2.50 24.2 SO.O 0.5 0.47 5.0 5.0
The chemical composition differs between the different mineralogy and so an estimate of the temperature dependent heat capacity was fit to shock induced melt data for carbonaceous chondrites from [7]. The specific heat used for chondrite comes from [7] Cv(Q)=a+bQ+c/Q2 where a=868 J/kg/K, b=0.1788 J/kg/K, and c=-1.83*107 J/kg/K2. The result is shown in Figure 2 where the red and blue lines show the liquidus and solidus boundaries, respectively, and calibrated model results from GEODYN are shown for a material starting at 10 K and 1600 K under uniaxial compression from 0-30GPa.
We also have relations from [33] for the evolution equations for the breakage (i.e. damage or comminution), solid compression, porosity and deformation tensor, B'e. In addition to the model described in [33], we have introduced
the equation of state and yield to account for high pressures. The evolution equations for breakage parameter B, the solid volumetric strain sv, porosity and the first invariant of B[ are:
B = TABU) (10a)
i>—T\{<t>,B)
D-r(I:A,)
(10b) (10c) (10d)
The expressions given in Eq. (10a)-(10d) are coupled and encompass the poroelastic deformation, compaction and bulking response commonly associated with granular materials under distortional deformation. In each Eulerian cell in a simulation, satisfaction of these four expressions requires iterations to find the parameter r, characterizing the extent of inelastic deformation. The yield surface / from [33] is a geologic cap model that accounts for both comminuting granular
material and material yield,y=y(B,s^</>,Bii,7).. Here, the pressure dependent yield surface for the solid material is,
Y = Y,
1 +Myp/Ys0 1 +MyMp (m, <5/[ 1+(m1 <£■„ -
Figure 3. Impact craters are compared between two materials with porosities of 0.1 (left side of each figure) and 0.4 (right side of each figure). The color distribution shown indicates the change in porosity from its initial value (red indicates 0.0). The impactor was a 10 ton steel projectile with an impact speed of 5 km/s.
The shear modulus for the porous materials is also a function of pressure,
G=(1 - a)g(1 - CB) (1 - ^ (1+s)s +ag0s, (13)
which could easily be cast as a function of temperature for future investigations. The parameters used in the current model are given in Table 2. The onset of comminution is a function of the solid compression, breakage and current porosity and changes with the current state of bulking and compaction of the material such that the material response is consistent through repeated loading and unloading conditions.
Table 2. Strength model parameters for the modified breakage model.
Parameter Value Units Parameter Value Units Parameter Value Units
ko 800 [GPa] £vH 0.4 0nax 0.41
go 500 [GPa] iic 0.004 B0 0.0
gOs 0 [GPa] My 0.0 [GPa] ßB 13.0
Ç 0.9 [GPa] Myh, 50.0 [GPa] b 0.99
TsO 2.0 [GPa] 00 0.4 m 0.5
mi 10.0 0mm 0.05
3. Standoff explosion and impact simulations
In [19] the estimated density of Bennu is listed as 1.26 g/cm3 and the porosity is estimated around 40%. However, density estimates range from 0.9-1.26 g/cm3 in recent literature illustrating the challenge of characterizing potentially hazardous asteroids in general. Due to the uncertainty in measurements, it is interesting to investigate the differences between plausible material responses based on the range of measured properties. Here we utilize the constitutive model introduced in section 2, which is capable of handling arbitrarily large porosities and the granular material may transition between a gaseous (when sv=p=0) or solid state depending on the value of jamming porosity, which indicates the unconsolidated density of the material. This model was used in GEODYN simulations of hypervelocity impact simulations where a 10 ton steel sphere was impacted onto a simulated regolith layer at 5 km/s. Figure 3 shows a comparison of impact craters two materials with porosities of </>= 0.1 (left side of each figure) and ij>= 0.4 (right side of each figure). The color distribution shown indicates the change in porosity from its initial value (red indicates ij> = 0.0), to illustrate the level of compaction. The comparison between the two porosity ranges shows a difference in the compaction front illustrated at 5 ms and 10 ms in Figure 3. Both materials have the same overall density of 1.26 g/cm3, such that the grain density is p = 1.26/(1-$)). The impedance differs between the two materials and the shock wave emanating from the impact site is highly dispersive due to the crushing of the highly porous material within the body. The width of the compaction bands shown in the simulations are similar at 10 ms, but it is clear that the 10% porous material has a much larger damaged volume.
Figure 4. Porosity and temperature distribution in 3D simulation with a 190 kT standoff explosion after 1.6 s elapsed. The initial impulse travelled through the porous material, but ceased crushing after traveling about 20m. The temperature distribution shows heated material ablating from the surface. The temperatures shown indicate that a significant portion of the surface material is still above the melt temperature, though the heated layer is less than 10m deep in some locations. This could be due to the fact that porous materials cannot effectively generate pressure upon heating until the material has expanded enough to fill the surrounding void space.
It is also interesting to consider the response of a highly porous asteroid to a standoff explosion. A shape model for Bennu from [40] was used in a three-dimensional GEODYN simulation. The deposited energy on the surface of the asteroid was 190 kT of energy from a distance of 160 m from the +x surface model of Bennu. The porosity and temperature distribution with a standoff explosion is shown in Figure 4 after 1.6 s elapsed. The initial impulse travelled through the porous material, but ceased crushing after traveling about 20 m. The temperature distribution shows heated material ablating from the surface. The temperatures shown indicate that a significant portion of the surface material is still above the melt temperature, though the heated layer is less than 10m deep in some locations. This could be due to the fact that porous materials cannot effectively generate pressure upon heating until the material has expanded enough to fill the surrounding void space.
Figure 5. The time history of the deflection velocity is shown for the center point of the simulated model of Bennu. The inset view is cut away by two planes and this point is shown at t = 1.6 s. The time history indicates that there is a rather sharp acceleration away from the standoff explosion followed by a pull-back and finally a gradual rise to the final deflection velocity of 0.1 m/s. The deposited energy on the surface of the asteroid was 190 kT of energy from a distance of 160m from the +x surface model of Bennu.
The time history of the deflection velocity is shown for the center point of the simulated model of Bennu. The inset view is cut away by two planes and this point is shown at t = 1.6 s. The time history indicates that there is a rather sharp acceleration away from the standoff explosion followed by a pull-back and finally a gradual rise to the final deflection velocity of 0.1 m/s. The velocity distribution shows a relatively late time residual velocity gradient across the entire body. The escape velocity is approximately 0.2 m/s for Bennu.
4. Conclusions
Bennu is a potentially hazardous asteroid is thought to be composed mainly of CI or CM carbonaceous chondrite covered in a layer of regolith. Impulsive strategies for disrupting or deflecting asteroids were investigated focusing on standoff explosion and hypervelocity impact. In the absence of the specific material composition and state on Bennu we introduced a continuum constitutive model based on the response of granular materials and provided impact and standoff explosion simulations to investigate the response of highly porous materials to these types of impulsive loading scenarios. The two different approaches differ in the amount of energy imparted to the asteroid as well as the resultant impulse propagation through the body. The width of the compaction bands is similar at 10 ms for hypervelocity impact simulations, but it is clear that lower porosity material has a much larger damaged volume. A shape model for Bennu was used in a three-dimensional GEODYN simulation. The initial impulse travelled through the porous material, but ceased crushing after traveling about 20 m. The temperatures shown indicate that a significant portion of the surface material is still above the melt temperature, though the heated layer is less than 10m deep in some locations. This could be due to the fact that porous materials cannot effectively generate pressure upon heating until the material has expanded enough to fill the surrounding void space.
Acknowledgements
This work was funded by the Laboratory Directed Research and Development Program at LLNL under project tracking code 12-ERD-005, and was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. EBH would like to thank M.B. Rubin for his contribution to the development of a similar constitutive model as described here.
References
[1] Lomov I, Herbold EB, Antoun TH, Miller P. Procedia Eng 2013;58:251.
[2] Sanchez P, Scheeres DJ. arXiv:13061622v2 [astro-phEP] 2014:1.
[3] Britt D, Consolmagno G. Icarus 2001;152:134.
[4] Britt DT, Yeomans D, Housen K, Consolmagno G. Asteroid Density , Porosity , and Structure, in:. Bottke Jr. WF, Cellino A, Paolicchi P, Binzel RP (Eds.). Ateroids III. Tucson: University of Arizona Press, Tucson; 1999.
[5] Benavidez PG, Durda DD, Enke BL, Bottke WF, Nesvorny D, Richardson DC, Asphaug E, Merline WJ. Icarus 2012;219:57.
[6] Consolmagno GJ, Britt DT, Macke RJ. Chemie der Erde - Geochemistry 2008;68:1.
[7] Asahara Y, Kubo T, Kondo T. Phys Earth Planet Inter 2004;143-144:421.
[8] Tyburczy JA, Frisch B, Ahrens TJ. Earth Planet Sci Lett 1986;80:201.
[9] Tomeoka K, Yamahana Y, Sekine T. Geochim Cosmochim Acta 1999;63:3683.
[10] Young ED. Philos Trans R Soc London Ser A 2001;359:2095.
[11] Horz F, Cintala MJ, See TH, Le L. Meteorit Planet Sci 2005;40:1329.
[12] Schmitt RT. Meteorit Planet Sci 2000;35:545.
[13] Tagle R, Hecht L. Meteorit Planet Sci 2006;41:1721.
[14] Lange MA, Lambert P, Ahrens TJ. Geochim Cosmochim Acta 1985;49:1715.
[15] Anderson WW, Ahrens TJ. Shock Wave Equations of State of Chondritic Meteorites, in:. Schmidt SC, Dandekar DP, Forbes JW (Eds.). Shock Compression Condens. Matter - 1997. Amherst, MA: American Institute of Physics; 1998.
[16] Tyburczy JA, Xu X, Ahrens TJ, Epstein S. Earth Planet Sci Lett 2001;192:23.
[17] Hergenrother CW, Barucci MA, Barnouin O, Bierhaus B, Binzel RP, Bottke WF, Chesley S, Clark BC, Clark BE, Cloutis E, D'Aubigny C, Delbo M, Emery J, Gaskell B, Howell E, Keller L, Kelley M, Marshall J, Michel P, Nolan M, Rizk B, Scheeres D, Takir D, Vokrouhlicky, Beshore E, Lauretta DS. The Design Reference Asteroid for the OSIRIS-REx Mission Target ( 101955 ) Bennu. arXiv:1409.4704 [asto-ph.EP]: 2014.
[18] Ghosh A, McSween HY. Meteorit Planet Sci 1999;34:121.
[19] Chesley SR, Farnocchia D, Nolan MC, Vokrouhlicky D, Chodas PW, Milani A, Spoto F, Rozitis B, Benner L a. M, Bottke WF, Busch MW, Emery JP, Howell ES, Lauretta DS, Margot J-L, Taylor P a. Icarus 2014;235:5.
[20] Chesley SR, Farnocchia D, Nolan MC, Vokrouhlicky D, Chodas PW, Milani A, Spoto F, Rozitis B, Benner LAM Bottke WF, Busch MW, Emery JP, Howell ES, Lauretta DS, Margot J, Taylor PA. Icarus 2014;235:5.
[21] Clark BE, Binzel RP, Howell ES, Cloutis EA, Ockert-bell M, Christensen P, Barucci MA, DeMeo F, Lauretta DS, Connolly H, Soderberg A, Hergenrother C, Lim L, Emery J, Mueller M. Icarus 2011;216:462.
[22] Tyburczy JA, Duffy TS, Ahrens TJ, Lange MA. J Geophys Res 1991;96:18011.
[23] Korycansky DG, Plesko CS. Acta Astronaut 2012;73:10.
[24] Herbold E, Lomov I, Miller P, Antoun T. Influence of Morphological and Mechanical Properties on Standoff Mitigation of Potentially Hazardous Asteroids, in:. 44th Lunar Planet. Sci. Conf., vol. 94550. The Woodlands: 2013.
[25] Antoun TH, Lomov I, Glenn L. Development and Application of a Strength and Damage Model for Rock under Dynamic Loading, in:. Elsworth, Tinucci J, Heasley K (Eds.). Proc. 38th U.S. Rock Mech. Symp. Rock Mech. Natl. Interes. A.A. Balkema; 2001.
[26] Lomov I, Rubin MB. Jounal dePhysique IV 2003;110:281.
[27] Hertel E, Bell R. An Improved Material Interface Reconstruction Algorithm for Eulerian Codes. 1992.
[28] Rubin MB, Lomov I. Int J Solids Struct 2003;40:4299.
[29] Vorobiev O. Int J Plast 2008;24:2221.
[30] Berger M, Colella P. J Comput Phys 1989;82:64.
[31] Einav I. Proc R Soc London Ser A 2007;463:3021.
[32] Einav I. J Mech Phys Solids 2007;55:1298.
[33] Rubin MB, Einav I. Int J Eng Sci 2011;49:1151.
[34] Mader CL, Gibbs TR, Hopson JW, Marsh SP, Hoyt MS, Thayer K V, Craig BG, Deal WE. LASL Shock Hugoniot Data. Berkeley, Los Angeles, London: University of California Press; 1980.
[35] Miehe C. 1995:243.
[36] Eckart C. Phys Rev 1948;73:373.
[37] Simo JC. Comput Methods Appl Mech Eng 1992;99:61.
[38] Rubin MB, Attia A. Int J Numer Methods Eng 1996;39:309.
[39] Bar-on E, Rubin MB, Yankelevsky DZ. Int J Solids Struct 2003;40:4519.
[40] Nolan MC, Magri C, Howell ES, Benner LAM, Giorgini JD, Hergenrother CW, Hudson RS, Lauretta DS, Margot JL, Ostro SJ, Scheeres DJ. 2013.
[41] Miller G, Pucket E. J Comput Phys 1996;128:134.
Appendix A. Godunov method implementation in GEODYN
A brief description of the modified Godunov scheme implemented in GEODYN is presented here. The conservation of mass, momentum and energy along with an equation for distortional deformation take the form
^(pF,) + V-(pvF^=pOp (Al)
where Fi denotes the advected history variables. The numerical scheme for a single fluid in a cell is based on the approach described in [41] with the additional capability to treat the full stress tensor associated with distortional deformation of solid materials. The n-dimensional equations (see Equation (A1)) are solved using an operator splitting technique, where the solution is found in each direction at a time,
К,,, |4o,o (-^0 ('4u, (-V. ( Г^)))))))] (A2)
where the operators Sijk,
■wra (A3)
are applied in a Strang-splitting lor second-order accuracy. Hie source term is always applied at the end ol the time step:
^ I'm"..! = • + лад, (A4)
Each directional operator is the update of the cell from one tune step to the next, where the fluxes are computed at the cells edges (2D) or faces (3D). The edge fluxes are computed with an upwind characteristic tracing [41] and Riemann solver in an acoustic approximation. The estimate of the velocity gradient L is computed in the Riemann solver step. The presence of mixed material cells requires treatment of the evolution of each volume fraction of material by the equation:
— + V-(£v) = ^iCV-v, (A5)
д1 ' 'Ka
where/a and Ka are the volume fractions and bulk modulus of each material a. The standard way of treating mixed cells in GEODYN is to have single valued velocity and stress for each material, but allow each material to evolve its constitutive history variables separately. The multimaterial cells update volume fractions according to self-consistent thermodynamics:
1 IK = YjaIKa, Tt=KjjJJKa (A6a)
1 / G = £/„ / G„, 7;,.: = Gj]faTSa / G„ (A6b)
where G and T are the shear modulus and stress tensor for each material and the velocity gradient distribution among
materials is performed in a similar manner: La = LG/Ga. A high-order interface reconstruction algorithm, which preserves linear interfaces during translation. A constraint on volume fractions is satisfied assuming fast equilibrium of partial pressures within the cell. The pressure relaxation algorithm consists of iterative adjustments of volume fractions where the averaged pressure is computed
HZ^i/s^l; (A7)
followed by the updated volume fraction iterate 5fa=pafc(pcrp)/Ka. The parameter p is chosen from physical considerations to allow each volume fraction to fall within a range between 0 and 1 within the cell.
