- Email: [email protected]
A Lagrangian model for debris cloud dynamics simulation
Int. ,L Impact Engng Vol.14, pp.229-240, 1993
0734-743X/93 $6.00+0.00 © 1993 Pergamon Press Ltd
Printed in Great Britain
A LAGRANGIAN MODEL FOR DEBRIS CLOUD DYNAMICS SIMULATION E.P. Fahrenthold Department of Mechanical Engineering, University of Texas Austin, TX 78712
ABSTRACT A new modeling approach has been developed for computer simulation of hypervelocity impacts on multi-plate orbital debris shields. This approach links an Eulerian finite difference code for shield perforation calculations to a Lagrangian finite element code for debris cloud evolution simulations. Mixture theory is used to account for the presence of void space in the debris cloud. INTRODUCTION Most numerical simulations of hypervelocity impact problems have employed Eulerian hydrocodes [e.g. CSQ (Thompson, 1990a) and CTH (McGlaun et al., 1990)], well suited to model perforation and erosion effects. However application of Eulerian codes in the design of multi-plate space debris shields (CourPalais and Crews, 1990 and Christiansen, 1990) suggests that they are not well suited to modeling debris cloud evolution, even in two dimensions. This is apparently due to basic mass dispersion problems associated with the use of Eulerian analysis methods to model very low density debris. Although rezoning, periodic deletion of low density debris, or other techniques may improve the efficiency of Eulerian codes in two dimensional debris modeling problems, such techniques can involve significant requirements for user intervention in the simulation. In addition, three dimensional Eulerian hydrocode analysis of debris cloud evolution appears to be impractical, given current supercomputer capabilities. Hence the present paper describes the development and application of a Lagrangian modeling approach to debris cloud dynamics simulation. Previous work comparing the accuracy of Lagrangian and Eulerian codes in one-dimensional modeling of debris cloud dynamics has indicated that Lagrangian models can provide an accurate description of experimental data (Asay and Trucano, 1990). However the extension of this modeling approach to general debris shielding design problems presents two major difficulties. First, the initial (perforation) portion of the impact problem is very difficult to model accurately with Lagrangian codes. Second, generating a conventional finite element mesh description of the debris cloud is extremely difficult or impractical. As a result, the author is not aware of any previous attempts to extrapolate the one-dimensional work of Asay and Trucano (1990) to the general case. Recognizing the inherent suitability of Eulerian codes for modeling perforation and and the apparent accuracy of Lagrangian codes in modeling debris cloud evolution, the present paper describes a systematic linking of the codes CTH (McGlaun et al., 1990) and DYNA2D (Hallquist, 1987) for the simulation of two dimensional space debris shield impact problems. In the approach described here, the mass and velocity distribution data obtained from a CTH model of initial perforation is post-processed to provide a DYNA2D model of the debris cloud behind the perforated plate. The Lagrangian debris cloud model is then used to simulate the transport of dispersing debris toward the next shield, or the shielded structure. 229
230
E.P.F.kIIRENIHOII)
To simulate the perforation of the next shield in a multi-plate assembly, the DYNA2D model of the expanded debris cloud is post-processed to initiate a new CTH simulation. This process is repeated to progress though the multi-plate shield system and finally model impact on the protected structure. Automated post-processing and mesh generation for the sequence of calculations is performed using newly developed routines written to interface the the DYNA2D and CTH codes. In order to solve the difficult problem of Lagrangian mesh generation for the impact debris cloud, the basic DYNA2D code has been augmented with a mixture theory based (Drumheller and Bedford, 1980) constitutive model of a solid or fluid medium containing voids, including a rate dependent law for void space evolution (Drumheller, 1987). The result is a thermodynamically consistent model of debris cloud evolution for use in space shield design applications. Simulation results show that that the outlined work improves upon existing capabilities for direct hydrocode modeling of such problems. METHODOLOGY The paragraphs which follow outline the methodology used in developing a combined EulerianLagrangian approach to debris cloud dynamics simulation. The Appendix discusses the analysis procedure in detail, while this and later sections present an overview and example simulations. The series of routines written to link the codes CTH and DYNA2D is referred to under the title DCT2D (Debris Cloud Translator 2-Dimensional). The series of routines written to augment the standard DYNA2D code for debris cloud simulation is referred to under the title DCA2D (Debris Cloud Augmentation 2Dimensional). All code development and analysis was performed on a Cray Y-MP/864. The CTH-DYNA interface and DYNA augmentations presented here are currently limited to two dimensional problems, although the general modeling methodology may be implemented in a three-dimensional form. To illustrate the analysis procedure, consider the representative problem of a 0.32 cm diameter aluminum sphere impacting upon a 0.081 cm thick aluminum plate at a velocity of 6.58 km/sec (Fig. la). At a standoff distance of 10.16 cm, this plate represents a ballistic limit Whipple shield for an aluminum wall of 0.127 cm thickness (Cour-Palais and Crews, 1990). CTH simulation of the initial impact in this problem was performed with mesh dimensions Ax = Ay = 0.004 cm. This mesh density exceeds by a factor of 6.25 that used by the CTH code development group in their published study of a canonical debris cloud problem (Trucano and McGlaun, 1990). Transmitting boundary conditions (McGlaun, 1982) were used to accommodate backsplash and debris transport outside the modeled region. The simulation employed the first momentum advection option (CONV=I, Trucano and McGlaun, 1990) in the CTH code. This option conserves momentum in mapping the deformed material mesh to the space-fixed mesh, while discarding any kinetic energy error associated with the remap, and is recommended by the CTH code development group (Trucano and McGlaun, 1990). The ANEOS (Thompson, 1990b) library equation of state for aluminum, with melting, was employed. The calculation required 3,135 CPU seconds.
Y
wall
wall
Y TL x
shield
Y =
O.Ocm
t
X
-
Y :
0.0 cm
plate #4
,
Y : -2.62 cm
plate#3
~
Y : -5.17 cm
plate #2
•
Y = -7.72 c m
plate #1
-
Y : -10.3 c m
Y = -10.3 cm
(~v: 6.58 kmlsec
Fig.la. Whipple shield impact problem
v=6.58k~s~
Fig. lb. Multi-plate shield impact problem
Debris
cloud
dynamics
231
',
-9.0
,
-9.0 -9.2 -9.2
-9.4
-9,4
-9,5
-9.6
-9.8 -9
,8
,~,~-10.0
o.-10.0 >..
-10.2
-10.2
-10.4
-10.4 -10.$
-10.6
-10.8
-10.8
-I1.0 -1 .0
I
I
I
I
-0.$
I
I
-0,2
t
-I1.0
I
0.6
0,2
.0
I .0
-0.6
-0,2
(~)
x
(a) -8.9
Y(cm)
Y(cm) •
* ~ ' . ~ . . ° . . I o
,-.,
-10.1
,;
" - ~
°
;. °.:o
4
-9.9 -0.4
X('cm)
X(cm)
-0.6
+0.4
+0.6
(d)
(c) 1.0
0.0
O
-~..y.g.~/~ ..., at ;.... ;~;~,....:~-
.... ~ ; :
.,*~ ~ "
"!i"., .:.-':" ~: " ,
::.
• .....-,.,
•
.
:
-
~,~
i
•
-~....-¢.,.',
~
-O.S Q
. .
"- .~;':"
"-
S
0.0
:,
;;.. ....
;.':~.::
!~: : • ...
!~
.
O:':::~ "'::''.
--J-. ..... ;; :
•'N~
:~3:,
~:~...~,~ . . , - . . . . -
-6.0
I .0
Co)
-9.5
Y(cm)
0.6
0,2
(~,,)
x
-1
0
.~ ;..
.
-1 .5
. : :-::~
::.:..
:
.'':
; -.
• .
-2.0
-2.5 • . ;
-3.0
. !
. i
i
i
i
i
i
i
X(cm)
1
. . . .
+3.0
-~.0
!
!
|
i
I
.
|
-2.0 -1.5 -1.0 -0.5
x (c,,,,)
(f)
(e) Fig. 2. Whipple shield impact simulation
•
0.0
•
,
0.5
,
F
1.0
,
t
l.S
•
2.0
232
E . P . FAHRENTrtOI.D
Figures 2a and 2b show the results of the CTH simulation, in the form of mesh density plots at impact and at one and one-half microseconds after impact. The plot in Fig. 2b shows a clearly defined debris cloud at one and one-half microseconds after impact, extending approximately over the region -10.0 cm < Y < -9.5 cm and 0.0 cm < X < 0.5 cm, where X is a radial coordinate in this axisymmetric problem, and Y = 0 defines the axial location of the shielded structure (not modeled in this calculation). At this point in time the CTH calculation was terminated and the post-processor CTHED was used to write information on the debris cloud region to a data file for use in constructing a DYNA2D debris model. Next the program DCT2D was run to prepare a DYNA2D input file describing the debris cloud mass, geometry, and velocity distribution, in a Lagrangian finite element form. All elements in the initial mesh containing over ninety-nine percent void space were deleted, resulting in the debris model shown in Fig. 2c. Comparison of the Lagrangian debris cloud model (Fig. 2c) with the debris cloud state at the end of the Eulerian simulation (Fig. 2b) illustrates the accuracy of the model translation procedure. Mass and kinetic energy error associated with the model translation process was approximately one percent. Since the Lagrangian elements contain variable amounts of void space, as determined by the CTH mass distribution data, caution should be exercised in interpreting the finite element geometry plot of Fig. 2c as a direct representation of the debris cloud mass. Given an initial density, velocity, and void fraction distribution obtained from the CTH simulation, the Lagrangian model (DYNA2D augmented by the DCA2D routines) was then integrated to propagate the debris. In this example the DYNA2D model of the debris was composed of 2,181 elements, and required 4,297 CPU seconds to simulate the first microsecond of debris cloud evolution. It is important to note that the Lagrangian calculation must be started with a low user-specified time step, since the default DYNA2D calculation for the initial time step will not account for void space effects. At two and one-half microseconds after impact, one microsecond after ending the Eulerian simulation, the debris cloud has evolved to the form shown in Fig. 2d (note the change in scale as compared to Fig. 2c). The radius of the debris cloud has increased significantly, and the leading edge of the debris cloud has translated along the impact centerline. The shell of the debris cloud has thinned as it expands, while the spreading of debris particles within the shell illustrates the presence of velocity variations across the debris cloud. The analysis reflects (qualitatively) results observed in impact experiments. The Lagrangian simulation of debris cloud evolution shown in Figures 2c and 2d assumed that adjacent finite elements composing the debris cloud were interconnected, i.e. that the debris cloud deforms as a cohesive body, retaining its mechanical strength after impact. An alternative modeling option incorporated in the DCT2D routines takes each finite element to represent a discrete debris cloud fragment, i.e. assumes that the debris cloud lacks any cohesive strength after impact. This is often a more realistic description of hypervelocity impact effects on Whipple shield structures. Consider again the representative Whipple shield impact problem depicted in Figures 2a and 2b. Starting with the initial Lagrangian debris cloud model (Fig. 2c), and neglecting any cohesive coupling of the motion of the deforming elements, the augmented version of DYNA2D was used to propagate the impact debris towards a wall plate located at axial position Y = 0. At 17.5 microseconds after impact the leading edge of the debris cloud has reached the wall plate, with the expanded debris cloud taking the form shown in Fig. 2e. This simulation required 13,320 CPU seconds. Comparing Figures 2c and 2e, the debris cloud has expanded radially by a factor of eight and axially by a factor of thirteen. The predicted angle of expansion of the cone of debris is approximately thirty-six degrees. Qualitatively the debris cloud shown in Fig. 2e is consistent with experimental observations of Whipple shield impacts. To complete the simulation and model debris impact on the wall plate, the mass, position, and velocity data from the DYNA2D simulation at 17.5 microseconds after impact (Fig. 2e) was input to DCT2D to generate a CTH input file. CTH was then used to simulate impact of the debris on the wall structure. The wall plate was modeled as aluminum with a yield strength of 0.966 x 109 dynes/cm2. Figure 2f shows the wall at 19.5 microseconds after initial impact on the bumper shield, or 2.0 microseconds after initial impact of the debris on the wall. This CTH simulation required 3,035 CPU seconds. Overall the simulation predicts perforation of the wall plate for this nominally ballistic limit shield configuration. Comparison of the preceding simulation results with the corresponding experiment (Cour-Palais and Crews, 1990) suggests that the analysis results are conservative, i.e. that the lethality of the debris cloud has been overestimated. However the results indicate that the proposed analysis procedure provides a
Debris cloud dynamics
233
stable and realistic transformation of the Eulerian impact modeling results to a Lagrangian form well adapted for debris cloud evolution calculations. General use of the simulation approach just discussed in debris shield design applications requires further evaluation of the methods employed. A later section describes a series of calculations conducted to simulate a multi-plate shield impact experiment. MIXTURE THEORY MODEL The preceding discussion noted the development of a new constitutive augmentation of the DYNA2D code, describing a solid or fluid continuum with voids, in order to make practical the use of this code in multi-dimensional debris modeling applications. The following paragraphs outline briefly the formulation and implementation of a mixture theory based model for debris cloud simulation. This model has been coded (as a user-defined, history dependent equation of state) in a Cray (UNICOS) version of DYNA2D. The model is currently implemented in an isothermal, hydrodynamic form, although extension to include non-isothermal and deviatoric stress effects is possible. For a discussion of mixture theory concepts, see Drumheller and Bedford (1980). Consider a solid or fluid continuum with voids, described by a bulk density p, a true density ~t, a void fraction ¢, a bulk pressure P, and a true pressure Ps" The parameters ~/and Ps represent the thermodynamic state of the material of interest, while the parameters P and p are corresponding average values for a Lagrangian f'mite element of fixed mass M and variable volume V. The preceding parameters are related by P = (1-¢) Ps ; P = (1-¢) ~/ ; p = M/V ; Ps = Ps (T)
(la,b,c,d)
where the expression (ld) represents the true isothermal equation of state for the material of interest. The history dependent element equation of state which is required for the Lagrangian finite element code has the functional form P = P (T, ¢)
(2)
In order to define the element state an additional equation is needed for ¢. Consistent with Drumheller's (1987) work on hypervelocity impact of mixtures, a rate equation will be employed here. The rate equation may be formulated with the aid of the general isothermal entropy inequality (Malvern, 1969) - + + (I/p) tr[TD] _>0
(3)
where W is the Helmholtz free energy density, T is the Cauchy stress tensor, D is the rate of deformation tensor, and "tr" is the trace operator. Considering only volumetric deformation, let T =-PI ; tr(D)=-p/p
(4a,b)
where I is the identity tensor. The entropy inequality then reduces to
- ~i'
+
(P/p) [l~IP] -> 0
(5a)
Using equations (1) and the conventional free energy density assumption tlS = W(T), the inequality (5a) reduces to [(1-¢)P/p 2 - 3W/~TI ~- (PT/p 2) ~ > 0
(5b)
If a rate law is assumed for ¢, as previously discussed, this implies constitutive relations of a general functional form
3~13~t = (1-¢)P/p 2 ; ~) = 0 [PTI92] ; PTIp 2 = Ps/p ; (1-(I))P/p2 = Ps]T2
(6a,b & 7a,b)
where the brackets in equation (6b) indicate functional dependence. Numerical implementation in DYNA2D of the mixture model just described assumed the following particular functional form for equation (6a)
234
E . P . F~HREN[HOLD
(8a,b)
Ps = Cl IX+ c 2 max(Ix,0)+ c 3 IX2 + (c4 + c5ix + c6 IX2) E ; IX= (T/'/o) - 1
where the parameters ci (i=1,2,...6) are constants, E is the internal energy, and To is the reference density. Consistent with the derived entropy inequality, the void fraction evolution equation was taken as = - [~/(To/p)2] ( P/[(1-~p)p] }
(8c)
where a is a constant which determines the rate of change of the void fraction in response to imposed pressure or deformation, The case (x -- 0 represents a constant void fraction, whereas for the "shifting equilibrium" case (Bowen, 1982) represented by (x ---) oo the void fraction adjusts instantaneously to changes in the bulk density. The latter case is representative of classical porous media models. A total of fourteen FORTRAN routines were added to DYNA2D or modified from their basic DYNA2D form in order to incorporate the mixture model just discussed into that code. MODEL TRANSLATION PROCEDURE The Eulerian-to-Lagrangian and Lagrangian-to-Eulerian model conversions used here require proper translation of mass, void fraction, and velocity distributions. Certain aspects of the model translation process warrant elaboration. In the Eulerian-to-Lagrangian (CTH to DYNA2D) model translation process, each Eulerian cell is mapped to a Lagrangian finite element of the same size and geometry. Of course this exact geometric correspondence holds only at the start of the Lagrangian calculation, since the finite elements move and deform during the simulation, while the Eulerian ceils are space fixed. The true density and void fraction for each cell at the end of the Eulerian simulation are assigned directly to a corresponding Lagrangian finite element. Translation of the velocity data is complicated slightly by the fact that the finite difference scheme used in CTH provides velocity data at the midpoints of the cell boundaries, while the finite element scheme in DYNA2D requires initial conditions data at nodal points corresponding to the Eulerian cell comers. Hence the following formulas were used to assign the initial Lagrangian nodal velocities . . . . (i) (j) Vx Od) and Vy Od) at the point located by the position coordinates x and y • (i,j) vx Vy
(ij-1) = [w x
(i,j)
where the v x
(i j)
(i,j-1) vx
(i-l,j)
= [ Wy
(i j)
and Vy
(i,j) + wx
(i-l,j)
Vy
(i j) vx
(i,j)
+ Wy
(i j - l ) ]/[w x
(id)
Vy
(ij) +w x
(i-l,j)
] / [ Wy
(i,j)
+ Wy
]
(9a)
]
(9b)
are the Eulerian velocity data, the superscripts (i), (j), and (i,j) denote logical
(index) coordinates in the CTH mesh (McGlaan et al., 1990), and wx
(i,j)
= (1/2) [ m
(i- 1 j)
+m
(i,j)
(ij)
] ; Wy
= (1/2) [ m
(id- 1)
with m (i'j) denoting the mass in cell (i j). The weighting factors w x
+m
(i j)
(id)
]
(9c,d)
(id)
and Wy
are appropriate for
the uniform mesh used here, and are required in order to account for spatial variations in the true density and void fraction. In the Lagrangian-to-Eulerian (DYNA2D to CTH) model translation process, one or two mass insertion packages for each finite element are written to the CTH input file. Triangular elements require one insertion package while quadrilateral elements require two, one for each of two triangular parts into which the quadrilateral is bisected. (The DCIED routines provide an option for the use of either quadrilateral or triangular finite elements.) The material in each insertion package is assigned a velocity equal to the average of the associated nodal velocities, and a density equal to the true density of the solid material in
Debris cloud dynamics
235
the associated finite element. The latter quantity is calculated using the element void fraction, which is maintained as a history variable in the Lagrangian simulation. Each triangular insert package is located in space as follows. First the coordinates of the centroid of the triangular insertion package (Xc,Y c) are calculated from the associated nodal coordinates at the end of the Lagrangian simulation(Xi,Yi). Then the coordinates of the corners of the triangular insertion package (xi,Yi) are calculated using the formulas: xi=cXi+(1.c)Xc
; Yi=cYi+(1-c)Yc
; c = ( 1 - 0 p ) l / 2 ; ie {1,2,3}
(10a,b,c)
The mapping of equations (10) represents a uniform dilatation, as seen in the element cross section, and was adopted in order to conserve both mass and true density in the model translation process. A subdivision of quadrilateral elements into two triangles is employed in the Lagrangian-to-Eulerian model translation process, so that in the case of radial symmetry the centroidal (cross section) coordinates of the element subdivisions and their associated insertion packages are identical. MULTI-PLATE SPACE DEBRIS SHIELDING This section describes a representative multi-plate shield impact modeling problem, simulated using the coordinated Eulerian-Lagrangian approach previously outlined. Parameters of the simulations are listed in Table 1. The problem involves normal impact of a 0.32 cm diameter sphere on a series of four bumper plates of thickness 0.0102 cm, followed by a wall of thickness 0.079 cm (Fig. lb). The plate-to-plate and plate-to-wall spacing was 2.54 cm. The corresponding experiment is described by Cour-Palais and Crews (1990). The CTH simulations were performed using the same mesh density, boundary conditions, and other modeling options previously discussed. The DYNA2D simulations were performed using models composed of up to 6,944 finite elements. The debris cloud was modeled as an isotropic, elastic-plastic hydrodynamic material (aluminum) with voids, and the following material properties: shear modulus = 0.250 g/(cm-~tsec2), bulk modulus = 6.52 x 10"lg/(cm-gsee2), yield strength = 3.45 x 10"3g/(cm-gsec2), plastic hardening modulus = 6.67 x 10"2g/(cm-gsec2), and values of t~ over the range 10-1 < o~< 10+3 (l.tsec/cm2) as indicated in Table 1. Variations in the parameter t~ were considered in order to investigate the results on the simulations. Analyses to date indicate that the results are insensitive to the choice of values for tx, although numerical stability requirements appear to place an upper limit on allowable values for tx in the explicit DYNA2D code. The following paragraphs briefly outline the series of nine Eulerian and Lagrangian simulations required to model perforation of all four shields and impact on the protected structure. The simulations are referred to by the codes listed in the fast column of Table 1. The total required CPU time was 18.8 hours, for a simulation time of 20.8 microseconds. The first CTH simulation ("plate #1" in Table 1) modeled perforation of the the first shield plate (Fig. 3a). Figure 3b shows a mass density plot at 1,5 microseconds after impact. At that time the CTH simulation was terminated and a DYNA2D model of the debris cloud was generated (Fig. 3c). This model ("debris #1" in Table 1) was used to simulate motion of the debris cloud towards the second plate, requiting a simulation time of 2.4 microseconds. The state of the debris cloud at 3.9 microseconds after impact, i.e. at the end of this Lagrangian simulation, is shown in Fig. 3d. At this point the results of the Lagrangian calculation were used to generate a new Eulerian model of impact on the second plate (Fig. 3e). The results of the second CTH simulation ("plate #2" in Table 1) are shown in Fig. 3f, a mass density plot at 5.5 microseconds after impact on the first shield. Comparison of the debris cloud models at the end of the first Lagrangian simulation (Fig. 3d) and the start of the second Eulerian simulation (Fig. 3e) illustrates the model translation process. Note that in the interest of reducing CPU time requirements, some of the widely dispersed debris present at the end of the first DYNA2D simulation (Fig. 3d) is neglected in the Lagrangian-to-Eulerian remap. CTH simulation of the second shield impact was followed by a second Lagrangian simulation of debris cloud evolution towards the third shield ("debris #2" in Table 1). As indicated in Table 1, this sequence of calculations was repeated to proceed through the entire multi-
236
~,. P ,
F M t P , E N I HOLD
-9.0
-9.0
-9.2
-9.2
-9.,4
-9.4
-9.6
-9.1~
-9.8
-S.O
~,.~- 1 0 . 0 >.-
~'- 10.0 >..
-10.2
-10•2
®
-10.4 -10.6
-10.4 -10,6
-10,8
-10,8
-11.0 - I .0
1
I
-0.$
I
l
I
i
-0.2
i
I
0.2
0.6
I
.0
-11,0 -1
t
!
i
i
-0.6
.0
i
f
i
i
0.2
-0.2
x (~)
x
(a)
i
0.G
.0
(~)
(b)
-9.2
-7.8
Y(cm)
• ~ -.,.. ~ . ~ . ~ : , . . , , . . . . ,
Y(cm)
. . . . . . .•. .. ;.:'.." ..
;,-,t . , , ' . . . ..,~,_,%:~. - " .... ..-.::.; ..
.',.. ,,:& ..
.° ,...
.
•
. t/
...
"t~ .
.. "
-9.1
-9.9
-0.3
X(cm)
+0.3
-0.7
X(cm)
+0.7
(d)
(c)
-8.6
-8.8
-6.8
-8.8
-7.0
-7,0
-7.2
-7.2
-7.4
-7.4 :
"~-7,6
T
~-
:.- - 7 . a
-7.B
-7.8
.b
(,
"
-8.0
-8 .0
• - ,~.
-8.2
:
,~ .('/.,
:./-'
--
-6,2
-8.4
-8 .a
-8.6
-8.6
-8.8
-8,1~ -I
.0
-O.B
-0.2
0.2 x
0.6
.0
-1.0
i
i
i
-0.6
i
i
-0.2
(~) (f)
(e)
Fig. 3. Multi-plate shield impact simulation
i
0.Z
l
~
0.6
i
t
1.0
Debris cloud dynamics
237
plate structure. To reduce CPU time requirements, only the central core of the debris cloud was propagated towards the wall plate. That is debris subject to wide radial dispersion was dropped from the simulation during the Eulerian-to-Lagrangian or Lagrangian-to-Eulerian rezones. No Lagrangian-toLagrangian rezones were required to complete the analysis. Final simulation of debris cloud impact on the wall structure is depicted in Fig. 4, which shows a mass density plot at 20.8 microseconds after initial impact on the f'n'st shield. As indicated by the plot in Fig. 4, the predicted result is a hole in the protected structure with an approximate diameter of 0.16 cm. This compares favorably with the experimental result of a torn wall plate, with damage dimensions approximately 0.2 cm x 0.5 cm (Cour-Palais and Crews, 1990). Simulation title
plate #1 debris #1 plate #2 debris #2 plate #3 debris #3 plate #4 debris #4 wall plate
Simulation Starttime End time CPU time alpha type (micro- (micro- (seconds) ( m i c r o E=Eulerian seconds) seconds) Cray secondsper L=Lagrangian Y-MP/864 em squared)
E L E L E L E L E
0.0 1.5 3.9 5.9 8.3 10.1 12.1 14.8 17.8
1.5 3.9 5.9 8.3 10.1 12.1 14.8 17.8 20.8
2,913 7,481 6,249 8,168 5,860 10,360 8,676 11,850 6,135
m a x i m u m Numberof element Eulerianor void Lagrangian fraction zones (%)
0.1
99
1000.0
98
1.0
99
1.0
99
250 x 500 2,423 250 x 500 6,944 250 x 750 6,510 250 x 750 4,462 250 x 750
Table 1. Simulation parameters: multi-plate shield impact problem CONCLUSIONS A number of analytical models have been developed for projectile impact on Whipple shields (see e.g. Grady and Passman, 1990). These models typically adopt a number of very basic assumptions regarding the debris cloud shape, mass distribution, velocity distribution, and other properties. These assumptions have been motivated in part by difficulties experienced in modeling debris cloud evolution with Eulerian hydrocodes. The modeling methodology presented here avoids both the major assumptions of analytical models and mass dispersion problems which may be encountered with purely Eulerian simulations. The results presented here suggest two conclusions regarding the use of a combined Eulerian-Lagrangian hydrocode modeling approach to space debris shield design : (1) Mixture theory based models for solid and fluid materials with voids provide a suitable means for extending the use of Lagrangian hydrocodes to multidimensional debris cloud modeling problems. (2) Lagrangian debris cloud models can offer acceptable CPU time requirements for direct computer simulation of multi-plate impact experiments, at least in two dimensions, while requiring minimal use~ intervention in the simulation. Future work can profitably focus on two areas. First, additional simulation work is needed, to further critique the modeling methodology used here against a range of hypervelocity impact experiments. Of particular importance is the effect of variations in material properties and constitutive equations on model predictions at various impact velocities. Second, additional software development work is needed to provide the capability for three dimensional simulation of oblique impact effects on proposed space debris shield designs. Extension of the modeling approach presented here to three dimensions is direct, and well motivated by the extremely large computer time requirements of three dimensional Eulerian hydrocode simulations.
238
E.P.
FAHRENTHOLD
dccth
,
0.50
,
,
,
,
i
,
,
i
,
.
©
,
0.25
0.00
-.i..
%8:
Eulerlan simulation of shield impact and perforation
rscth
"
..:::~::i..
-0.25
!~" .. :;( ~;:-, 4-'~.
-0.50
-0.75
"7 ,~%- I .
-
dcout
O0
-I
. 2S
-1
.50
-1
.75
,
,o
dcout -2.
O0
-2.25 i
-2.50 -I
.5
-1
.0
i -0.5
m
i 0.5
0.0
t
i 1.0
1.5
dcllnkout
x (cm)
dcg2dinp~
Fig. 4. Wall impact for the multi-plate shield impact simulation
Eulerian to Lagranglan rezone
stats.out dcg2dout rvl2ddat statsout + dcbase + dcg2dout dcdyna
Fig. A-1. Eulerian simulation and rezone procedure
rvl2d dat dcdyna (rnodIMed f I l LLNLcode)
1
Lagrangian simulationOf debris cloud evolution
dcrezdat dcrezdat rvl2ddat dcrezInp .-.I~ dcrezf I dcg2d.out
Lagrangian to Eulerian
1
dcrezout
statsout dcg2dout
rezone
dccthl + dcrezout + dccth2 1
dccth
Fragment model wlth quadrllateralelements statsp.out dcgp2dout rvl2ddat . i i , ~
' 1
dccth
©
Pre-processing for a Lagrangiento Lagranglan rezone
rscth Fragment model wlth triangularelements
statsq.out
dcgq2dout rvlq2ddat
Fig. A-2. Fragmented debris cloud model generadon procedure
1
dcout
_
Fig. A-3. Lagrangian simulation and rezone procedure
Debris cloud dynamics
239
ACKNOWLEDGEMENTS This work was supported under the NASA Regional Universities Grant Program. The assistance of Eric L. Christiansen (NASA Technical Officer) and Jeanne Lee Crews of Johnson Space Center is gratefully acknowledged. Additional funding was provided by Cray Research, Inc. under the Cray University Research and Development Grant Program. Computer time was provided by the University of Texas System Center for High Performance Computing. REFERENCES Asay, J.R. and T.C. Trucano (1990). Studies of Density Distributions in One-Dimensional ShockInduced Debris Clouds. Int. J. of Impact Engineering, 10, 35-50. Bowen, R.M. (1982). Compressible Porous Media Models by Use of the Theory of Mixtures. Int. J. of Engineering Science, 20, 697-735. Christiansen, E.L. (1990). Advanced Meteoroid and Debris Shielding Concepts. AIAA 90-1336, AIAA/NASA/DOD Orbital Debris Conference, Baltimore, April 16-19. Cour-Palais, B.G. and J.L. Crews (1990). A Multi-Shock Concept for Spacecraft Shielding. Int. J. of Impact Engineering, 10, 135-146. Drumheller, D.S. (1987). Hypervelocity Impact of Mixtures. Int. J. of Impact Engineering, 5, 261-268. Drumheller, D.S. and A. Bedford (1980). A Thermomechanical Theory for Reacting Immiscible Mixtures. Archives for Rational Mechanics and Analysis, 73, 257-284. Grady, D.E. and S.L. Passman (1990). Stability and Fragmentation of Ejecta in Hypervelocity Impact. Int. J. of Impact Engineering, 10, 197-212. Hallquist, J.O. (1987). User's Manual for DYNA2D. UCID 18756 (Rev. 3), Lawrence Livermore National Laboratory, Livermore, California. Malvern, L.E. (1969). Introduction to the Mechanics of a Continuous Medium. Prentice-Hall, Englewood Cliffs, New Jersey. McGlaun. J.M. (1982). Improvements in CSQ :A Transmitting Boundary Condition. SAND82-1248, Sandia National Laboratories, Albuquerque, New Mexico. McGlaun, J.M., S.L. Thompson and M.G. Elrick (1990). CTH: A Three Dimensional Shock Wave Physics Code. Int. J. of Impact Engineering, 10, 351-360. Thompson, S.L. (1990a). CSQ III: An Eulerian Finite Difference Program for Two-Dimensional Material Response: User's Manual. SAND87-2763, Sandia National Laboratories, Albuquerque, New Mexico. Thompson, S.L. (1990b). ANEOS - Analytic Equations of State for Shock Physics Codes: Input Manual. SAND89-295 I, Sandia National Laboratories, Albuquerque, New Mexico. Trucano, T.G. and J.M. McGlaun (1990). Hypervelocity Impact Simulations Using CTH: Case Studies. Int. J. of Impact Engineering, 10, 601-614. APPENDIX This appendix provides a detailed description of the structure and use of the software written to link CTH and DYNA2D. The general analysis procedure is represented by the flow charts shown in Figures A-1 through A-3. As presently formulated, the procedure models two-dimensional (axisymmetric) impact of like materials. Mass, kinetic energy, void fraction, and density are conserved in Eulerian-to-Lagrangian or Lagrangian-to-Eulerian translations of debris cloud data. Thermal energy is discarded. Generalization of this analysis procedure is possible, given additional development work. A typical impact simulation proceeds as follows. A standard Eulerian model of the projectile impact on the first debris shield is formulated using CTI-I (Fig. A-I). Geometry, material properties, initial velocities, and other input data are specified in the input file dccth. The simulation is halted after perforation of the plate and initial formation of the debris cloud, but before the bulk debris cloud density is reduced to a level at which the Eulerian mesh density becomes inadequate. The output f'de (rscth) from the CTH simulation is then input to the post-processor CTHED. The standard CTH post-processor CTHED is then used to write an output data file which describes the state (mass, velocity, etc.) of all Eulerian cells in a rectangular region of space behind the shield. (The user selects the region of space containing the debris which must be propagated to the next shield, by viewing the CTH simulation
240
E . P . FAHRENI'ItOtD
results.) The mass, velocity, and pressure distribution in the debris cloud is contained in the postprocessor output file dcout. Next (Fig. A-I) the routines dclinkfand dcg2dfare used to generate a Lagrangian model of the debris cloud (dcdyna), given the Eulerian simulation results (dcout) and the geometry, material, and control data contained in the user input files dcg2d.inp and dcbase. Specifically, the output data file dcout, which is quite extensive and is written in a specific CTHED format, is read and screened by dclink.f, which writes an output data file dclink.out that contains only the position, mass, velocity, void fraction, and other data needed to construct the Lagrangian finite element mesh with appropriate initial conditions. Then the program dcg2df is invoked, which reads the data file dclink.out and performs the following three functions: (1) Defines the nodes, elements, element connectivities, and initial velocities for the Lagrangian finite element mesh used to model the debris cloud. Initial velocities are assigned to the finite element nodes based on interpolation of the CTH cell velocity data. The user may choose to delete all element interconnections, essentially simulating the debris cloud as a collection of noninteracting deforming elements ("fragments") with voids (Fig. A-2). The use of this option, well suited to many hypervelocity problems, is illustrated in the example simulations discussed in the text. (2) Screens the preliminary finite element mesh to identify and delete all elements with void space percentages above a user-specified level. This necessitates a somewhat complex renumbering of nodes and elements to satisfy DYNA2D input format requirements. This step is important to avoid the inclusion of very low mass, high void space elements which may greatly slow the computation while representing very little debris mass. The preceding results are written to dcg2d.out in a format suitable for direct inclusion in a DYNA2D input file. (3) Defines the initial void space associated with each finite element in the Lagrangian model, writing the results to rvl2d.dat which becomes an auxiliary input data file for the DYNA2D code. Combined with the new constitutive modeling routines previously discussed, this procedure provides a mass and kinetic energy consistent interface between the CTH and DYNA2D codes without requiring that each debris particle be explicitly modeled by the Lagrangian mesh. Such a requirement would in general be impractical from both a model generation and CPU time consumption point of view. The output file dcg2d.out is appended to a short user-prepared control file dcdata (written in a DYNA2D specified forma0 which provides standard information on material properties, time step size, etc. needed to perform the Lagrangian simulation. The resulting file dcdyna is input to an augmented version of the Lagrangian code DYNA2D (Fig. A-3), which incorporates a mixture equation of state for a continuum with voids, as well as a void fraction evolution equation. The auxiliary input file rvl2d.dat created by dcg2df is required to properly initialize the DYNA2D calculation. A Lagrangian simulation is then used to propagate the debris to the next shield. An auxiliary output file dcrez.dat is created by the modified version of DYNA2D for use in initializing the next Eulerian shield impact calculation. The results of the Lagrangian simulation of the debris cloud evolution (dcrez.dat) are input to the routine dcrezf(Fig. A-3), to develop an Eulerian description of the debris cloud impacting the next shield. The routine dcrez.f also requires as input certain geometry, mesh connectivity, and mass distribution data contained in the files dcrez.inp, rvl2d.dat, and dcg2d.out, the latter two having been previously generated by dcg2d.f. The output from dcrez.f, contained in the file dcrez.out, is combined directly with standard CTH specified geometry, material, and control data (dccth.1 and dccth.2) describing the next shield (or wall) impact calculation to be performed. The resulting CTH input file (dccth) is used to initiate a repeat cycle of calculations. In some cases, one or more rezones of the Lagrangian mesh may be required to propagate the debris between two adjacent shields. In this case (Fig. A-3), the file dccth is created as previously discussed, but without a new shield model, and then input to the CTH pre-processor CTHGEN. The resulting output file (rscth) is then processed as previously described in order to create a rezoned Lagrangian model. This Lagrangian rezone procedure is well suited to mixture theory based debris models, although a special (direc0 rezone routine could be written for this purpose. Note that the standard DYNA2D rezoner cannot be used in this case, due to the presence of void space in the Lagrangian finite elements. In any case, conventional Lagrangian rezone procedures are ill-suited to this application, due to the extremely complex geometry of the debris clouds. In summary, the programs DCT2D and DCA2D provide a highly automated coupling of CTH and DYNA2D for use in debris cloud modeling problems.
0734-743X/93 $6.00+0.00 © 1993 Pergamon Press Ltd
Printed in Great Britain
A LAGRANGIAN MODEL FOR DEBRIS CLOUD DYNAMICS SIMULATION E.P. Fahrenthold Department of Mechanical Engineering, University of Texas Austin, TX 78712
ABSTRACT A new modeling approach has been developed for computer simulation of hypervelocity impacts on multi-plate orbital debris shields. This approach links an Eulerian finite difference code for shield perforation calculations to a Lagrangian finite element code for debris cloud evolution simulations. Mixture theory is used to account for the presence of void space in the debris cloud. INTRODUCTION Most numerical simulations of hypervelocity impact problems have employed Eulerian hydrocodes [e.g. CSQ (Thompson, 1990a) and CTH (McGlaun et al., 1990)], well suited to model perforation and erosion effects. However application of Eulerian codes in the design of multi-plate space debris shields (CourPalais and Crews, 1990 and Christiansen, 1990) suggests that they are not well suited to modeling debris cloud evolution, even in two dimensions. This is apparently due to basic mass dispersion problems associated with the use of Eulerian analysis methods to model very low density debris. Although rezoning, periodic deletion of low density debris, or other techniques may improve the efficiency of Eulerian codes in two dimensional debris modeling problems, such techniques can involve significant requirements for user intervention in the simulation. In addition, three dimensional Eulerian hydrocode analysis of debris cloud evolution appears to be impractical, given current supercomputer capabilities. Hence the present paper describes the development and application of a Lagrangian modeling approach to debris cloud dynamics simulation. Previous work comparing the accuracy of Lagrangian and Eulerian codes in one-dimensional modeling of debris cloud dynamics has indicated that Lagrangian models can provide an accurate description of experimental data (Asay and Trucano, 1990). However the extension of this modeling approach to general debris shielding design problems presents two major difficulties. First, the initial (perforation) portion of the impact problem is very difficult to model accurately with Lagrangian codes. Second, generating a conventional finite element mesh description of the debris cloud is extremely difficult or impractical. As a result, the author is not aware of any previous attempts to extrapolate the one-dimensional work of Asay and Trucano (1990) to the general case. Recognizing the inherent suitability of Eulerian codes for modeling perforation and and the apparent accuracy of Lagrangian codes in modeling debris cloud evolution, the present paper describes a systematic linking of the codes CTH (McGlaun et al., 1990) and DYNA2D (Hallquist, 1987) for the simulation of two dimensional space debris shield impact problems. In the approach described here, the mass and velocity distribution data obtained from a CTH model of initial perforation is post-processed to provide a DYNA2D model of the debris cloud behind the perforated plate. The Lagrangian debris cloud model is then used to simulate the transport of dispersing debris toward the next shield, or the shielded structure. 229
230
E.P.F.kIIRENIHOII)
To simulate the perforation of the next shield in a multi-plate assembly, the DYNA2D model of the expanded debris cloud is post-processed to initiate a new CTH simulation. This process is repeated to progress though the multi-plate shield system and finally model impact on the protected structure. Automated post-processing and mesh generation for the sequence of calculations is performed using newly developed routines written to interface the the DYNA2D and CTH codes. In order to solve the difficult problem of Lagrangian mesh generation for the impact debris cloud, the basic DYNA2D code has been augmented with a mixture theory based (Drumheller and Bedford, 1980) constitutive model of a solid or fluid medium containing voids, including a rate dependent law for void space evolution (Drumheller, 1987). The result is a thermodynamically consistent model of debris cloud evolution for use in space shield design applications. Simulation results show that that the outlined work improves upon existing capabilities for direct hydrocode modeling of such problems. METHODOLOGY The paragraphs which follow outline the methodology used in developing a combined EulerianLagrangian approach to debris cloud dynamics simulation. The Appendix discusses the analysis procedure in detail, while this and later sections present an overview and example simulations. The series of routines written to link the codes CTH and DYNA2D is referred to under the title DCT2D (Debris Cloud Translator 2-Dimensional). The series of routines written to augment the standard DYNA2D code for debris cloud simulation is referred to under the title DCA2D (Debris Cloud Augmentation 2Dimensional). All code development and analysis was performed on a Cray Y-MP/864. The CTH-DYNA interface and DYNA augmentations presented here are currently limited to two dimensional problems, although the general modeling methodology may be implemented in a three-dimensional form. To illustrate the analysis procedure, consider the representative problem of a 0.32 cm diameter aluminum sphere impacting upon a 0.081 cm thick aluminum plate at a velocity of 6.58 km/sec (Fig. la). At a standoff distance of 10.16 cm, this plate represents a ballistic limit Whipple shield for an aluminum wall of 0.127 cm thickness (Cour-Palais and Crews, 1990). CTH simulation of the initial impact in this problem was performed with mesh dimensions Ax = Ay = 0.004 cm. This mesh density exceeds by a factor of 6.25 that used by the CTH code development group in their published study of a canonical debris cloud problem (Trucano and McGlaun, 1990). Transmitting boundary conditions (McGlaun, 1982) were used to accommodate backsplash and debris transport outside the modeled region. The simulation employed the first momentum advection option (CONV=I, Trucano and McGlaun, 1990) in the CTH code. This option conserves momentum in mapping the deformed material mesh to the space-fixed mesh, while discarding any kinetic energy error associated with the remap, and is recommended by the CTH code development group (Trucano and McGlaun, 1990). The ANEOS (Thompson, 1990b) library equation of state for aluminum, with melting, was employed. The calculation required 3,135 CPU seconds.
Y
wall
wall
Y TL x
shield
Y =
O.Ocm
t
X
-
Y :
0.0 cm
plate #4
,
Y : -2.62 cm
plate#3
~
Y : -5.17 cm
plate #2
•
Y = -7.72 c m
plate #1
-
Y : -10.3 c m
Y = -10.3 cm
(~v: 6.58 kmlsec
Fig.la. Whipple shield impact problem
v=6.58k~s~
Fig. lb. Multi-plate shield impact problem
Debris
cloud
dynamics
231
',
-9.0
,
-9.0 -9.2 -9.2
-9.4
-9,4
-9,5
-9.6
-9.8 -9
,8
,~,~-10.0
o.-10.0 >..
-10.2
-10.2
-10.4
-10.4 -10.$
-10.6
-10.8
-10.8
-I1.0 -1 .0
I
I
I
I
-0.$
I
I
-0,2
t
-I1.0
I
0.6
0,2
.0
I .0
-0.6
-0,2
(~)
x
(a) -8.9
Y(cm)
Y(cm) •
* ~ ' . ~ . . ° . . I o
,-.,
-10.1
,;
" - ~
°
;. °.:o
4
-9.9 -0.4
X('cm)
X(cm)
-0.6
+0.4
+0.6
(d)
(c) 1.0
0.0
O
-~..y.g.~/~ ..., at ;.... ;~;~,....:~-
.... ~ ; :
.,*~ ~ "
"!i"., .:.-':" ~: " ,
::.
• .....-,.,
•
.
:
-
~,~
i
•
-~....-¢.,.',
~
-O.S Q
. .
"- .~;':"
"-
S
0.0
:,
;;.. ....
;.':~.::
!~: : • ...
!~
.
O:':::~ "'::''.
--J-. ..... ;; :
•'N~
:~3:,
~:~...~,~ . . , - . . . . -
-6.0
I .0
Co)
-9.5
Y(cm)
0.6
0,2
(~,,)
x
-1
0
.~ ;..
.
-1 .5
. : :-::~
::.:..
:
.'':
; -.
• .
-2.0
-2.5 • . ;
-3.0
. !
. i
i
i
i
i
i
i
X(cm)
1
. . . .
+3.0
-~.0
!
!
|
i
I
.
|
-2.0 -1.5 -1.0 -0.5
x (c,,,,)
(f)
(e) Fig. 2. Whipple shield impact simulation
•
0.0
•
,
0.5
,
F
1.0
,
t
l.S
•
2.0
232
E . P . FAHRENTrtOI.D
Figures 2a and 2b show the results of the CTH simulation, in the form of mesh density plots at impact and at one and one-half microseconds after impact. The plot in Fig. 2b shows a clearly defined debris cloud at one and one-half microseconds after impact, extending approximately over the region -10.0 cm < Y < -9.5 cm and 0.0 cm < X < 0.5 cm, where X is a radial coordinate in this axisymmetric problem, and Y = 0 defines the axial location of the shielded structure (not modeled in this calculation). At this point in time the CTH calculation was terminated and the post-processor CTHED was used to write information on the debris cloud region to a data file for use in constructing a DYNA2D debris model. Next the program DCT2D was run to prepare a DYNA2D input file describing the debris cloud mass, geometry, and velocity distribution, in a Lagrangian finite element form. All elements in the initial mesh containing over ninety-nine percent void space were deleted, resulting in the debris model shown in Fig. 2c. Comparison of the Lagrangian debris cloud model (Fig. 2c) with the debris cloud state at the end of the Eulerian simulation (Fig. 2b) illustrates the accuracy of the model translation procedure. Mass and kinetic energy error associated with the model translation process was approximately one percent. Since the Lagrangian elements contain variable amounts of void space, as determined by the CTH mass distribution data, caution should be exercised in interpreting the finite element geometry plot of Fig. 2c as a direct representation of the debris cloud mass. Given an initial density, velocity, and void fraction distribution obtained from the CTH simulation, the Lagrangian model (DYNA2D augmented by the DCA2D routines) was then integrated to propagate the debris. In this example the DYNA2D model of the debris was composed of 2,181 elements, and required 4,297 CPU seconds to simulate the first microsecond of debris cloud evolution. It is important to note that the Lagrangian calculation must be started with a low user-specified time step, since the default DYNA2D calculation for the initial time step will not account for void space effects. At two and one-half microseconds after impact, one microsecond after ending the Eulerian simulation, the debris cloud has evolved to the form shown in Fig. 2d (note the change in scale as compared to Fig. 2c). The radius of the debris cloud has increased significantly, and the leading edge of the debris cloud has translated along the impact centerline. The shell of the debris cloud has thinned as it expands, while the spreading of debris particles within the shell illustrates the presence of velocity variations across the debris cloud. The analysis reflects (qualitatively) results observed in impact experiments. The Lagrangian simulation of debris cloud evolution shown in Figures 2c and 2d assumed that adjacent finite elements composing the debris cloud were interconnected, i.e. that the debris cloud deforms as a cohesive body, retaining its mechanical strength after impact. An alternative modeling option incorporated in the DCT2D routines takes each finite element to represent a discrete debris cloud fragment, i.e. assumes that the debris cloud lacks any cohesive strength after impact. This is often a more realistic description of hypervelocity impact effects on Whipple shield structures. Consider again the representative Whipple shield impact problem depicted in Figures 2a and 2b. Starting with the initial Lagrangian debris cloud model (Fig. 2c), and neglecting any cohesive coupling of the motion of the deforming elements, the augmented version of DYNA2D was used to propagate the impact debris towards a wall plate located at axial position Y = 0. At 17.5 microseconds after impact the leading edge of the debris cloud has reached the wall plate, with the expanded debris cloud taking the form shown in Fig. 2e. This simulation required 13,320 CPU seconds. Comparing Figures 2c and 2e, the debris cloud has expanded radially by a factor of eight and axially by a factor of thirteen. The predicted angle of expansion of the cone of debris is approximately thirty-six degrees. Qualitatively the debris cloud shown in Fig. 2e is consistent with experimental observations of Whipple shield impacts. To complete the simulation and model debris impact on the wall plate, the mass, position, and velocity data from the DYNA2D simulation at 17.5 microseconds after impact (Fig. 2e) was input to DCT2D to generate a CTH input file. CTH was then used to simulate impact of the debris on the wall structure. The wall plate was modeled as aluminum with a yield strength of 0.966 x 109 dynes/cm2. Figure 2f shows the wall at 19.5 microseconds after initial impact on the bumper shield, or 2.0 microseconds after initial impact of the debris on the wall. This CTH simulation required 3,035 CPU seconds. Overall the simulation predicts perforation of the wall plate for this nominally ballistic limit shield configuration. Comparison of the preceding simulation results with the corresponding experiment (Cour-Palais and Crews, 1990) suggests that the analysis results are conservative, i.e. that the lethality of the debris cloud has been overestimated. However the results indicate that the proposed analysis procedure provides a
Debris cloud dynamics
233
stable and realistic transformation of the Eulerian impact modeling results to a Lagrangian form well adapted for debris cloud evolution calculations. General use of the simulation approach just discussed in debris shield design applications requires further evaluation of the methods employed. A later section describes a series of calculations conducted to simulate a multi-plate shield impact experiment. MIXTURE THEORY MODEL The preceding discussion noted the development of a new constitutive augmentation of the DYNA2D code, describing a solid or fluid continuum with voids, in order to make practical the use of this code in multi-dimensional debris modeling applications. The following paragraphs outline briefly the formulation and implementation of a mixture theory based model for debris cloud simulation. This model has been coded (as a user-defined, history dependent equation of state) in a Cray (UNICOS) version of DYNA2D. The model is currently implemented in an isothermal, hydrodynamic form, although extension to include non-isothermal and deviatoric stress effects is possible. For a discussion of mixture theory concepts, see Drumheller and Bedford (1980). Consider a solid or fluid continuum with voids, described by a bulk density p, a true density ~t, a void fraction ¢, a bulk pressure P, and a true pressure Ps" The parameters ~/and Ps represent the thermodynamic state of the material of interest, while the parameters P and p are corresponding average values for a Lagrangian f'mite element of fixed mass M and variable volume V. The preceding parameters are related by P = (1-¢) Ps ; P = (1-¢) ~/ ; p = M/V ; Ps = Ps (T)
(la,b,c,d)
where the expression (ld) represents the true isothermal equation of state for the material of interest. The history dependent element equation of state which is required for the Lagrangian finite element code has the functional form P = P (T, ¢)
(2)
In order to define the element state an additional equation is needed for ¢. Consistent with Drumheller's (1987) work on hypervelocity impact of mixtures, a rate equation will be employed here. The rate equation may be formulated with the aid of the general isothermal entropy inequality (Malvern, 1969) - + + (I/p) tr[TD] _>0
(3)
where W is the Helmholtz free energy density, T is the Cauchy stress tensor, D is the rate of deformation tensor, and "tr" is the trace operator. Considering only volumetric deformation, let T =-PI ; tr(D)=-p/p
(4a,b)
where I is the identity tensor. The entropy inequality then reduces to
- ~i'
+
(P/p) [l~IP] -> 0
(5a)
Using equations (1) and the conventional free energy density assumption tlS = W(T), the inequality (5a) reduces to [(1-¢)P/p 2 - 3W/~TI ~- (PT/p 2) ~ > 0
(5b)
If a rate law is assumed for ¢, as previously discussed, this implies constitutive relations of a general functional form
3~13~t = (1-¢)P/p 2 ; ~) = 0 [PTI92] ; PTIp 2 = Ps/p ; (1-(I))P/p2 = Ps]T2
(6a,b & 7a,b)
where the brackets in equation (6b) indicate functional dependence. Numerical implementation in DYNA2D of the mixture model just described assumed the following particular functional form for equation (6a)
234
E . P . F~HREN[HOLD
(8a,b)
Ps = Cl IX+ c 2 max(Ix,0)+ c 3 IX2 + (c4 + c5ix + c6 IX2) E ; IX= (T/'/o) - 1
where the parameters ci (i=1,2,...6) are constants, E is the internal energy, and To is the reference density. Consistent with the derived entropy inequality, the void fraction evolution equation was taken as = - [~/(To/p)2] ( P/[(1-~p)p] }
(8c)
where a is a constant which determines the rate of change of the void fraction in response to imposed pressure or deformation, The case (x -- 0 represents a constant void fraction, whereas for the "shifting equilibrium" case (Bowen, 1982) represented by (x ---) oo the void fraction adjusts instantaneously to changes in the bulk density. The latter case is representative of classical porous media models. A total of fourteen FORTRAN routines were added to DYNA2D or modified from their basic DYNA2D form in order to incorporate the mixture model just discussed into that code. MODEL TRANSLATION PROCEDURE The Eulerian-to-Lagrangian and Lagrangian-to-Eulerian model conversions used here require proper translation of mass, void fraction, and velocity distributions. Certain aspects of the model translation process warrant elaboration. In the Eulerian-to-Lagrangian (CTH to DYNA2D) model translation process, each Eulerian cell is mapped to a Lagrangian finite element of the same size and geometry. Of course this exact geometric correspondence holds only at the start of the Lagrangian calculation, since the finite elements move and deform during the simulation, while the Eulerian ceils are space fixed. The true density and void fraction for each cell at the end of the Eulerian simulation are assigned directly to a corresponding Lagrangian finite element. Translation of the velocity data is complicated slightly by the fact that the finite difference scheme used in CTH provides velocity data at the midpoints of the cell boundaries, while the finite element scheme in DYNA2D requires initial conditions data at nodal points corresponding to the Eulerian cell comers. Hence the following formulas were used to assign the initial Lagrangian nodal velocities . . . . (i) (j) Vx Od) and Vy Od) at the point located by the position coordinates x and y • (i,j) vx Vy
(ij-1) = [w x
(i,j)
where the v x
(i j)
(i,j-1) vx
(i-l,j)
= [ Wy
(i j)
and Vy
(i,j) + wx
(i-l,j)
Vy
(i j) vx
(i,j)
+ Wy
(i j - l ) ]/[w x
(id)
Vy
(ij) +w x
(i-l,j)
] / [ Wy
(i,j)
+ Wy
]
(9a)
]
(9b)
are the Eulerian velocity data, the superscripts (i), (j), and (i,j) denote logical
(index) coordinates in the CTH mesh (McGlaan et al., 1990), and wx
(i,j)
= (1/2) [ m
(i- 1 j)
+m
(i,j)
(ij)
] ; Wy
= (1/2) [ m
(id- 1)
with m (i'j) denoting the mass in cell (i j). The weighting factors w x
+m
(i j)
(id)
]
(9c,d)
(id)
and Wy
are appropriate for
the uniform mesh used here, and are required in order to account for spatial variations in the true density and void fraction. In the Lagrangian-to-Eulerian (DYNA2D to CTH) model translation process, one or two mass insertion packages for each finite element are written to the CTH input file. Triangular elements require one insertion package while quadrilateral elements require two, one for each of two triangular parts into which the quadrilateral is bisected. (The DCIED routines provide an option for the use of either quadrilateral or triangular finite elements.) The material in each insertion package is assigned a velocity equal to the average of the associated nodal velocities, and a density equal to the true density of the solid material in
Debris cloud dynamics
235
the associated finite element. The latter quantity is calculated using the element void fraction, which is maintained as a history variable in the Lagrangian simulation. Each triangular insert package is located in space as follows. First the coordinates of the centroid of the triangular insertion package (Xc,Y c) are calculated from the associated nodal coordinates at the end of the Lagrangian simulation(Xi,Yi). Then the coordinates of the corners of the triangular insertion package (xi,Yi) are calculated using the formulas: xi=cXi+(1.c)Xc
; Yi=cYi+(1-c)Yc
; c = ( 1 - 0 p ) l / 2 ; ie {1,2,3}
(10a,b,c)
The mapping of equations (10) represents a uniform dilatation, as seen in the element cross section, and was adopted in order to conserve both mass and true density in the model translation process. A subdivision of quadrilateral elements into two triangles is employed in the Lagrangian-to-Eulerian model translation process, so that in the case of radial symmetry the centroidal (cross section) coordinates of the element subdivisions and their associated insertion packages are identical. MULTI-PLATE SPACE DEBRIS SHIELDING This section describes a representative multi-plate shield impact modeling problem, simulated using the coordinated Eulerian-Lagrangian approach previously outlined. Parameters of the simulations are listed in Table 1. The problem involves normal impact of a 0.32 cm diameter sphere on a series of four bumper plates of thickness 0.0102 cm, followed by a wall of thickness 0.079 cm (Fig. lb). The plate-to-plate and plate-to-wall spacing was 2.54 cm. The corresponding experiment is described by Cour-Palais and Crews (1990). The CTH simulations were performed using the same mesh density, boundary conditions, and other modeling options previously discussed. The DYNA2D simulations were performed using models composed of up to 6,944 finite elements. The debris cloud was modeled as an isotropic, elastic-plastic hydrodynamic material (aluminum) with voids, and the following material properties: shear modulus = 0.250 g/(cm-~tsec2), bulk modulus = 6.52 x 10"lg/(cm-gsee2), yield strength = 3.45 x 10"3g/(cm-gsec2), plastic hardening modulus = 6.67 x 10"2g/(cm-gsec2), and values of t~ over the range 10-1 < o~< 10+3 (l.tsec/cm2) as indicated in Table 1. Variations in the parameter t~ were considered in order to investigate the results on the simulations. Analyses to date indicate that the results are insensitive to the choice of values for tx, although numerical stability requirements appear to place an upper limit on allowable values for tx in the explicit DYNA2D code. The following paragraphs briefly outline the series of nine Eulerian and Lagrangian simulations required to model perforation of all four shields and impact on the protected structure. The simulations are referred to by the codes listed in the fast column of Table 1. The total required CPU time was 18.8 hours, for a simulation time of 20.8 microseconds. The first CTH simulation ("plate #1" in Table 1) modeled perforation of the the first shield plate (Fig. 3a). Figure 3b shows a mass density plot at 1,5 microseconds after impact. At that time the CTH simulation was terminated and a DYNA2D model of the debris cloud was generated (Fig. 3c). This model ("debris #1" in Table 1) was used to simulate motion of the debris cloud towards the second plate, requiting a simulation time of 2.4 microseconds. The state of the debris cloud at 3.9 microseconds after impact, i.e. at the end of this Lagrangian simulation, is shown in Fig. 3d. At this point the results of the Lagrangian calculation were used to generate a new Eulerian model of impact on the second plate (Fig. 3e). The results of the second CTH simulation ("plate #2" in Table 1) are shown in Fig. 3f, a mass density plot at 5.5 microseconds after impact on the first shield. Comparison of the debris cloud models at the end of the first Lagrangian simulation (Fig. 3d) and the start of the second Eulerian simulation (Fig. 3e) illustrates the model translation process. Note that in the interest of reducing CPU time requirements, some of the widely dispersed debris present at the end of the first DYNA2D simulation (Fig. 3d) is neglected in the Lagrangian-to-Eulerian remap. CTH simulation of the second shield impact was followed by a second Lagrangian simulation of debris cloud evolution towards the third shield ("debris #2" in Table 1). As indicated in Table 1, this sequence of calculations was repeated to proceed through the entire multi-
236
~,. P ,
F M t P , E N I HOLD
-9.0
-9.0
-9.2
-9.2
-9.,4
-9.4
-9.6
-9.1~
-9.8
-S.O
~,.~- 1 0 . 0 >.-
~'- 10.0 >..
-10.2
-10•2
®
-10.4 -10.6
-10.4 -10,6
-10,8
-10,8
-11.0 - I .0
1
I
-0.$
I
l
I
i
-0.2
i
I
0.2
0.6
I
.0
-11,0 -1
t
!
i
i
-0.6
.0
i
f
i
i
0.2
-0.2
x (~)
x
(a)
i
0.G
.0
(~)
(b)
-9.2
-7.8
Y(cm)
• ~ -.,.. ~ . ~ . ~ : , . . , , . . . . ,
Y(cm)
. . . . . . .•. .. ;.:'.." ..
;,-,t . , , ' . . . ..,~,_,%:~. - " .... ..-.::.; ..
.',.. ,,:& ..
.° ,...
.
•
. t/
...
"t~ .
.. "
-9.1
-9.9
-0.3
X(cm)
+0.3
-0.7
X(cm)
+0.7
(d)
(c)
-8.6
-8.8
-6.8
-8.8
-7.0
-7,0
-7.2
-7.2
-7.4
-7.4 :
"~-7,6
T
~-
:.- - 7 . a
-7.B
-7.8
.b
(,
"
-8.0
-8 .0
• - ,~.
-8.2
:
,~ .('/.,
:./-'
--
-6,2
-8.4
-8 .a
-8.6
-8.6
-8.8
-8,1~ -I
.0
-O.B
-0.2
0.2 x
0.6
.0
-1.0
i
i
i
-0.6
i
i
-0.2
(~) (f)
(e)
Fig. 3. Multi-plate shield impact simulation
i
0.Z
l
~
0.6
i
t
1.0
Debris cloud dynamics
237
plate structure. To reduce CPU time requirements, only the central core of the debris cloud was propagated towards the wall plate. That is debris subject to wide radial dispersion was dropped from the simulation during the Eulerian-to-Lagrangian or Lagrangian-to-Eulerian rezones. No Lagrangian-toLagrangian rezones were required to complete the analysis. Final simulation of debris cloud impact on the wall structure is depicted in Fig. 4, which shows a mass density plot at 20.8 microseconds after initial impact on the f'n'st shield. As indicated by the plot in Fig. 4, the predicted result is a hole in the protected structure with an approximate diameter of 0.16 cm. This compares favorably with the experimental result of a torn wall plate, with damage dimensions approximately 0.2 cm x 0.5 cm (Cour-Palais and Crews, 1990). Simulation title
plate #1 debris #1 plate #2 debris #2 plate #3 debris #3 plate #4 debris #4 wall plate
Simulation Starttime End time CPU time alpha type (micro- (micro- (seconds) ( m i c r o E=Eulerian seconds) seconds) Cray secondsper L=Lagrangian Y-MP/864 em squared)
E L E L E L E L E
0.0 1.5 3.9 5.9 8.3 10.1 12.1 14.8 17.8
1.5 3.9 5.9 8.3 10.1 12.1 14.8 17.8 20.8
2,913 7,481 6,249 8,168 5,860 10,360 8,676 11,850 6,135
m a x i m u m Numberof element Eulerianor void Lagrangian fraction zones (%)
0.1
99
1000.0
98
1.0
99
1.0
99
250 x 500 2,423 250 x 500 6,944 250 x 750 6,510 250 x 750 4,462 250 x 750
Table 1. Simulation parameters: multi-plate shield impact problem CONCLUSIONS A number of analytical models have been developed for projectile impact on Whipple shields (see e.g. Grady and Passman, 1990). These models typically adopt a number of very basic assumptions regarding the debris cloud shape, mass distribution, velocity distribution, and other properties. These assumptions have been motivated in part by difficulties experienced in modeling debris cloud evolution with Eulerian hydrocodes. The modeling methodology presented here avoids both the major assumptions of analytical models and mass dispersion problems which may be encountered with purely Eulerian simulations. The results presented here suggest two conclusions regarding the use of a combined Eulerian-Lagrangian hydrocode modeling approach to space debris shield design : (1) Mixture theory based models for solid and fluid materials with voids provide a suitable means for extending the use of Lagrangian hydrocodes to multidimensional debris cloud modeling problems. (2) Lagrangian debris cloud models can offer acceptable CPU time requirements for direct computer simulation of multi-plate impact experiments, at least in two dimensions, while requiring minimal use~ intervention in the simulation. Future work can profitably focus on two areas. First, additional simulation work is needed, to further critique the modeling methodology used here against a range of hypervelocity impact experiments. Of particular importance is the effect of variations in material properties and constitutive equations on model predictions at various impact velocities. Second, additional software development work is needed to provide the capability for three dimensional simulation of oblique impact effects on proposed space debris shield designs. Extension of the modeling approach presented here to three dimensions is direct, and well motivated by the extremely large computer time requirements of three dimensional Eulerian hydrocode simulations.
238
E.P.
FAHRENTHOLD
dccth
,
0.50
,
,
,
,
i
,
,
i
,
.
©
,
0.25
0.00
-.i..
%8:
Eulerlan simulation of shield impact and perforation
rscth
"
..:::~::i..
-0.25
!~" .. :;( ~;:-, 4-'~.
-0.50
-0.75
"7 ,~%- I .
-
dcout
O0
-I
. 2S
-1
.50
-1
.75
,
,o
dcout -2.
O0
-2.25 i
-2.50 -I
.5
-1
.0
i -0.5
m
i 0.5
0.0
t
i 1.0
1.5
dcllnkout
x (cm)
dcg2dinp~
Fig. 4. Wall impact for the multi-plate shield impact simulation
Eulerian to Lagranglan rezone
stats.out dcg2dout rvl2ddat statsout + dcbase + dcg2dout dcdyna
Fig. A-1. Eulerian simulation and rezone procedure
rvl2d dat dcdyna (rnodIMed f I l LLNLcode)
1
Lagrangian simulationOf debris cloud evolution
dcrezdat dcrezdat rvl2ddat dcrezInp .-.I~ dcrezf I dcg2d.out
Lagrangian to Eulerian
1
dcrezout
statsout dcg2dout
rezone
dccthl + dcrezout + dccth2 1
dccth
Fragment model wlth quadrllateralelements statsp.out dcgp2dout rvl2ddat . i i , ~
' 1
dccth
©
Pre-processing for a Lagrangiento Lagranglan rezone
rscth Fragment model wlth triangularelements
statsq.out
dcgq2dout rvlq2ddat
Fig. A-2. Fragmented debris cloud model generadon procedure
1
dcout
_
Fig. A-3. Lagrangian simulation and rezone procedure
Debris cloud dynamics
239
ACKNOWLEDGEMENTS This work was supported under the NASA Regional Universities Grant Program. The assistance of Eric L. Christiansen (NASA Technical Officer) and Jeanne Lee Crews of Johnson Space Center is gratefully acknowledged. Additional funding was provided by Cray Research, Inc. under the Cray University Research and Development Grant Program. Computer time was provided by the University of Texas System Center for High Performance Computing. REFERENCES Asay, J.R. and T.C. Trucano (1990). Studies of Density Distributions in One-Dimensional ShockInduced Debris Clouds. Int. J. of Impact Engineering, 10, 35-50. Bowen, R.M. (1982). Compressible Porous Media Models by Use of the Theory of Mixtures. Int. J. of Engineering Science, 20, 697-735. Christiansen, E.L. (1990). Advanced Meteoroid and Debris Shielding Concepts. AIAA 90-1336, AIAA/NASA/DOD Orbital Debris Conference, Baltimore, April 16-19. Cour-Palais, B.G. and J.L. Crews (1990). A Multi-Shock Concept for Spacecraft Shielding. Int. J. of Impact Engineering, 10, 135-146. Drumheller, D.S. (1987). Hypervelocity Impact of Mixtures. Int. J. of Impact Engineering, 5, 261-268. Drumheller, D.S. and A. Bedford (1980). A Thermomechanical Theory for Reacting Immiscible Mixtures. Archives for Rational Mechanics and Analysis, 73, 257-284. Grady, D.E. and S.L. Passman (1990). Stability and Fragmentation of Ejecta in Hypervelocity Impact. Int. J. of Impact Engineering, 10, 197-212. Hallquist, J.O. (1987). User's Manual for DYNA2D. UCID 18756 (Rev. 3), Lawrence Livermore National Laboratory, Livermore, California. Malvern, L.E. (1969). Introduction to the Mechanics of a Continuous Medium. Prentice-Hall, Englewood Cliffs, New Jersey. McGlaun. J.M. (1982). Improvements in CSQ :A Transmitting Boundary Condition. SAND82-1248, Sandia National Laboratories, Albuquerque, New Mexico. McGlaun, J.M., S.L. Thompson and M.G. Elrick (1990). CTH: A Three Dimensional Shock Wave Physics Code. Int. J. of Impact Engineering, 10, 351-360. Thompson, S.L. (1990a). CSQ III: An Eulerian Finite Difference Program for Two-Dimensional Material Response: User's Manual. SAND87-2763, Sandia National Laboratories, Albuquerque, New Mexico. Thompson, S.L. (1990b). ANEOS - Analytic Equations of State for Shock Physics Codes: Input Manual. SAND89-295 I, Sandia National Laboratories, Albuquerque, New Mexico. Trucano, T.G. and J.M. McGlaun (1990). Hypervelocity Impact Simulations Using CTH: Case Studies. Int. J. of Impact Engineering, 10, 601-614. APPENDIX This appendix provides a detailed description of the structure and use of the software written to link CTH and DYNA2D. The general analysis procedure is represented by the flow charts shown in Figures A-1 through A-3. As presently formulated, the procedure models two-dimensional (axisymmetric) impact of like materials. Mass, kinetic energy, void fraction, and density are conserved in Eulerian-to-Lagrangian or Lagrangian-to-Eulerian translations of debris cloud data. Thermal energy is discarded. Generalization of this analysis procedure is possible, given additional development work. A typical impact simulation proceeds as follows. A standard Eulerian model of the projectile impact on the first debris shield is formulated using CTI-I (Fig. A-I). Geometry, material properties, initial velocities, and other input data are specified in the input file dccth. The simulation is halted after perforation of the plate and initial formation of the debris cloud, but before the bulk debris cloud density is reduced to a level at which the Eulerian mesh density becomes inadequate. The output f'de (rscth) from the CTH simulation is then input to the post-processor CTHED. The standard CTH post-processor CTHED is then used to write an output data file which describes the state (mass, velocity, etc.) of all Eulerian cells in a rectangular region of space behind the shield. (The user selects the region of space containing the debris which must be propagated to the next shield, by viewing the CTH simulation
240
E . P . FAHRENI'ItOtD
results.) The mass, velocity, and pressure distribution in the debris cloud is contained in the postprocessor output file dcout. Next (Fig. A-I) the routines dclinkfand dcg2dfare used to generate a Lagrangian model of the debris cloud (dcdyna), given the Eulerian simulation results (dcout) and the geometry, material, and control data contained in the user input files dcg2d.inp and dcbase. Specifically, the output data file dcout, which is quite extensive and is written in a specific CTHED format, is read and screened by dclink.f, which writes an output data file dclink.out that contains only the position, mass, velocity, void fraction, and other data needed to construct the Lagrangian finite element mesh with appropriate initial conditions. Then the program dcg2df is invoked, which reads the data file dclink.out and performs the following three functions: (1) Defines the nodes, elements, element connectivities, and initial velocities for the Lagrangian finite element mesh used to model the debris cloud. Initial velocities are assigned to the finite element nodes based on interpolation of the CTH cell velocity data. The user may choose to delete all element interconnections, essentially simulating the debris cloud as a collection of noninteracting deforming elements ("fragments") with voids (Fig. A-2). The use of this option, well suited to many hypervelocity problems, is illustrated in the example simulations discussed in the text. (2) Screens the preliminary finite element mesh to identify and delete all elements with void space percentages above a user-specified level. This necessitates a somewhat complex renumbering of nodes and elements to satisfy DYNA2D input format requirements. This step is important to avoid the inclusion of very low mass, high void space elements which may greatly slow the computation while representing very little debris mass. The preceding results are written to dcg2d.out in a format suitable for direct inclusion in a DYNA2D input file. (3) Defines the initial void space associated with each finite element in the Lagrangian model, writing the results to rvl2d.dat which becomes an auxiliary input data file for the DYNA2D code. Combined with the new constitutive modeling routines previously discussed, this procedure provides a mass and kinetic energy consistent interface between the CTH and DYNA2D codes without requiring that each debris particle be explicitly modeled by the Lagrangian mesh. Such a requirement would in general be impractical from both a model generation and CPU time consumption point of view. The output file dcg2d.out is appended to a short user-prepared control file dcdata (written in a DYNA2D specified forma0 which provides standard information on material properties, time step size, etc. needed to perform the Lagrangian simulation. The resulting file dcdyna is input to an augmented version of the Lagrangian code DYNA2D (Fig. A-3), which incorporates a mixture equation of state for a continuum with voids, as well as a void fraction evolution equation. The auxiliary input file rvl2d.dat created by dcg2df is required to properly initialize the DYNA2D calculation. A Lagrangian simulation is then used to propagate the debris to the next shield. An auxiliary output file dcrez.dat is created by the modified version of DYNA2D for use in initializing the next Eulerian shield impact calculation. The results of the Lagrangian simulation of the debris cloud evolution (dcrez.dat) are input to the routine dcrezf(Fig. A-3), to develop an Eulerian description of the debris cloud impacting the next shield. The routine dcrez.f also requires as input certain geometry, mesh connectivity, and mass distribution data contained in the files dcrez.inp, rvl2d.dat, and dcg2d.out, the latter two having been previously generated by dcg2d.f. The output from dcrez.f, contained in the file dcrez.out, is combined directly with standard CTH specified geometry, material, and control data (dccth.1 and dccth.2) describing the next shield (or wall) impact calculation to be performed. The resulting CTH input file (dccth) is used to initiate a repeat cycle of calculations. In some cases, one or more rezones of the Lagrangian mesh may be required to propagate the debris between two adjacent shields. In this case (Fig. A-3), the file dccth is created as previously discussed, but without a new shield model, and then input to the CTH pre-processor CTHGEN. The resulting output file (rscth) is then processed as previously described in order to create a rezoned Lagrangian model. This Lagrangian rezone procedure is well suited to mixture theory based debris models, although a special (direc0 rezone routine could be written for this purpose. Note that the standard DYNA2D rezoner cannot be used in this case, due to the presence of void space in the Lagrangian finite elements. In any case, conventional Lagrangian rezone procedures are ill-suited to this application, due to the extremely complex geometry of the debris clouds. In summary, the programs DCT2D and DCA2D provide a highly automated coupling of CTH and DYNA2D for use in debris cloud modeling problems.