- Email: [email protected]
Penetration of a rigid projectile into a finite thickness elastic–plastic target — comparison between theory and numerical computations
International Journal of Impact Engineering 25 (2001) 265}290
Penetration of a rigid projectile into a "nite thickness elastic}plastic target * comparison between theory and numerical computations G. Yossifon, M.B. Rubin*, A.L. Yarin Faculty of Mechanical Engineering, Technion, Israel Institute of Technology, 32000 Haifa, Israel Received 28 January 2000; received in revised form 19 June 2000
Abstract Recently, two di!erent analytical solutions have been developed for penetration of a rigid projectile into an elastic}plastic target of "nite thickness. One solution satis"es the balance of linear momentum pointwise in target region, but it satis"es the free-surface boundary conditions only in an integrated sense. The second solution satis"es the balance of linear momentum only along a "nite number of instantaneous streamlines, but it satis"es the boundary conditions exactly at each intersection of the streamlines with the boundary. The "rst solution is valid only for normal penetration whereas the second solution has been generalized for oblique penetration. The main objective of the present paper is to compare the predictions of these two theoretical approaches with a numerical solution obtained using the hydrocode Autodyn2D. By modifying the boundary condition used in the "rst method it is possible to obtain a reasonably accurate description of the penetration process (including the exit stage) with the computational time reduced from a few hours for Autodyn2D to only a few minutes for the analytical solution. 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction A number of theoretical models have been developed for normal penetration of a rigid projectile into a "nite thickness elastic}plastic target, e.g. [1,2] and references therein. One theoretical approach is to use the kinematics of an inviscid #uid in order to generate an approximate velocity "eld in the target. The origins of this approach date back to [3}5]. Speci"cally, the velocity "eld is taken to be an irrotational "eld that is derivable from a velocity potential. In [2,6] the method of * Corresponding author. 0734-743X/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 7 3 4 - 7 4 3 X ( 0 0 ) 0 0 0 4 0 - 3
266
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Nomenclature A B C D F G H I L M P R R R R R R * R U ; Y c e,e, h e,e, h f g n r rH p s t t v w, w G x z z
e e n
virtual mass of the target coe$cient of the form drag plastic drag symmetric part of the velocity gradient (the rate of deformation tensor) resultant axial drag force acting on the projectile function in Eq. (21) target thickness the unit tensor projectile's length projectile's mass penetration depth radius of the projectile projectile's radius at the intersection with the rear surface of the target projectile's radius at the separation point projectile's radius at the intersection with the front surface of the target projectile's radius at the point on its surface which satis"es condition (32) projectile's radius at its tail asymptotic value of the projectile's radius velocity of the projectile initial velocity of the projectile constant yield strength in uniaxial stress speed of sound base vectors of a cylindrical polar coordinate system base vectors of the cylindrical polar coordinate system associated with the projectile function in Eq. (13) function in Eq. (16) unit vector normal to the elastic}plastic boundary or to the projectile's surface radial coordinate radial position of the point of intersection of the elastic}plastic boundary with the rear surface pressure shock parameter time traction (stress) vector absolute velocity of a material point in the target weighting functions axial coordinate of the projectile's tip axial coordinate intersection of the elastic}plastic boundary with the axis r"0.
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
267
Greek letters C e f , f h k m m , m o r, r s
t
Gruneisen parameter linear strain tensor dimensionless coordinates in Eqs. (17) and (22) angular coordinate shear modulus axial coordinate measured relative to the projectile dimensionless coordinates in Eqs. (17) and (22) mass density Cauchy stress tensor and its deviatoric part unit tangent vector to the projectile's surface velocity potential stream function
Subscripts A P n r z q
indicates indicates indicates indicates indicates indicates
that that that that that that
the the the the the the
parameter parameter parameter parameter parameter parameter
corresponds corresponds corresponds corresponds corresponds corresponds
to to to to to to
point A (see Fig. 3) the projectile the normal component the radial component the axial component the tangential component
Superscripts (e) (p) 䢇
-
indicates that the parameter corresponds to the elastic "eld indicates that the parameter corresponds to the plastic "eld denotes time di!erentiation denotes dimensionless parameters
singularities has been used to develop a velocity "eld that exactly satis"es the conditions of incompressibility in the target region and impenetrability of the projectile for any reasonable projectile shape. In [2] this method was applied to normal penetration and in [6] it was applied to oblique penetration. More speci"cally, when the projectile is taken to be an ovoid of Rankine (Fig. 1) and the impact is normal, the associated velocity "eld is quite simple and is known to be characterized by a combination of a single source and uniform #ow. The solution procedure assumes that the target material is incompressible. Then, the target region is divided into an elastic region ahead of the projectile, where the strains remain in"nitesimal, and a rigid-plastic region near the projectile's tip, where the strains can be very large. Furthermore, the rigid-plastic response is taken to be rate-independent with constant yield strength.
268
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 1. Actual scaled projectile's tip shaped as an ovoid of Rankine.
For normal penetration of an ovoid of Rankine [2], the velocity "eld is relatively simple and the divergence of the deviatoric stress becomes the gradient of a scalar potential. Consequently, the balance of linear momentum can be solved exactly in both the elastic and rigid-plastic regions to determine an expression for the pressure "eld and the stresses. Then, the e!ects of the free front and rear surfaces of the target, and the continuity of the surface tractions at the elastic/rigid-plastic boundary are satis"ed in an approximate integrated sense. This is because the assumed velocity "eld is not exact for the problem of penetration into a "nite target. Speci"cally, the continuity condition at the elastic}plastic boundary is satis"ed only at the intersection with the symmetry axis. The condition that the rear surface remains traction free is satis"ed either at the symmetry axis when the rear surface remains elastic, or in an integrated sense when it becomes plastic. The decelerating force applied to the projectile by the target is calculated by analytically integrating the surface tractions over the portion of the projectile's surface that is in contact with the target. Then, the resulting rigid-body equation of motion of the projectile is solved numerically to determine the penetration depth and velocity as functions of time. The theoretical predictions were compared with numerous experiments [2,6,7], and the agreement was found to be good even though the formulation uses no adjustable parameters. Only the mechanical properties of the materials, and projectile and target geometrical dimensions are used. For oblique penetration [6], the method of singularities yields a more complicated velocity "eld since an additional distribution of doublets is required to account for the projectile's rotation during the penetration process. For this case, the divergence of the deviatoric stress is no longer equal to the gradient of a scalar. Consequently, the balance of linear momentum cannot be solved exactly for the assumed velocity "eld. However, the balance of linear momentum can be solved in an approximate sense by integrating pointwise along a "nite set of instantaneous streamlines. Using this procedure it is possible to satisfy the continuity conditions on stress at the elastic}plastic boundary and the stress-free boundary conditions at portions of the target's front and rear surfaces, both at the intersections of each streamline with these surfaces. The method of [6] can also be applied to the simpler case of normal impact, as is discussed below. Here, method I denotes the procedure which satis"es the balance of linear momentum exactly, but the continuity conditions and boundary conditions only approximately [2]. Also, method II (the method of [6] applied to normal impact) denotes the procedure which satis"es the balance of linear momentum and the continuity and boundary conditions only along instantaneous stream-
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
269
lines. Consequently, the balance of linear momentum is satis"ed more accurately in method I than in method II. Whereas, the continuity and boundary conditions are satis"ed more accurately in method II than in method I. Thus, it is not clear for a particular problem of normal impact, which of these methods will predict more accurate values of the speci"c quantities of interest. The main objective of the present paper is to examine the accuracy of these two methods for normal penetration of a rigid projectile into rigid-plastic target. This is accomplished by comparing the predictions of the two approximate theoretical approaches with numerical simulations using the hydrocode Autodyn2D [8]. The Autodyn2D calculations include the compressibility of the target material and they determine the velocity "eld which is a numerical approximation of the exact solution of the balance of linear momentum and the continuity and boundary conditions. Consequently, it is possible to compare the theoretical and numerical predictions of such global quantities as the drag force and such local quantities as the stress acting on the projectile, the location of the elastic}plastic boundary and the velocity "eld. Also, a modi"cation of the stress-free boundary condition at the target's rear surface is proposed following the work in [9] to improve the theoretical predictions of method I. An outline of the paper is as follows: Section 2 reviews the basic equations of the two theoretical methods; Section 3 describes the boundary and continuity conditions for method I; Section 4 describes the boundary and continuity conditions for method II; Section 5 discusses additional conditions for both methods I and II; Section 6 presents the main results of example problems which compare the theoretical predictions of methods I and II with the results obtained in calculations using Autodyn2D; and Section 7 presents the main conclusions. 2. Basic equations 2.1. The coordinate system and the governing equations Consider normal penetration of a rigid projectile into an incompressible elastic}plastic target of thickness H, mass density o, shear modulus k and constant yield stress Y (in uniaxial stress). The projectile is assumed to have the shape of an ovoid of Rankine with radius R , length L, and mass M. Given these geometrical and material parameters, as well as the projectile's initial position and velocity, the theoretical model predicts the resulting penetration process. In particular, it predicts the instantaneous position and velocity of the projectile, as well as the "nal depth of penetration if the target defeats the projectile, or the residual velocity if the projectile perforates the target. Due to the axial symmetry of the problem, it is convenient to consider a cylindrical polar coordinate system with coordinates +r, h, z,, base vectors +e , e , e , and a "xed origin located at P F X the point of the "rst contact of the projectile with the target (see Fig. 2). Furthermore, it is convenient to consider another coordinate system attached to the projectile with coordinates +r, h, m,, base vectors +e , e , e ,, and a moving origin located at the distance of R /2 from the P F K projectile's tip. Here m"z!x(t)!R /2, R is the asymptotic value of the projectile's radius and the projectile's tip is located at z"x(t). More speci"cally, the projectile is taken to be an ovoid of Rankine of "nite length L which is characterized by the asymptotic radius R . In the following analysis, an equation of motion is derived of the form: [M#A(x, x )]xK "B(x, x )x #C(x, x ),
(1)
270
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 2. De"nition of the coordinate systems.
where the functions A(x, x ), B(x, x ) and C(x, x ) are determined analytically. Here a superposed dot denotes ordinary time di!erentiation. Also, the functions A, B and C correspond to what is commonly called the virtual mass of the target, form drag (quadratic in velocity), and plastic drag, respectively. Eq. (1) is solved numerically subject to the initial conditions x(0)"0,
(2a)
xR (0)"!; (2b) where ; ('0) is the projectile's initial velocity. In this solution x(t) is a monotonically decreasing function of time and integration is terminated when x vanishes. 2.2. Stress xeld in the target Since the target material is assumed to be incompressible, the conservation of mass and the balance of linear momentum equations become div v"0,
(3a)
ovR "! p#div r ,
(3b)
where the gradient operator and the divergence operator div are de"ned with respect to the present position of a material point, v is the absolute velocity of a material point in the target, and the Cauchy stress tensor r has been separated into its pressure p and its deviatoric part r , such that r"!pI#r ,
(4a)
r ) I"0.
(4b)
Also, I denotes the unit tensor, and the notation A ) B"tr(AB2) denotes the inner product between two second-order tensors A and B. At the instant that the projectile "rst touches the target's front surface, the target material begins to deform elastically. At some instant during the penetration process plastic yield initiates and an elastic}plastic boundary propagates away from the tip of the projectile. In the elastic region the
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
271
strains remain small, whereas in the plastic region the strains can be very large. Furthermore, for simplicity, it is assumed that the material response in the plastic region is rate-insensitive and rigid-plastic (with elastic strains being neglected). Under these conditions, the deviatoric stress can be approximated by the constitutive equation
2 > D, 3 (D ) D)
r "
(5)
where D is the symmetric part of the velocity gradient. The relationship of (5) to more general constitutive equations for elastic}viscoplastic metals has been discussed in [10]. Also, in the elastic region, it is assumed that the material response is linear elastic and isotropic so that for isochoric motion the deviatoric stress is related to the linear strain e by the expression r "2ke.
(6)
2.3. Velocity xeld in the target For axisymmetric irrotational motion, the velocity "eld is determined by a potential function
(r, z, t) such that v"
.
(7)
Consequently, the equations of motion (3) reduce to
"0,
o
(8a)
* 1 #
)
#p "div r , *t 2
(8b)
Following [2], the projectile is taken to be an ovoid of Rankine because the velocity "eld which satis"es the continuity equation and the condition of impenetrability at the projectile's surface is characterized by a simple combination of a single source and uniform #ow. Speci"cally, the ovoid is a body of revolution of length L whose lateral surface is de"ned by setting the radial coordinate r equal to the function R(m) which is de"ned by
R(m)"
R !m m , # (m#2R ) 2 2
(9a)
m"z!x(t)!R /2, (9b) where R is a constant that controls the radius of the projectile and m is measured relative to a "xed material point in the projectile (see Figs. 1 and 2), which locates the virtual source. For this case, the velocity potential , the radial component v , and the axial component v of the velocity "eld, P X which satisfy the continuity equation (8a), as well as the impenetrability condition at the ovoid surface (9) are given by
"xR
R 1 , 4 (m#r)
(10a)
272
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
R r v "!xR , P 4(m#r)
(10b)
v "0, F
(10c)
R m . v "!xR X 4(m#r)
(10d)
The model assumes that this velocity "eld is valid throughout the target for all of the di!erent penetration stages. Generally, the perforation process of a single plate can be divided into several stages, each described by a di!erent velocity "eld, as follows: embedding of the projectile in the target, dynamic penetration, bulge formation, bulge advancement, plug formation and exit. Since the velocity "eld (10) does not satisfy the boundary conditions at the free surfaces of the target it is most suited to describe the dynamic penetration stage. For relatively thick and ductile metal targets that are of interest in the present study, the dominant stage of longest duration is the dynamic penetration stage. Consequently, for such targets the other penetration stages are unlikely to signi"cantly e!ect the prediction of integral quantities such as the penetration depth and the residual velocity. This potential velocity "eld does not describe the detailed kinematics of shear near the projectile surface. Nevertheless, predictions of integral quantities such as the drag force acting on the projectile using this velocity "eld are expected to be good. This is mainly due to the fact that the decelerating force of the target is dominated by the normal stresses applied to the projectile's tip and not by the shear stresses. This aspect of the penetration process is also incorporated in many hydrocode programs since they do not take into account the shear (friction) forces between adjacent cells of the target and the projectile material. 2.4. Stress xeld in the rigid-plastic region The stress "eld in the rigid-plastic region is determined by the deviatoric stress and the pressure. Using the velocity "eld (10), the deviatoric stress tensor (5) can be written in the form [2]
r ">
1 [(m!2r)e e #(m#r)e e P P F F 3(m#r)
#(r!2m)e e !3rm(e e #e e )], (11) X X P X X P where denotes the tensor product. Furthermore, it can be shown using (11) that the divergence of the deviatoric stress is determined by the derivative of a scalar potential [2]
m#r . (12) R This simple result indicates that the balance of linear momentum (8b) can be integrated exactly pointwise in the rigid-plastic region to obtain an expression for the pressure of the form [2] div r " !> ln
R R m R p"f (t)#ox ! ! !oxK R 4(m#r) 4(m#r) 32(m#r) !> ln
m#r . R
(13)
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
273
Here f (t) is a function of time only, which must be determined by proper boundary and continuity conditions associated with the rigid-plastic region. It is important to emphasize that the simple result (12) is a consequence of the particular velocity "eld (10) that has been assumed for a projectile shaped as an ovoid of Rankine. For general projectile shapes, the simple result (12) does not hold and a more general description can be found in [2]. 2.5. Stress xeld in the elastic region In the elastic region of the target the strains remain small. The velocity "eld (10) is consistent with this fact because it possesses the correct asymptotic behavior, namely the velocity vanishes as r R or zP!R. Thus, the in"nitesimal strain e can be determined by integrating the approximate equation *e(r, z, t) "D, *t
(14)
where D is determined by the velocity "eld (10). Next, using expression (6) for the deviatoric stress and noting that div D"0, it follows that in the elastic region div r "2k div e"0.
(15)
Consequently, the balance of linear momentum (8b) yields an expression for the pressure in the elastic region of the form [2]
R R m R ! !oxK R , p"g(t)#ox ! 4(m#r) 4(m#r) 32(m#r)
(16)
where g(t) is a function of time only which must be determined by proper boundary and continuity conditions associated with the elastic region. 2.6. Method for determining the functions f (t) and g(t) From the previous subsections it is clear that in order to determine the pressure "eld in the rigid-plastic and elastic regions it is necessary to determine the functions f (t) and g(t), respectively. Once the pressure "eld is known, the stress "eld is determined, and the components of the drag force imposed on the projectile by the target can be calculated. Then, integrating over the portion of the projectile's surface that is in contact with the target material at a given instant, it is possible to determine the drag force acting on the projectile. 3. Boundary and continuity conditions for method I 3.1. General This method is similar to the one used in [2]. Speci"cally, the balance of linear momentum is satis"ed pointwise in the target region, whereas, the boundary and continuity conditions are only satis"ed approximately.
274
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
3.2. Stress-free boundary condition at the target's rear surface With the use of the assumed velocity "eld (10) it is not possible to satisfy the boundary condition that the target's rear surface remains stress-free pointwise. However, this condition can be approximated by requiring the axial stress component p , to vanish at r"0 and z"!H. Thus, it XX follows from (6), (14) and (16) that in the case of elastically deforming material at the rear surface of the target, the function g(t) is determined by the expression [2]
1 1 1 1 1 k xK #o ! # x # ! # , g(t)"!oR 4mM 4mM 32mM mM fM 2
(17a)
x 1 H ! ! , mM "! R 2 R
(17b)
H 1 fM "! ! . R 2
(17c)
For the case where the elastic}plastic boundary has reached the target's rear surface, the boundary condition for f (t) is proposed in a modi"ed weak form, which requires p to satisfy the XX condition
2n
PH
0
w(r)p (r,!H, t)r dr"0 XX
(18)
for the general case when the projectile protrudes through the rear surface. Here rH is the radial position of the point of intersection of the elastic}plastic boundary with the rear surface, R is the radius of the projectile at the intersection with the rear surface of the target, and w(r) is a weighting function which is described in more detail below. From Eqs. (4), (11) and (13), the axial stress in the plastic region on the rear surface of the target can be written in the form
R R m R ! #oxK R p (r,!H, t)"!f (t)!ox ! 4(m #r) XX 4(m #r) 32(m #r)
#> ln
r!2m m #r . # 3(m #r) R
(19)
Here, f (t) is a function of time only that must be determined by the boundary condition (18) at the target's rear surface. As can be seen from (19), given the velocity x and acceleration xK of the projectile at a certain time t, the function f (t) produces a uniform shift of the magnitude of the axial stress for any point on the rear surface.
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
275
For the special case when w(r)"1, the boundary condition (18) reduces to that used in [2], which requires the resultant force applied on the plastic region of the target's rear surface to vanish. From (11), (13) and (18) it follows that w(r ) ) r mM 1 dr oR P MH oP MH w(r ) ) r # dr 0 4(mM #r 0 4(mM #r ) 32(mM #r ) xK # f (t)" H x PM w(r ) ) r dr P MH w(r ) ) r dr 0 0 r !2mM dr >P MH w(r ) ) r ln(mM #r )# 3(mM #r ) 0 ! , (20) P MH w(r ) ) r dr 0 where a superposed bar denotes that the variable is a nondimensional variable which has been divided by R . Moreover, if the elastic}plastic boundary has reached the target's rear surface but the projectile's tip has not, then the value of R in (18) and (20) vanishes.
3.3. Continuity of traction on the elastic}plastic boundary The elastic}plastic boundary is determined by the values of r and z which cause the solution in the elastic region of the target to satisfy the von Mises yield condition G(r, z, x)"[ r ) r ] !>"k[6e ) e] !>"0.
(21)
For an exact solution it is necessary to require continuity of the traction vector at each point of the elastic}plastic boundary. However, for the approximate velocity "eld used here, this condition is approximated by requiring continuity of the axial stress at the intersection z"z of the elastic}plastic boundary with the axis r"0. Thus, for a given value of x(t), condition (21) evaluated at r"0 and z"z , yields the following equation [2]:
1 1 4> ! " , mM fM 3k
(22a)
1 x z , fM " ! "mM # R R 2
(22b)
z x 1 mM " ! ! . R R 2
(22c)
Next, equating the axial stress on both sides of the elastic}plastic boundary at r"0 and z"z , it can be shown that the function f (t) in the plastic region is given by
2 k 1 1 f (t)"g(t)#> ln(mM )! ! ! # . 3 2 mM fM
(23)
276
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
4. Boundary and continuity conditions for method II 4.1. General This method is similar to the one used in [6]. Speci"cally, the balance of linear momentum and the continuity and boundary conditions are satis"ed only along instantaneous streamlines. For the case of oblique penetration it was shown in [6] that for the assumed velocity "eld, the divergence of the deviatoric stress "eld is not the gradient of a scalar function as in the case of a normal impact (12). Nevertheless, as an approximation, the solution proposed in [6] determined the values of f (t) and g(t) by integration along the instantaneous streamline from a point on the surface of the projectile. This procedure causes di!erent values f (t) and g(t) to be obtained for each instantaneous streamline so that in this approximation these functions actually become dependent on both space and time instead of being dependent on time only, as in the solution of method I. In the present work, this procedure of integrating along instantaneous streamlines is used for the case of normal penetration. The advantage of this method is that the boundary and continuity conditions can be satis"ed pointwise at the intersections of each streamline with appropriate boundaries. However, the disadvantage of this method is that the linear momentum equation (3b) is no longer satis"ed pointwise in the entire elastic and plastic regions of the target, but only along instantaneous streamlines. This inability to satisfy both the linear momentum equation and the boundary and continuity conditions pointwise is due to the fact that the speci"c velocity "eld is approximate. To this end, the stream function t associated with the assumed axisymmetric velocity "eld (10) can be written in the form 1 *t , v "! P r *z
(24a)
1 *t , v " X r *r
(24b)
m xR R #1 , t(m, r)" 4 (m#r
(24c)
where the value of t has been taken to be zero on the axis (r"0) ahead of the projectile's tip where m(0. Since the value of the stream function t is constant along a streamline, a functional expression for the streamline that originates from an arbitrary point (for example point A in Fig. 3) on the projectile surface located at R can be expressed in the equivalent parametric forms
(25a)
4t !1 , x R
(25b)
4t !1 r(m, t )"m 1! x R
4t !1 m(r, t )"r x R
1!
4t !1 , x R
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
277
Fig. 3. Instantaneous streamlines relative to the "xed coordinate system.
t "t(m , R ),
(25c)
2R !R , m "m(R )" (25d) 2(R !R) where m is the value of m corresponding to point A located on the surface of the projectile at R , and t is value of the stream function t associated with the streamline that intersects point A. A plot of the instantaneous streamlines obtained using relations (25a) or (25b) is shown in Fig. 3. Note that (25d) follows from (9a) with R"R and m"m . 4.2. Boundary and continuity conditions In general, it is necessary to determine both of the functions f (t) and g(t) along each streamline. Typically, the function g(t) associated with the outer elastic region is determined by satisfying the boundary condition of zero traction on the target's front and rear surfaces, or by the condition that the stress asymptotically vanishes if the streamline approaches in"nity. Also, the function f (t) associated with the rigid}plastic region is determined by the continuity condition at the elastic}plastic boundary. For a semi-in"nite target it can be shown with the help of (4a), (6) and (16) that since the stress must vanish at in"nity, the function g(t) becomes g(t)"0.
(26)
A number of calculations for "nite thickness targets were performed in which g(t) was determined by the exact boundary conditions at the intersections of the streamlines and the boundaries. It was found that the simpli"ed form (26) produces reasonably accurate results as long as the thickness of the target is not too small. This is because the stresses in the elastic region are much smaller than those in the plastic region. Therefore, for simplicity, assumption (26) is used even for "nite thickness targets. This assumption was also used in [6] for the case of oblique penetration. Here, the function f (t) associated with a section of a speci"c instantaneous streamline contained in the plastic region is determined. This is accomplished by satisfying boundary and continuity
278
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 4. Approximations of the boundary and continuity conditions. The shaded regions are plastic and the unshaded regions are elastic or rigid.
conditions along this streamline. If the streamline intersects the elastic}plastic boundary (for example, see the streamline t in Fig. 4), then the continuity of the normal component of the traction vector on that surface is enforced (point Q in Fig. 4). Given the stress tensor r and the unit vector n normal to a speci"c surface, an expression for the normal component of the traction vector can be written using (4a) in the form p "n ) (rn)"!p#n ) (r n). LL
(27)
By equating the normal stress (27) on both sides of the elastic}plastic boundary at point Q , it is possible to obtain an expression for f (t) associated with this section of the instantaneous streamline t f (t)"> ln
m#r #n ) (r Nn)!n ) (r Cn). R
(28)
In this equation, m and r are the coordinates of point Q , n is a unit vector normal to the elastic}plastic boundary at point Q , and r N and r C are, respectively, the deviatoric stress tensors in the plastic and elastic regions at point Q . If the streamline in the plastic region intersects the stress-free rear surface (z"!H) [see the instantaneous streamline t in Fig. 4] or the front surface (z"0) [see the instantaneous streamline
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
279
t in Fig. 4] of the target, then the value of f (t) is determined by the condition that the normal component of the traction vector (p ) on that surface vanishes, so that XX m#r R m R f (t)"e ) (r Ne )#> ln !ox ! ! X X R 4(m#r) 32(m#r) R #oxK R (29) 4(m#r)
must be satis"ed at points Q or Q . Moreover, using (11) the expression for f (t) corresponding to the instantaneous streamlines t or t becomes R m R 1 ! #oxK R f (t)"!ox ! 4(m#r) 32(m#r) 4(m#r
r!2m m#r # , (30) 3(m#r) R where the coordinates m"m(t) and r"r(t) are evaluated at the points Q and Q , respectively. #> ln
5. Additional conditions for both methods I and II 5.1. Separation of the target material from the projectile This subsection de"nes the point where the target material separates from the projectile's surface. This separation point is determined by the radius R de"ned by R "min[R , R , R ], (31) * where R is the projectile's radius at the front surface, R is the projectile's radius at its tail. Also, * R is the radial position of the point on the projectile's surface satisfying the condition p (R )"0, (32) LL where p is the normal component of the traction acting on the projectile's surface (see Fig. 5). LL 5.2. Equation of motion of the projectile The equation of motion of the projectile in the axial direction is obtained from rigid body dynamics and can be written in the form F"MxK ,
(33a)
2 no R 2 M" N (1!RM ) !(1!RM ) #(1!RM )\ ! , * * * 3 3 2
(33b)
where M is the mass of the projectile, o its density and F is the resultant axial drag force acting in N the positive e direction [cf. Eq. (1)]. This drag force can be calculated by integrating the traction X vector t acting on the projectile's surface, which is given by the expression t"rn"p (R)n#p (R)s, LL LO
(34)
280
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 5. De"nition of the position of the separation point.
where n and s are the unit vectors in the directions of the outward normal and tangent to the projectile's surface, respectively. As in [2], the e!ect of the shear stress p on the value of F is LO neglected. Also, it is recalled that the normal stress p is compressive over the region R 4R4R LL and it is predicted quite accurately using the method of [2]. Thus,
L 0
p (R)R dR dh"!2n LL
0
p (R)R dR. (35) LL 0 0 In this equation, R denotes the projectile's radius at the rear surface, and R denotes the radius of the separation point of target material from the projectile's surface. Next, using the expression for p (obtained in [2]) for an arbitrary point on the projectile's surface, it follows that LL (R !R ox (R !3R)(R !R) #oxK p "!f (t)! LL 2 2 R (!8R #9R) R !R !> ln 4 #> . (36) 3(4R !3R) R For method I, f is a function of time only which is determined in the plastic region by the procedure discussed in Section 3. For method II, f is a function of time and space f (t, R) that is determined by the procedure described in Section 4. F" t ) e da+! X
6. Results Using the formulation described in the previous sections a computer program has been developed to numerically integrate (33), subject to the initial conditions (2) with ; being the impact velocity. The test case chosen for purposes of comparison between the di!erent solution schemes is that of a projectile (¸"56.5 mm, R "5.65 mm, M"86.6 g) made from Tungsten alloy, and a target made from aluminum (o"2703 kg/m, >"400 MPa, k"27.6 GPa). Two cases were examined: a semi-in"nite target; and a "nite thickness target (H"120 mm). Except where otherwise stated,
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
281
all of the analytical results use the weighting function w which is discussed in the last subsection. In the hydrocode computations, the rigid penetrator is modeled by a Lagrangian numerical processor while the target is modeled by an Eulerian numerical processor because of the large deformations that the target experiences. The target material is modeled using a Mie}Gruneisen shock equation for the pressure and an elastic perfectly plastic (constant yield strength) constitutive equation for deviatoric stress with the von Mises associated #ow rule. The Mie}Gruneisen equation constants of the target material used in all the hydrocode Autodyn2D computations are: c"5240 m/s (speed of sound), s"1.4 (shock parameter), C"1.97 (Gruneisen parameter). 6.1. Semi-inxnite target In order to eliminate the e!ect of the target's rear surface, and thus focus attention on the quasi-steady stage of the penetration process, it is conve nient to examine the case of penetration into a semi-in"nite target. Moreover, if the projectile penetrates into the target deep enough, then the e!ect of the target's front surface also becomes negligible. Fig. 6a shows the velocity of the
Fig. 6. Computational results for the case of a rigid projectile penetrating a semi-in"nite homogenous target. (a) Projectile's velocity versus penetration depth. (b) Drag force applied to the projectile versus penetration depth. Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
282
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
projectile versus the depth of penetration P and Fig. 6b shows the drag force F versus P. From Fig. 6b it can be seen that the drag force calculated by method II is signi"cantly higher than that calculated by Autodyn2D, whereas the value calculated by method I is very close to the Autodyn2D prediction. Figs. 7}9 show the comparison of various quantities predicted by the analytical model (methods I and II) and the code Autodyn2D. In each of these "gures the results are shown for the same impact conditions and the same penetration depth P"20R .
Fig. 7. Stress distribution along the projectile's surface. (a) Pressure versus the normalized radial distance. (b) Normal component of the stress versus the normalized radial distance. The penetration depth is P"20R . Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
283
Fig. 8. Predictions of (a) the axial velocity and (b) the pressure pro"les in the target along the symmetry axis beginning from the projectile's tip. The penetration depth is P"20R . Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
Figs. 7a and b show the pressure and the normal component of the traction vector imposed by the target on the projectile's surface, respectively, as calculated by the di!erent solution schemes. Here, the stress tensor is taken to be positive in tension, whereas the pressure is positive in compression. Consequently, the di!erences between the pressure p and the normal stress p are LL due to both the sign convention and the fact that p includes deviatoric stress. The values of these LL parameters in Autodyn2D were obtained by examining the target cells attached to the surface of the projectile. From these "gures it is quite clear that the main di!erence in the drag force is due to the di!erences in the pressure distributions calculated along the surface of the projectile. These di!erences increase with increasing radial distance from the axis. In the Autodyn2D computation, where the target is modeled with Eulerian cells and the projectile is modeled with Lagrangian cells, the cells interact across an Euler}Lagrange interface. In an Euler}Lagrange interface, the Lagrangian cells act as a geometric constraint on the Eulerian cells, while the Eulerian cells provide a pressure boundary condition for the Lagrangian cells (the tangential component of the stress vector applied on the projectile's surface is zero). At such an interface, slip can occur between the target material and the surface of the projectile. Thus, the only component of the stress vector contributing to the calculation of the drag force in Autodyn2D is
284
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 9. Comparison between the predictions of the analytical method I (which are the same as those of method II) and of Autodyn2d for: (a) the velocity "eld in the target and (b) the elastic}plastic boundary. The penetration depth is P"20R . Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
the normal component. In reality, there is also a contribution of the tangential component of the stress vector to this drag force. This approximation which neglects the tangential stress is common in most computation schemes. It results mainly from the di$culty in determining a proper value for the e!ective friction coe$cient associated with the interaction between the surfaces of the
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
285
projectile and the target. A partial justi"cation for such an approximation is the small contribution of the tangential stress component to the overall drag force relative to that of the normal component. Figs. 8a and b show the corresponding magnitude of the axial velocity and the pressure pro"les in the target along the symmetry axis beginning from the tip of the projectile. Furthermore, it is noted that on the axis r"0, the axial velocities of both methods I and II coincide. From Fig. 8a it is seen that the velocity at the projectile's tip is somewhat di!erent for each of these cases because the drag force as a function of time is di!erent. Consequently, the same penetration depth occurs at di!erent times for the di!erent simulations. More speci"cally, Fig. 8a shows that the analytical velocity "eld has a steeper gradient than the numerical one. Also, from Fig. 8b it is seen that the pressure "elds are nearly the same, except that Autodyn2D predicts a lower pressure near the projectile's tip. The simplifying assumption of incompressibility was used in the development of the analytical model. In contrast, the dilatational behavior of materials in Autodyn2D is modeled by specifying the Mie}Gruneisen equation of state which permits volume changes. In order to study the e!ect of the incompressibility assumption, another Autodyn2D computation (referred as `incompressiblea) was done using pseudo-aluminum having a bulk modulus 5 times greater than actual aluminum. From Fig. 8a it is seen that the velocity "eld is nearly insensitive to the later change although the resulting pressures in the target are somewhat higher (Fig. 8b) for the incompressible case. This leads to a higher drag force and a lower penetration depth. However, the di!erence is not more than 5%. Consequently, the incompressibility assumption used in the analytical model is reasonable. Fig. 9a shows a comparison of the velocity "elds predicted by Autodyn2D and the assumed "eld in the analytical models I or II at a speci"ed penetration depth. As mentioned previously, both Autodyn2D and the assumed analytical "eld enforce the impenetrability condition. Therefore in both cases the normal velocity component at the projectile surface is identical for the same projectile velocity. However, the projectile velocity predicted by Autodyn2D and the analytical methods are di!erent at the same penetration depth. Consequently, for comparison purposes the analytical "eld shown in Fig. 7a has been scaled so that the magnitude of the velocity vector at the tip of the projectile coincides with that in the Autodyn2D calculation. Moreover, it is noted that no restriction is placed on the tangential component of the velocity in the numerical calculation. Consequently, the tangential components of both velocity "elds di!er signi"cantly. Speci"cally, it can be seen that the analytical velocity "eld causes more radial motion than the numerical "eld. This is consistent with the fact that the numerical "eld includes some compressibility (unlike the analytical "eld) and that the analytical "eld has a steeper axial velocity gradient than the numerical one (Fig. 8a). In spite of these di!erences in the velocity "eld, it can be seen from Fig. 9b that the predicted elastic}plastic boundaries are in very close agreement. In this regard, it is noted that the analytical elastic}plastic boundary is a function only of the penetration depth and thus is the same for both methods I and II. 6.2. Finite thickness target Neither the analytical model nor the Autodyn2D calculation includes constitutive equations for modeling the various failure mechanisms that occur near the front and rear surfaces of actual
286
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
targets. Consequently, the theoretical predictions are actually limited to relatively thick ductile metal plates for which these failure modes are not expected to signi"cantly e!ect the prediction of integral quantities like penetration depth and residual velocity. The main di!erence between the Autodyn2D calculation and the theoretical method I is associated with the determination of the velocity "eld. The Autodyn2D calculation determines the velocity "eld that satis"es the boundary conditions of the free front and rear surfaces of the target, whereas the method I assumes a velocity "eld which does not satisfy these boundary conditions. The boundary condition (18) has been formulated in a weak form where the weighting function w can take di!erent forms. The objective of this subsection is to analyze the in#uence of three di!erent weighting functions w "1, 1 w " , r
(37a) (37b)
"p (r,!H, t)!p (rH,!H, t)" XX (37c) w " XX "p (0,!H, t)" XX For simplicity, only the theoretical method I is considered. The physical meaning of these weighting functions is as follows: For the case w"w , the boundary condition (18) corresponds to the physical condition that the total force acting on the plastic region of the rear surface vanishes, which is identical to the condition used in [2]. For the case w"w , the boundary condition (18) averages the axial stress by integrating only with respect of the radial coordinate instead of averaging with respect to the area, as is done for the weighting function w"w . Moreover, the weighting function w"w strongly emphasizes the magnitude of the axial stress near the axis (r"0) in the averaged boundary condition (18). The calculations associated with the weighting functions (37a)}(37c) are denoted by +w , w , w ,. Figs. 10a and b show the projectile velocity U and the drag force F, respectively, as functions of the penetration depth. The experimental value of the residual velocity of the projectile measured in this case was 537 m/s. According to the boundary condition (18), the weighting functions for the rear boundary begin to e!ect the results only once the elastic}plastic boundary reaches the rear target surface. The in#uence of these functions is clearly seen in Fig. 10a when the elastic}plastic boundary reaches the target's rear surface and the penetration depth is approximately 84 mm. Until this penetration depth, the solutions for the three weighting functions coincide. As can be seen from the solution for w"w in Fig. 10b, this condition causes the unphysical feature that the drag force exhibits a temporary increase when the elastic}plastic boundary reaches the rear surface. Only when the projectile starts perforating the rear surface does the drag force decrease sharply to zero. This should be contrasted with the predictions of Autodyn2D where the drag force monotonically decreases as the projectile approaches the rear surface. This unphysical behavior of the model motivated a search for other weighting functions that would yield an improved physical response of the drag force. Fig. 10b shows that the weighting functions w"w and w"w both satisfy this requirement. Fig. 11 plots the normalized radius versus the axial stress on the rear surface of the target for calculations using the three di!erent weighting functions (37a)}(37c). The sign of the stress
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
287
Fig. 10. Computational results for the case of a rigid projectile penetrating into a "nite thickness homogenous target. (a) The projectile's velocity versus penetration depth. (b) The drag force applied to the projectile versus penetration depth. The experimental value of the residual velocity of the projectile was 537 m/s. Projectile: M"86.6 g, ¸"56.5 mm, ; "845 m/s. Ovoid of Rankine: R "5.65 mm. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"120 mm.
corresponds to the usual convention that axial stress is positive in tension and negative in compression. All of the axial stress distribution curves in Fig. 11 are based on the calculation with w"w for the values of x and xK obtained at the moment the projectile's penetration depth is 113 mm (t"162 ls). Using these values in (19), as well as the boundary condition (18), di!erent values of f (t) are obtained for the di!erent weighting functions. As a result, the di!erent stress distributions shown in Fig. 11 are merely shifted by a constant value. It can be seen from Fig. 11 that the axial stress is negative (compression) on the axis r"0 and changes sign (to tension) at some radial distance from it, with all of these points being inside the plastic region on the target's rear surface. This fact corresponds to the requirement that a weighted average (18) of the axial stress vanishes. Furthermore, ampli"cation of the axial stress values close
288
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 11. Axial stress distribution on the target's rear surface (z"!H) resulting from di!erent weighting functions used in the boundary condition (18) of the stress free rear target surface. The penetration depth is P"113 mm. Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"120 mm. x , xK for all the distributions are obtained from the calculation using w
to the r"0 axis by the weighting function w"w results in a shift in the curve to the right (Fig. 11). This causes a decrease in the value of the compressive axial stresses near the axis r"0. This also causes an increase of the residual velocity of the projectile (Fig. 10a) as well as a reduction of the drag force imposed on the projectile (Fig. 10b). In particular, it can be seen from Fig. 11 that the e!ect of the axial stress on the axis (r"0) using the boundary condition (18) is emphasized more by the weighting function w than by w , which itself emphasizes this term more than the weighting function w . This is consistent with the results in Fig. 10b which show that the functions w and w reduce the drag force during the "nal stage of penetration more than the function w . Also, it can be seen from Fig. 10b that the drag forces calculated by method I with (w"w , w , w ) are higher than those calculated by Autodyn2D for the initial stage of penetra tion. The reason for this deviation is the in#uence of the modi"ed velocity "eld near the target's free front surface in the Autodyn2D calculation. Furthermore, it is seen that at the moment the projectile perforates the target's thickness (P"0.12 m), the drag force calculated by method I decreases sharply to zero, whereas that calculated by Autodyn2D decreases more moderately. This is due to bulging of the target's rear surface that is modeled by the Autodyn2D calculation and not in method I (for which the boundary condition is evaluated on the undeformed target's rear surface). Moreover, even if the boundary condition in method I was applied on the target's deformed rear surface, the results would not be signi"cantly di!erent because the assumed velocity "eld does not model this bulging phenomenon. These results indicate that although the target geometry is "xed in method I and the velocity "eld does not describe the exit stage correctly, it is possible to use the weighting function w to obtain
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
289
a reasonably good description of the drag force for this exit stage. Speci"cally, for this choice of the weighting function, the deviation (relative to Autodyn2D) of the residual velocity of the projectile is less than 0.5% (Fig. 10a). Also, it can be seen from Fig. 10b that method II produces the correct qualitative physical behavior of the drag force acting on the projectile by the target. This is due to the fact that method II satis"es the boundary and continuity conditions pointwise at the intersections of these boundaries with the streamlines. However, as was observed previously with the semi-in"nite target, method II overestimates the drag force acting on the projectile.
7. Conclusions Two di!erent theoretical approaches (method I [2] and method II [6]) were examined and compared with numerical calculations of Autodyn2D. Method I is based on satisfying the balance of linear momentum equation pointwise in the target, while approximating the boundary and continuity conditions only along the symmetry axis. On the other hand, method II satis"es the balance of linear momentum and the boundary and continuity conditions pointwise only along instantaneous streamlines relative to the target. For method I it was necessary to introduce approximate forms of the boundary and continuity conditions using appropriate weighting functions in order to obtain physically reasonable response during the exit phase of the projectile from a "nite thickness target. The physical implications of these approximations were investigated by comparing the results with numerical solutions of the exact formulation of the penetration problem using the hydrocode Autodyn2D. For the theoretical method I it follows that E The model produces the correct value of the drag force for penetration into a semi-in"nite target. E The solution for the exit stage of a "nite thickness target is sensitive to choice of the weighting function. E The behavior of the drag force applied to the projectile during the exit stage can be signi"cantly improved by a proper speci"cation of the weighting function. Moreover, for the theoretical method II it follows that E The model overpredicts the value of the drag force for penetration into a semi-in"nite target. E The model yields physically plausible response of the drag force during the exit stage of the penetration process since it satis"es the boundary and continuity conditions pointwise along streamlines. It is also important to emphasize that E The typical calculation time for a penetration process is reduced from several hours for Autodyn2D to only a few minutes for the analytical model. In summary, it can be concluded that both analytical methods I and II yield reliable results for integral parameters of the process such as the penetration depth, the ballistic limit and the residual velocity of the projectile. Moreover, it was shown in [2] for method I and in [6] for method II that these methods compare well with numerous experimental data. The results of the present work have shown that these methods are also capable of predicting such "ne characteristics as the shape
290
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
of the elastic}plastic boundary and the velocity "eld in the target. Furthermore, the temporal evolution of the penetration velocity and the drag force predicted by these methods agree quite well with the data obtained from simulations of Autodyn2D. Also, during the period when the projectile penetrates the target's rear surface, the time dependence of the drag force can be predicted by the analytical model when use is made of the weighting function w"w which emphasizes the value of axial stress on the symmetry axis. Furthermore, it is noted that the results presented for the particular case in this work are characteristic of those obtained for other values of the parameters.
Acknowledgements The work of M.B. Rubin and A.L. Yarin was partially supported by the fund for the promotion of research at Technion.
References [1] Backman ME, Goldsmith W. The mechanics of penetration of projectiles into targets. Int J Engng Sci 1978;16:1}99. [2] Yarin AL, Rubin MB, Roisman IV. Penetration of a rigid projectile into an elastic}plastic target of "nite thickness. Int J Impact Engng 1995;16:801}31. [3] Hill R. Cavitation and the in#uence of headshape in attack of thick targets by non-deforming projectiles. J Mech Phys Solids 1980;28:249}63. [4] Tate A. A comment on a paper by Awerbuch and Bodner concerning the mechanics of plate perforation by a projectile. Int J Engng Sci 1979;17:341}4. [5] Ravid M, Bodner SR. Dynamic perforation of viscoplastic plates by rigid projectiles. Int J Engng Sci 1983;21:577}91. [6] Roisman IV, Yarin AL, Rubin MB. Oblique penetration of a rigid projectile into an elastic}plastic target. Int J Impact Engng 1997;19:769}95. [7] Roisman IV, Weber K, Yarin AL, Hohler V, Rubin MB. Oblique penetration of a rigid projectile into a thick elastic}plastic target: theory and experiment. Int J Impact Engng 1999;22:707}26. [8] Autodyn2D. Interactive non-linear dynamic analysis software, theory manual. Revision 4.0, Century Dynamics, 1998. [9] Yossifon G. Penetration of a rigid projectile into a metal multilayer target and characterization of the debris cloud. MSc thesis, Technion * Israel Institute of Technology, Haifa, Israel, May 1999. [10] Rubin MB, Yarin AL. On the relationship between phenomenological models for elastic}viscoplastic metals and polymeric liquids. J Non-Newton Fluid Mech 1993;50:79}88 and 1995;57:321.
Penetration of a rigid projectile into a "nite thickness elastic}plastic target * comparison between theory and numerical computations G. Yossifon, M.B. Rubin*, A.L. Yarin Faculty of Mechanical Engineering, Technion, Israel Institute of Technology, 32000 Haifa, Israel Received 28 January 2000; received in revised form 19 June 2000
Abstract Recently, two di!erent analytical solutions have been developed for penetration of a rigid projectile into an elastic}plastic target of "nite thickness. One solution satis"es the balance of linear momentum pointwise in target region, but it satis"es the free-surface boundary conditions only in an integrated sense. The second solution satis"es the balance of linear momentum only along a "nite number of instantaneous streamlines, but it satis"es the boundary conditions exactly at each intersection of the streamlines with the boundary. The "rst solution is valid only for normal penetration whereas the second solution has been generalized for oblique penetration. The main objective of the present paper is to compare the predictions of these two theoretical approaches with a numerical solution obtained using the hydrocode Autodyn2D. By modifying the boundary condition used in the "rst method it is possible to obtain a reasonably accurate description of the penetration process (including the exit stage) with the computational time reduced from a few hours for Autodyn2D to only a few minutes for the analytical solution. 2001 Elsevier Science Ltd. All rights reserved.
1. Introduction A number of theoretical models have been developed for normal penetration of a rigid projectile into a "nite thickness elastic}plastic target, e.g. [1,2] and references therein. One theoretical approach is to use the kinematics of an inviscid #uid in order to generate an approximate velocity "eld in the target. The origins of this approach date back to [3}5]. Speci"cally, the velocity "eld is taken to be an irrotational "eld that is derivable from a velocity potential. In [2,6] the method of * Corresponding author. 0734-743X/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 7 3 4 - 7 4 3 X ( 0 0 ) 0 0 0 4 0 - 3
266
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Nomenclature A B C D F G H I L M P R R R R R R * R U ; Y c e,e, h e,e, h f g n r rH p s t t v w, w G x z z
e e n
virtual mass of the target coe$cient of the form drag plastic drag symmetric part of the velocity gradient (the rate of deformation tensor) resultant axial drag force acting on the projectile function in Eq. (21) target thickness the unit tensor projectile's length projectile's mass penetration depth radius of the projectile projectile's radius at the intersection with the rear surface of the target projectile's radius at the separation point projectile's radius at the intersection with the front surface of the target projectile's radius at the point on its surface which satis"es condition (32) projectile's radius at its tail asymptotic value of the projectile's radius velocity of the projectile initial velocity of the projectile constant yield strength in uniaxial stress speed of sound base vectors of a cylindrical polar coordinate system base vectors of the cylindrical polar coordinate system associated with the projectile function in Eq. (13) function in Eq. (16) unit vector normal to the elastic}plastic boundary or to the projectile's surface radial coordinate radial position of the point of intersection of the elastic}plastic boundary with the rear surface pressure shock parameter time traction (stress) vector absolute velocity of a material point in the target weighting functions axial coordinate of the projectile's tip axial coordinate intersection of the elastic}plastic boundary with the axis r"0.
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
267
Greek letters C e f , f h k m m , m o r, r s
t
Gruneisen parameter linear strain tensor dimensionless coordinates in Eqs. (17) and (22) angular coordinate shear modulus axial coordinate measured relative to the projectile dimensionless coordinates in Eqs. (17) and (22) mass density Cauchy stress tensor and its deviatoric part unit tangent vector to the projectile's surface velocity potential stream function
Subscripts A P n r z q
indicates indicates indicates indicates indicates indicates
that that that that that that
the the the the the the
parameter parameter parameter parameter parameter parameter
corresponds corresponds corresponds corresponds corresponds corresponds
to to to to to to
point A (see Fig. 3) the projectile the normal component the radial component the axial component the tangential component
Superscripts (e) (p) 䢇
-
indicates that the parameter corresponds to the elastic "eld indicates that the parameter corresponds to the plastic "eld denotes time di!erentiation denotes dimensionless parameters
singularities has been used to develop a velocity "eld that exactly satis"es the conditions of incompressibility in the target region and impenetrability of the projectile for any reasonable projectile shape. In [2] this method was applied to normal penetration and in [6] it was applied to oblique penetration. More speci"cally, when the projectile is taken to be an ovoid of Rankine (Fig. 1) and the impact is normal, the associated velocity "eld is quite simple and is known to be characterized by a combination of a single source and uniform #ow. The solution procedure assumes that the target material is incompressible. Then, the target region is divided into an elastic region ahead of the projectile, where the strains remain in"nitesimal, and a rigid-plastic region near the projectile's tip, where the strains can be very large. Furthermore, the rigid-plastic response is taken to be rate-independent with constant yield strength.
268
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 1. Actual scaled projectile's tip shaped as an ovoid of Rankine.
For normal penetration of an ovoid of Rankine [2], the velocity "eld is relatively simple and the divergence of the deviatoric stress becomes the gradient of a scalar potential. Consequently, the balance of linear momentum can be solved exactly in both the elastic and rigid-plastic regions to determine an expression for the pressure "eld and the stresses. Then, the e!ects of the free front and rear surfaces of the target, and the continuity of the surface tractions at the elastic/rigid-plastic boundary are satis"ed in an approximate integrated sense. This is because the assumed velocity "eld is not exact for the problem of penetration into a "nite target. Speci"cally, the continuity condition at the elastic}plastic boundary is satis"ed only at the intersection with the symmetry axis. The condition that the rear surface remains traction free is satis"ed either at the symmetry axis when the rear surface remains elastic, or in an integrated sense when it becomes plastic. The decelerating force applied to the projectile by the target is calculated by analytically integrating the surface tractions over the portion of the projectile's surface that is in contact with the target. Then, the resulting rigid-body equation of motion of the projectile is solved numerically to determine the penetration depth and velocity as functions of time. The theoretical predictions were compared with numerous experiments [2,6,7], and the agreement was found to be good even though the formulation uses no adjustable parameters. Only the mechanical properties of the materials, and projectile and target geometrical dimensions are used. For oblique penetration [6], the method of singularities yields a more complicated velocity "eld since an additional distribution of doublets is required to account for the projectile's rotation during the penetration process. For this case, the divergence of the deviatoric stress is no longer equal to the gradient of a scalar. Consequently, the balance of linear momentum cannot be solved exactly for the assumed velocity "eld. However, the balance of linear momentum can be solved in an approximate sense by integrating pointwise along a "nite set of instantaneous streamlines. Using this procedure it is possible to satisfy the continuity conditions on stress at the elastic}plastic boundary and the stress-free boundary conditions at portions of the target's front and rear surfaces, both at the intersections of each streamline with these surfaces. The method of [6] can also be applied to the simpler case of normal impact, as is discussed below. Here, method I denotes the procedure which satis"es the balance of linear momentum exactly, but the continuity conditions and boundary conditions only approximately [2]. Also, method II (the method of [6] applied to normal impact) denotes the procedure which satis"es the balance of linear momentum and the continuity and boundary conditions only along instantaneous stream-
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
269
lines. Consequently, the balance of linear momentum is satis"ed more accurately in method I than in method II. Whereas, the continuity and boundary conditions are satis"ed more accurately in method II than in method I. Thus, it is not clear for a particular problem of normal impact, which of these methods will predict more accurate values of the speci"c quantities of interest. The main objective of the present paper is to examine the accuracy of these two methods for normal penetration of a rigid projectile into rigid-plastic target. This is accomplished by comparing the predictions of the two approximate theoretical approaches with numerical simulations using the hydrocode Autodyn2D [8]. The Autodyn2D calculations include the compressibility of the target material and they determine the velocity "eld which is a numerical approximation of the exact solution of the balance of linear momentum and the continuity and boundary conditions. Consequently, it is possible to compare the theoretical and numerical predictions of such global quantities as the drag force and such local quantities as the stress acting on the projectile, the location of the elastic}plastic boundary and the velocity "eld. Also, a modi"cation of the stress-free boundary condition at the target's rear surface is proposed following the work in [9] to improve the theoretical predictions of method I. An outline of the paper is as follows: Section 2 reviews the basic equations of the two theoretical methods; Section 3 describes the boundary and continuity conditions for method I; Section 4 describes the boundary and continuity conditions for method II; Section 5 discusses additional conditions for both methods I and II; Section 6 presents the main results of example problems which compare the theoretical predictions of methods I and II with the results obtained in calculations using Autodyn2D; and Section 7 presents the main conclusions. 2. Basic equations 2.1. The coordinate system and the governing equations Consider normal penetration of a rigid projectile into an incompressible elastic}plastic target of thickness H, mass density o, shear modulus k and constant yield stress Y (in uniaxial stress). The projectile is assumed to have the shape of an ovoid of Rankine with radius R , length L, and mass M. Given these geometrical and material parameters, as well as the projectile's initial position and velocity, the theoretical model predicts the resulting penetration process. In particular, it predicts the instantaneous position and velocity of the projectile, as well as the "nal depth of penetration if the target defeats the projectile, or the residual velocity if the projectile perforates the target. Due to the axial symmetry of the problem, it is convenient to consider a cylindrical polar coordinate system with coordinates +r, h, z,, base vectors +e , e , e , and a "xed origin located at P F X the point of the "rst contact of the projectile with the target (see Fig. 2). Furthermore, it is convenient to consider another coordinate system attached to the projectile with coordinates +r, h, m,, base vectors +e , e , e ,, and a moving origin located at the distance of R /2 from the P F K projectile's tip. Here m"z!x(t)!R /2, R is the asymptotic value of the projectile's radius and the projectile's tip is located at z"x(t). More speci"cally, the projectile is taken to be an ovoid of Rankine of "nite length L which is characterized by the asymptotic radius R . In the following analysis, an equation of motion is derived of the form: [M#A(x, x )]xK "B(x, x )x #C(x, x ),
(1)
270
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 2. De"nition of the coordinate systems.
where the functions A(x, x ), B(x, x ) and C(x, x ) are determined analytically. Here a superposed dot denotes ordinary time di!erentiation. Also, the functions A, B and C correspond to what is commonly called the virtual mass of the target, form drag (quadratic in velocity), and plastic drag, respectively. Eq. (1) is solved numerically subject to the initial conditions x(0)"0,
(2a)
xR (0)"!; (2b) where ; ('0) is the projectile's initial velocity. In this solution x(t) is a monotonically decreasing function of time and integration is terminated when x vanishes. 2.2. Stress xeld in the target Since the target material is assumed to be incompressible, the conservation of mass and the balance of linear momentum equations become div v"0,
(3a)
ovR "! p#div r ,
(3b)
where the gradient operator and the divergence operator div are de"ned with respect to the present position of a material point, v is the absolute velocity of a material point in the target, and the Cauchy stress tensor r has been separated into its pressure p and its deviatoric part r , such that r"!pI#r ,
(4a)
r ) I"0.
(4b)
Also, I denotes the unit tensor, and the notation A ) B"tr(AB2) denotes the inner product between two second-order tensors A and B. At the instant that the projectile "rst touches the target's front surface, the target material begins to deform elastically. At some instant during the penetration process plastic yield initiates and an elastic}plastic boundary propagates away from the tip of the projectile. In the elastic region the
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
271
strains remain small, whereas in the plastic region the strains can be very large. Furthermore, for simplicity, it is assumed that the material response in the plastic region is rate-insensitive and rigid-plastic (with elastic strains being neglected). Under these conditions, the deviatoric stress can be approximated by the constitutive equation
2 > D, 3 (D ) D)
r "
(5)
where D is the symmetric part of the velocity gradient. The relationship of (5) to more general constitutive equations for elastic}viscoplastic metals has been discussed in [10]. Also, in the elastic region, it is assumed that the material response is linear elastic and isotropic so that for isochoric motion the deviatoric stress is related to the linear strain e by the expression r "2ke.
(6)
2.3. Velocity xeld in the target For axisymmetric irrotational motion, the velocity "eld is determined by a potential function
(r, z, t) such that v"
.
(7)
Consequently, the equations of motion (3) reduce to
"0,
o
(8a)
* 1 #
)
#p "div r , *t 2
(8b)
Following [2], the projectile is taken to be an ovoid of Rankine because the velocity "eld which satis"es the continuity equation and the condition of impenetrability at the projectile's surface is characterized by a simple combination of a single source and uniform #ow. Speci"cally, the ovoid is a body of revolution of length L whose lateral surface is de"ned by setting the radial coordinate r equal to the function R(m) which is de"ned by
R(m)"
R !m m , # (m#2R ) 2 2
(9a)
m"z!x(t)!R /2, (9b) where R is a constant that controls the radius of the projectile and m is measured relative to a "xed material point in the projectile (see Figs. 1 and 2), which locates the virtual source. For this case, the velocity potential , the radial component v , and the axial component v of the velocity "eld, P X which satisfy the continuity equation (8a), as well as the impenetrability condition at the ovoid surface (9) are given by
"xR
R 1 , 4 (m#r)
(10a)
272
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
R r v "!xR , P 4(m#r)
(10b)
v "0, F
(10c)
R m . v "!xR X 4(m#r)
(10d)
The model assumes that this velocity "eld is valid throughout the target for all of the di!erent penetration stages. Generally, the perforation process of a single plate can be divided into several stages, each described by a di!erent velocity "eld, as follows: embedding of the projectile in the target, dynamic penetration, bulge formation, bulge advancement, plug formation and exit. Since the velocity "eld (10) does not satisfy the boundary conditions at the free surfaces of the target it is most suited to describe the dynamic penetration stage. For relatively thick and ductile metal targets that are of interest in the present study, the dominant stage of longest duration is the dynamic penetration stage. Consequently, for such targets the other penetration stages are unlikely to signi"cantly e!ect the prediction of integral quantities such as the penetration depth and the residual velocity. This potential velocity "eld does not describe the detailed kinematics of shear near the projectile surface. Nevertheless, predictions of integral quantities such as the drag force acting on the projectile using this velocity "eld are expected to be good. This is mainly due to the fact that the decelerating force of the target is dominated by the normal stresses applied to the projectile's tip and not by the shear stresses. This aspect of the penetration process is also incorporated in many hydrocode programs since they do not take into account the shear (friction) forces between adjacent cells of the target and the projectile material. 2.4. Stress xeld in the rigid-plastic region The stress "eld in the rigid-plastic region is determined by the deviatoric stress and the pressure. Using the velocity "eld (10), the deviatoric stress tensor (5) can be written in the form [2]
r ">
1 [(m!2r)e e #(m#r)e e P P F F 3(m#r)
#(r!2m)e e !3rm(e e #e e )], (11) X X P X X P where denotes the tensor product. Furthermore, it can be shown using (11) that the divergence of the deviatoric stress is determined by the derivative of a scalar potential [2]
m#r . (12) R This simple result indicates that the balance of linear momentum (8b) can be integrated exactly pointwise in the rigid-plastic region to obtain an expression for the pressure of the form [2] div r " !> ln
R R m R p"f (t)#ox ! ! !oxK R 4(m#r) 4(m#r) 32(m#r) !> ln
m#r . R
(13)
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
273
Here f (t) is a function of time only, which must be determined by proper boundary and continuity conditions associated with the rigid-plastic region. It is important to emphasize that the simple result (12) is a consequence of the particular velocity "eld (10) that has been assumed for a projectile shaped as an ovoid of Rankine. For general projectile shapes, the simple result (12) does not hold and a more general description can be found in [2]. 2.5. Stress xeld in the elastic region In the elastic region of the target the strains remain small. The velocity "eld (10) is consistent with this fact because it possesses the correct asymptotic behavior, namely the velocity vanishes as r R or zP!R. Thus, the in"nitesimal strain e can be determined by integrating the approximate equation *e(r, z, t) "D, *t
(14)
where D is determined by the velocity "eld (10). Next, using expression (6) for the deviatoric stress and noting that div D"0, it follows that in the elastic region div r "2k div e"0.
(15)
Consequently, the balance of linear momentum (8b) yields an expression for the pressure in the elastic region of the form [2]
R R m R ! !oxK R , p"g(t)#ox ! 4(m#r) 4(m#r) 32(m#r)
(16)
where g(t) is a function of time only which must be determined by proper boundary and continuity conditions associated with the elastic region. 2.6. Method for determining the functions f (t) and g(t) From the previous subsections it is clear that in order to determine the pressure "eld in the rigid-plastic and elastic regions it is necessary to determine the functions f (t) and g(t), respectively. Once the pressure "eld is known, the stress "eld is determined, and the components of the drag force imposed on the projectile by the target can be calculated. Then, integrating over the portion of the projectile's surface that is in contact with the target material at a given instant, it is possible to determine the drag force acting on the projectile. 3. Boundary and continuity conditions for method I 3.1. General This method is similar to the one used in [2]. Speci"cally, the balance of linear momentum is satis"ed pointwise in the target region, whereas, the boundary and continuity conditions are only satis"ed approximately.
274
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
3.2. Stress-free boundary condition at the target's rear surface With the use of the assumed velocity "eld (10) it is not possible to satisfy the boundary condition that the target's rear surface remains stress-free pointwise. However, this condition can be approximated by requiring the axial stress component p , to vanish at r"0 and z"!H. Thus, it XX follows from (6), (14) and (16) that in the case of elastically deforming material at the rear surface of the target, the function g(t) is determined by the expression [2]
1 1 1 1 1 k xK #o ! # x # ! # , g(t)"!oR 4mM 4mM 32mM mM fM 2
(17a)
x 1 H ! ! , mM "! R 2 R
(17b)
H 1 fM "! ! . R 2
(17c)
For the case where the elastic}plastic boundary has reached the target's rear surface, the boundary condition for f (t) is proposed in a modi"ed weak form, which requires p to satisfy the XX condition
2n
PH
0
w(r)p (r,!H, t)r dr"0 XX
(18)
for the general case when the projectile protrudes through the rear surface. Here rH is the radial position of the point of intersection of the elastic}plastic boundary with the rear surface, R is the radius of the projectile at the intersection with the rear surface of the target, and w(r) is a weighting function which is described in more detail below. From Eqs. (4), (11) and (13), the axial stress in the plastic region on the rear surface of the target can be written in the form
R R m R ! #oxK R p (r,!H, t)"!f (t)!ox ! 4(m #r) XX 4(m #r) 32(m #r)
#> ln
r!2m m #r . # 3(m #r) R
(19)
Here, f (t) is a function of time only that must be determined by the boundary condition (18) at the target's rear surface. As can be seen from (19), given the velocity x and acceleration xK of the projectile at a certain time t, the function f (t) produces a uniform shift of the magnitude of the axial stress for any point on the rear surface.
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
275
For the special case when w(r)"1, the boundary condition (18) reduces to that used in [2], which requires the resultant force applied on the plastic region of the target's rear surface to vanish. From (11), (13) and (18) it follows that w(r ) ) r mM 1 dr oR P MH oP MH w(r ) ) r # dr 0 4(mM #r 0 4(mM #r ) 32(mM #r ) xK # f (t)" H x PM w(r ) ) r dr P MH w(r ) ) r dr 0 0 r !2mM dr >P MH w(r ) ) r ln(mM #r )# 3(mM #r ) 0 ! , (20) P MH w(r ) ) r dr 0 where a superposed bar denotes that the variable is a nondimensional variable which has been divided by R . Moreover, if the elastic}plastic boundary has reached the target's rear surface but the projectile's tip has not, then the value of R in (18) and (20) vanishes.
3.3. Continuity of traction on the elastic}plastic boundary The elastic}plastic boundary is determined by the values of r and z which cause the solution in the elastic region of the target to satisfy the von Mises yield condition G(r, z, x)"[ r ) r ] !>"k[6e ) e] !>"0.
(21)
For an exact solution it is necessary to require continuity of the traction vector at each point of the elastic}plastic boundary. However, for the approximate velocity "eld used here, this condition is approximated by requiring continuity of the axial stress at the intersection z"z of the elastic}plastic boundary with the axis r"0. Thus, for a given value of x(t), condition (21) evaluated at r"0 and z"z , yields the following equation [2]:
1 1 4> ! " , mM fM 3k
(22a)
1 x z , fM " ! "mM # R R 2
(22b)
z x 1 mM " ! ! . R R 2
(22c)
Next, equating the axial stress on both sides of the elastic}plastic boundary at r"0 and z"z , it can be shown that the function f (t) in the plastic region is given by
2 k 1 1 f (t)"g(t)#> ln(mM )! ! ! # . 3 2 mM fM
(23)
276
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
4. Boundary and continuity conditions for method II 4.1. General This method is similar to the one used in [6]. Speci"cally, the balance of linear momentum and the continuity and boundary conditions are satis"ed only along instantaneous streamlines. For the case of oblique penetration it was shown in [6] that for the assumed velocity "eld, the divergence of the deviatoric stress "eld is not the gradient of a scalar function as in the case of a normal impact (12). Nevertheless, as an approximation, the solution proposed in [6] determined the values of f (t) and g(t) by integration along the instantaneous streamline from a point on the surface of the projectile. This procedure causes di!erent values f (t) and g(t) to be obtained for each instantaneous streamline so that in this approximation these functions actually become dependent on both space and time instead of being dependent on time only, as in the solution of method I. In the present work, this procedure of integrating along instantaneous streamlines is used for the case of normal penetration. The advantage of this method is that the boundary and continuity conditions can be satis"ed pointwise at the intersections of each streamline with appropriate boundaries. However, the disadvantage of this method is that the linear momentum equation (3b) is no longer satis"ed pointwise in the entire elastic and plastic regions of the target, but only along instantaneous streamlines. This inability to satisfy both the linear momentum equation and the boundary and continuity conditions pointwise is due to the fact that the speci"c velocity "eld is approximate. To this end, the stream function t associated with the assumed axisymmetric velocity "eld (10) can be written in the form 1 *t , v "! P r *z
(24a)
1 *t , v " X r *r
(24b)
m xR R #1 , t(m, r)" 4 (m#r
(24c)
where the value of t has been taken to be zero on the axis (r"0) ahead of the projectile's tip where m(0. Since the value of the stream function t is constant along a streamline, a functional expression for the streamline that originates from an arbitrary point (for example point A in Fig. 3) on the projectile surface located at R can be expressed in the equivalent parametric forms
(25a)
4t !1 , x R
(25b)
4t !1 r(m, t )"m 1! x R
4t !1 m(r, t )"r x R
1!
4t !1 , x R
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
277
Fig. 3. Instantaneous streamlines relative to the "xed coordinate system.
t "t(m , R ),
(25c)
2R !R , m "m(R )" (25d) 2(R !R) where m is the value of m corresponding to point A located on the surface of the projectile at R , and t is value of the stream function t associated with the streamline that intersects point A. A plot of the instantaneous streamlines obtained using relations (25a) or (25b) is shown in Fig. 3. Note that (25d) follows from (9a) with R"R and m"m . 4.2. Boundary and continuity conditions In general, it is necessary to determine both of the functions f (t) and g(t) along each streamline. Typically, the function g(t) associated with the outer elastic region is determined by satisfying the boundary condition of zero traction on the target's front and rear surfaces, or by the condition that the stress asymptotically vanishes if the streamline approaches in"nity. Also, the function f (t) associated with the rigid}plastic region is determined by the continuity condition at the elastic}plastic boundary. For a semi-in"nite target it can be shown with the help of (4a), (6) and (16) that since the stress must vanish at in"nity, the function g(t) becomes g(t)"0.
(26)
A number of calculations for "nite thickness targets were performed in which g(t) was determined by the exact boundary conditions at the intersections of the streamlines and the boundaries. It was found that the simpli"ed form (26) produces reasonably accurate results as long as the thickness of the target is not too small. This is because the stresses in the elastic region are much smaller than those in the plastic region. Therefore, for simplicity, assumption (26) is used even for "nite thickness targets. This assumption was also used in [6] for the case of oblique penetration. Here, the function f (t) associated with a section of a speci"c instantaneous streamline contained in the plastic region is determined. This is accomplished by satisfying boundary and continuity
278
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 4. Approximations of the boundary and continuity conditions. The shaded regions are plastic and the unshaded regions are elastic or rigid.
conditions along this streamline. If the streamline intersects the elastic}plastic boundary (for example, see the streamline t in Fig. 4), then the continuity of the normal component of the traction vector on that surface is enforced (point Q in Fig. 4). Given the stress tensor r and the unit vector n normal to a speci"c surface, an expression for the normal component of the traction vector can be written using (4a) in the form p "n ) (rn)"!p#n ) (r n). LL
(27)
By equating the normal stress (27) on both sides of the elastic}plastic boundary at point Q , it is possible to obtain an expression for f (t) associated with this section of the instantaneous streamline t f (t)"> ln
m#r #n ) (r Nn)!n ) (r Cn). R
(28)
In this equation, m and r are the coordinates of point Q , n is a unit vector normal to the elastic}plastic boundary at point Q , and r N and r C are, respectively, the deviatoric stress tensors in the plastic and elastic regions at point Q . If the streamline in the plastic region intersects the stress-free rear surface (z"!H) [see the instantaneous streamline t in Fig. 4] or the front surface (z"0) [see the instantaneous streamline
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
279
t in Fig. 4] of the target, then the value of f (t) is determined by the condition that the normal component of the traction vector (p ) on that surface vanishes, so that XX m#r R m R f (t)"e ) (r Ne )#> ln !ox ! ! X X R 4(m#r) 32(m#r) R #oxK R (29) 4(m#r)
must be satis"ed at points Q or Q . Moreover, using (11) the expression for f (t) corresponding to the instantaneous streamlines t or t becomes R m R 1 ! #oxK R f (t)"!ox ! 4(m#r) 32(m#r) 4(m#r
r!2m m#r # , (30) 3(m#r) R where the coordinates m"m(t) and r"r(t) are evaluated at the points Q and Q , respectively. #> ln
5. Additional conditions for both methods I and II 5.1. Separation of the target material from the projectile This subsection de"nes the point where the target material separates from the projectile's surface. This separation point is determined by the radius R de"ned by R "min[R , R , R ], (31) * where R is the projectile's radius at the front surface, R is the projectile's radius at its tail. Also, * R is the radial position of the point on the projectile's surface satisfying the condition p (R )"0, (32) LL where p is the normal component of the traction acting on the projectile's surface (see Fig. 5). LL 5.2. Equation of motion of the projectile The equation of motion of the projectile in the axial direction is obtained from rigid body dynamics and can be written in the form F"MxK ,
(33a)
2 no R 2 M" N (1!RM ) !(1!RM ) #(1!RM )\ ! , * * * 3 3 2
(33b)
where M is the mass of the projectile, o its density and F is the resultant axial drag force acting in N the positive e direction [cf. Eq. (1)]. This drag force can be calculated by integrating the traction X vector t acting on the projectile's surface, which is given by the expression t"rn"p (R)n#p (R)s, LL LO
(34)
280
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 5. De"nition of the position of the separation point.
where n and s are the unit vectors in the directions of the outward normal and tangent to the projectile's surface, respectively. As in [2], the e!ect of the shear stress p on the value of F is LO neglected. Also, it is recalled that the normal stress p is compressive over the region R 4R4R LL and it is predicted quite accurately using the method of [2]. Thus,
L 0
p (R)R dR dh"!2n LL
0
p (R)R dR. (35) LL 0 0 In this equation, R denotes the projectile's radius at the rear surface, and R denotes the radius of the separation point of target material from the projectile's surface. Next, using the expression for p (obtained in [2]) for an arbitrary point on the projectile's surface, it follows that LL (R !R ox (R !3R)(R !R) #oxK p "!f (t)! LL 2 2 R (!8R #9R) R !R !> ln 4 #> . (36) 3(4R !3R) R For method I, f is a function of time only which is determined in the plastic region by the procedure discussed in Section 3. For method II, f is a function of time and space f (t, R) that is determined by the procedure described in Section 4. F" t ) e da+! X
6. Results Using the formulation described in the previous sections a computer program has been developed to numerically integrate (33), subject to the initial conditions (2) with ; being the impact velocity. The test case chosen for purposes of comparison between the di!erent solution schemes is that of a projectile (¸"56.5 mm, R "5.65 mm, M"86.6 g) made from Tungsten alloy, and a target made from aluminum (o"2703 kg/m, >"400 MPa, k"27.6 GPa). Two cases were examined: a semi-in"nite target; and a "nite thickness target (H"120 mm). Except where otherwise stated,
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
281
all of the analytical results use the weighting function w which is discussed in the last subsection. In the hydrocode computations, the rigid penetrator is modeled by a Lagrangian numerical processor while the target is modeled by an Eulerian numerical processor because of the large deformations that the target experiences. The target material is modeled using a Mie}Gruneisen shock equation for the pressure and an elastic perfectly plastic (constant yield strength) constitutive equation for deviatoric stress with the von Mises associated #ow rule. The Mie}Gruneisen equation constants of the target material used in all the hydrocode Autodyn2D computations are: c"5240 m/s (speed of sound), s"1.4 (shock parameter), C"1.97 (Gruneisen parameter). 6.1. Semi-inxnite target In order to eliminate the e!ect of the target's rear surface, and thus focus attention on the quasi-steady stage of the penetration process, it is conve nient to examine the case of penetration into a semi-in"nite target. Moreover, if the projectile penetrates into the target deep enough, then the e!ect of the target's front surface also becomes negligible. Fig. 6a shows the velocity of the
Fig. 6. Computational results for the case of a rigid projectile penetrating a semi-in"nite homogenous target. (a) Projectile's velocity versus penetration depth. (b) Drag force applied to the projectile versus penetration depth. Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
282
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
projectile versus the depth of penetration P and Fig. 6b shows the drag force F versus P. From Fig. 6b it can be seen that the drag force calculated by method II is signi"cantly higher than that calculated by Autodyn2D, whereas the value calculated by method I is very close to the Autodyn2D prediction. Figs. 7}9 show the comparison of various quantities predicted by the analytical model (methods I and II) and the code Autodyn2D. In each of these "gures the results are shown for the same impact conditions and the same penetration depth P"20R .
Fig. 7. Stress distribution along the projectile's surface. (a) Pressure versus the normalized radial distance. (b) Normal component of the stress versus the normalized radial distance. The penetration depth is P"20R . Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
283
Fig. 8. Predictions of (a) the axial velocity and (b) the pressure pro"les in the target along the symmetry axis beginning from the projectile's tip. The penetration depth is P"20R . Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
Figs. 7a and b show the pressure and the normal component of the traction vector imposed by the target on the projectile's surface, respectively, as calculated by the di!erent solution schemes. Here, the stress tensor is taken to be positive in tension, whereas the pressure is positive in compression. Consequently, the di!erences between the pressure p and the normal stress p are LL due to both the sign convention and the fact that p includes deviatoric stress. The values of these LL parameters in Autodyn2D were obtained by examining the target cells attached to the surface of the projectile. From these "gures it is quite clear that the main di!erence in the drag force is due to the di!erences in the pressure distributions calculated along the surface of the projectile. These di!erences increase with increasing radial distance from the axis. In the Autodyn2D computation, where the target is modeled with Eulerian cells and the projectile is modeled with Lagrangian cells, the cells interact across an Euler}Lagrange interface. In an Euler}Lagrange interface, the Lagrangian cells act as a geometric constraint on the Eulerian cells, while the Eulerian cells provide a pressure boundary condition for the Lagrangian cells (the tangential component of the stress vector applied on the projectile's surface is zero). At such an interface, slip can occur between the target material and the surface of the projectile. Thus, the only component of the stress vector contributing to the calculation of the drag force in Autodyn2D is
284
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 9. Comparison between the predictions of the analytical method I (which are the same as those of method II) and of Autodyn2d for: (a) the velocity "eld in the target and (b) the elastic}plastic boundary. The penetration depth is P"20R . Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"(semi-in"nite target).
the normal component. In reality, there is also a contribution of the tangential component of the stress vector to this drag force. This approximation which neglects the tangential stress is common in most computation schemes. It results mainly from the di$culty in determining a proper value for the e!ective friction coe$cient associated with the interaction between the surfaces of the
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
285
projectile and the target. A partial justi"cation for such an approximation is the small contribution of the tangential stress component to the overall drag force relative to that of the normal component. Figs. 8a and b show the corresponding magnitude of the axial velocity and the pressure pro"les in the target along the symmetry axis beginning from the tip of the projectile. Furthermore, it is noted that on the axis r"0, the axial velocities of both methods I and II coincide. From Fig. 8a it is seen that the velocity at the projectile's tip is somewhat di!erent for each of these cases because the drag force as a function of time is di!erent. Consequently, the same penetration depth occurs at di!erent times for the di!erent simulations. More speci"cally, Fig. 8a shows that the analytical velocity "eld has a steeper gradient than the numerical one. Also, from Fig. 8b it is seen that the pressure "elds are nearly the same, except that Autodyn2D predicts a lower pressure near the projectile's tip. The simplifying assumption of incompressibility was used in the development of the analytical model. In contrast, the dilatational behavior of materials in Autodyn2D is modeled by specifying the Mie}Gruneisen equation of state which permits volume changes. In order to study the e!ect of the incompressibility assumption, another Autodyn2D computation (referred as `incompressiblea) was done using pseudo-aluminum having a bulk modulus 5 times greater than actual aluminum. From Fig. 8a it is seen that the velocity "eld is nearly insensitive to the later change although the resulting pressures in the target are somewhat higher (Fig. 8b) for the incompressible case. This leads to a higher drag force and a lower penetration depth. However, the di!erence is not more than 5%. Consequently, the incompressibility assumption used in the analytical model is reasonable. Fig. 9a shows a comparison of the velocity "elds predicted by Autodyn2D and the assumed "eld in the analytical models I or II at a speci"ed penetration depth. As mentioned previously, both Autodyn2D and the assumed analytical "eld enforce the impenetrability condition. Therefore in both cases the normal velocity component at the projectile surface is identical for the same projectile velocity. However, the projectile velocity predicted by Autodyn2D and the analytical methods are di!erent at the same penetration depth. Consequently, for comparison purposes the analytical "eld shown in Fig. 7a has been scaled so that the magnitude of the velocity vector at the tip of the projectile coincides with that in the Autodyn2D calculation. Moreover, it is noted that no restriction is placed on the tangential component of the velocity in the numerical calculation. Consequently, the tangential components of both velocity "elds di!er signi"cantly. Speci"cally, it can be seen that the analytical velocity "eld causes more radial motion than the numerical "eld. This is consistent with the fact that the numerical "eld includes some compressibility (unlike the analytical "eld) and that the analytical "eld has a steeper axial velocity gradient than the numerical one (Fig. 8a). In spite of these di!erences in the velocity "eld, it can be seen from Fig. 9b that the predicted elastic}plastic boundaries are in very close agreement. In this regard, it is noted that the analytical elastic}plastic boundary is a function only of the penetration depth and thus is the same for both methods I and II. 6.2. Finite thickness target Neither the analytical model nor the Autodyn2D calculation includes constitutive equations for modeling the various failure mechanisms that occur near the front and rear surfaces of actual
286
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
targets. Consequently, the theoretical predictions are actually limited to relatively thick ductile metal plates for which these failure modes are not expected to signi"cantly e!ect the prediction of integral quantities like penetration depth and residual velocity. The main di!erence between the Autodyn2D calculation and the theoretical method I is associated with the determination of the velocity "eld. The Autodyn2D calculation determines the velocity "eld that satis"es the boundary conditions of the free front and rear surfaces of the target, whereas the method I assumes a velocity "eld which does not satisfy these boundary conditions. The boundary condition (18) has been formulated in a weak form where the weighting function w can take di!erent forms. The objective of this subsection is to analyze the in#uence of three di!erent weighting functions w "1, 1 w " , r
(37a) (37b)
"p (r,!H, t)!p (rH,!H, t)" XX (37c) w " XX "p (0,!H, t)" XX For simplicity, only the theoretical method I is considered. The physical meaning of these weighting functions is as follows: For the case w"w , the boundary condition (18) corresponds to the physical condition that the total force acting on the plastic region of the rear surface vanishes, which is identical to the condition used in [2]. For the case w"w , the boundary condition (18) averages the axial stress by integrating only with respect of the radial coordinate instead of averaging with respect to the area, as is done for the weighting function w"w . Moreover, the weighting function w"w strongly emphasizes the magnitude of the axial stress near the axis (r"0) in the averaged boundary condition (18). The calculations associated with the weighting functions (37a)}(37c) are denoted by +w , w , w ,. Figs. 10a and b show the projectile velocity U and the drag force F, respectively, as functions of the penetration depth. The experimental value of the residual velocity of the projectile measured in this case was 537 m/s. According to the boundary condition (18), the weighting functions for the rear boundary begin to e!ect the results only once the elastic}plastic boundary reaches the rear target surface. The in#uence of these functions is clearly seen in Fig. 10a when the elastic}plastic boundary reaches the target's rear surface and the penetration depth is approximately 84 mm. Until this penetration depth, the solutions for the three weighting functions coincide. As can be seen from the solution for w"w in Fig. 10b, this condition causes the unphysical feature that the drag force exhibits a temporary increase when the elastic}plastic boundary reaches the rear surface. Only when the projectile starts perforating the rear surface does the drag force decrease sharply to zero. This should be contrasted with the predictions of Autodyn2D where the drag force monotonically decreases as the projectile approaches the rear surface. This unphysical behavior of the model motivated a search for other weighting functions that would yield an improved physical response of the drag force. Fig. 10b shows that the weighting functions w"w and w"w both satisfy this requirement. Fig. 11 plots the normalized radius versus the axial stress on the rear surface of the target for calculations using the three di!erent weighting functions (37a)}(37c). The sign of the stress
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
287
Fig. 10. Computational results for the case of a rigid projectile penetrating into a "nite thickness homogenous target. (a) The projectile's velocity versus penetration depth. (b) The drag force applied to the projectile versus penetration depth. The experimental value of the residual velocity of the projectile was 537 m/s. Projectile: M"86.6 g, ¸"56.5 mm, ; "845 m/s. Ovoid of Rankine: R "5.65 mm. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"120 mm.
corresponds to the usual convention that axial stress is positive in tension and negative in compression. All of the axial stress distribution curves in Fig. 11 are based on the calculation with w"w for the values of x and xK obtained at the moment the projectile's penetration depth is 113 mm (t"162 ls). Using these values in (19), as well as the boundary condition (18), di!erent values of f (t) are obtained for the di!erent weighting functions. As a result, the di!erent stress distributions shown in Fig. 11 are merely shifted by a constant value. It can be seen from Fig. 11 that the axial stress is negative (compression) on the axis r"0 and changes sign (to tension) at some radial distance from it, with all of these points being inside the plastic region on the target's rear surface. This fact corresponds to the requirement that a weighted average (18) of the axial stress vanishes. Furthermore, ampli"cation of the axial stress values close
288
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
Fig. 11. Axial stress distribution on the target's rear surface (z"!H) resulting from di!erent weighting functions used in the boundary condition (18) of the stress free rear target surface. The penetration depth is P"113 mm. Projectile: M"86.6 g, ¸"56.5 mm, R "5.65 mm, ; "845 m/s. Target: o"2703 kg/m, >"400 MPa, k"27.6 GPa, H"120 mm. x , xK for all the distributions are obtained from the calculation using w
to the r"0 axis by the weighting function w"w results in a shift in the curve to the right (Fig. 11). This causes a decrease in the value of the compressive axial stresses near the axis r"0. This also causes an increase of the residual velocity of the projectile (Fig. 10a) as well as a reduction of the drag force imposed on the projectile (Fig. 10b). In particular, it can be seen from Fig. 11 that the e!ect of the axial stress on the axis (r"0) using the boundary condition (18) is emphasized more by the weighting function w than by w , which itself emphasizes this term more than the weighting function w . This is consistent with the results in Fig. 10b which show that the functions w and w reduce the drag force during the "nal stage of penetration more than the function w . Also, it can be seen from Fig. 10b that the drag forces calculated by method I with (w"w , w , w ) are higher than those calculated by Autodyn2D for the initial stage of penetra tion. The reason for this deviation is the in#uence of the modi"ed velocity "eld near the target's free front surface in the Autodyn2D calculation. Furthermore, it is seen that at the moment the projectile perforates the target's thickness (P"0.12 m), the drag force calculated by method I decreases sharply to zero, whereas that calculated by Autodyn2D decreases more moderately. This is due to bulging of the target's rear surface that is modeled by the Autodyn2D calculation and not in method I (for which the boundary condition is evaluated on the undeformed target's rear surface). Moreover, even if the boundary condition in method I was applied on the target's deformed rear surface, the results would not be signi"cantly di!erent because the assumed velocity "eld does not model this bulging phenomenon. These results indicate that although the target geometry is "xed in method I and the velocity "eld does not describe the exit stage correctly, it is possible to use the weighting function w to obtain
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
289
a reasonably good description of the drag force for this exit stage. Speci"cally, for this choice of the weighting function, the deviation (relative to Autodyn2D) of the residual velocity of the projectile is less than 0.5% (Fig. 10a). Also, it can be seen from Fig. 10b that method II produces the correct qualitative physical behavior of the drag force acting on the projectile by the target. This is due to the fact that method II satis"es the boundary and continuity conditions pointwise at the intersections of these boundaries with the streamlines. However, as was observed previously with the semi-in"nite target, method II overestimates the drag force acting on the projectile.
7. Conclusions Two di!erent theoretical approaches (method I [2] and method II [6]) were examined and compared with numerical calculations of Autodyn2D. Method I is based on satisfying the balance of linear momentum equation pointwise in the target, while approximating the boundary and continuity conditions only along the symmetry axis. On the other hand, method II satis"es the balance of linear momentum and the boundary and continuity conditions pointwise only along instantaneous streamlines relative to the target. For method I it was necessary to introduce approximate forms of the boundary and continuity conditions using appropriate weighting functions in order to obtain physically reasonable response during the exit phase of the projectile from a "nite thickness target. The physical implications of these approximations were investigated by comparing the results with numerical solutions of the exact formulation of the penetration problem using the hydrocode Autodyn2D. For the theoretical method I it follows that E The model produces the correct value of the drag force for penetration into a semi-in"nite target. E The solution for the exit stage of a "nite thickness target is sensitive to choice of the weighting function. E The behavior of the drag force applied to the projectile during the exit stage can be signi"cantly improved by a proper speci"cation of the weighting function. Moreover, for the theoretical method II it follows that E The model overpredicts the value of the drag force for penetration into a semi-in"nite target. E The model yields physically plausible response of the drag force during the exit stage of the penetration process since it satis"es the boundary and continuity conditions pointwise along streamlines. It is also important to emphasize that E The typical calculation time for a penetration process is reduced from several hours for Autodyn2D to only a few minutes for the analytical model. In summary, it can be concluded that both analytical methods I and II yield reliable results for integral parameters of the process such as the penetration depth, the ballistic limit and the residual velocity of the projectile. Moreover, it was shown in [2] for method I and in [6] for method II that these methods compare well with numerous experimental data. The results of the present work have shown that these methods are also capable of predicting such "ne characteristics as the shape
290
G. Yossifon et al. / International Journal of Impact Engineering 25 (2001) 265}290
of the elastic}plastic boundary and the velocity "eld in the target. Furthermore, the temporal evolution of the penetration velocity and the drag force predicted by these methods agree quite well with the data obtained from simulations of Autodyn2D. Also, during the period when the projectile penetrates the target's rear surface, the time dependence of the drag force can be predicted by the analytical model when use is made of the weighting function w"w which emphasizes the value of axial stress on the symmetry axis. Furthermore, it is noted that the results presented for the particular case in this work are characteristic of those obtained for other values of the parameters.
Acknowledgements The work of M.B. Rubin and A.L. Yarin was partially supported by the fund for the promotion of research at Technion.
References [1] Backman ME, Goldsmith W. The mechanics of penetration of projectiles into targets. Int J Engng Sci 1978;16:1}99. [2] Yarin AL, Rubin MB, Roisman IV. Penetration of a rigid projectile into an elastic}plastic target of "nite thickness. Int J Impact Engng 1995;16:801}31. [3] Hill R. Cavitation and the in#uence of headshape in attack of thick targets by non-deforming projectiles. J Mech Phys Solids 1980;28:249}63. [4] Tate A. A comment on a paper by Awerbuch and Bodner concerning the mechanics of plate perforation by a projectile. Int J Engng Sci 1979;17:341}4. [5] Ravid M, Bodner SR. Dynamic perforation of viscoplastic plates by rigid projectiles. Int J Engng Sci 1983;21:577}91. [6] Roisman IV, Yarin AL, Rubin MB. Oblique penetration of a rigid projectile into an elastic}plastic target. Int J Impact Engng 1997;19:769}95. [7] Roisman IV, Weber K, Yarin AL, Hohler V, Rubin MB. Oblique penetration of a rigid projectile into a thick elastic}plastic target: theory and experiment. Int J Impact Engng 1999;22:707}26. [8] Autodyn2D. Interactive non-linear dynamic analysis software, theory manual. Revision 4.0, Century Dynamics, 1998. [9] Yossifon G. Penetration of a rigid projectile into a metal multilayer target and characterization of the debris cloud. MSc thesis, Technion * Israel Institute of Technology, Haifa, Israel, May 1999. [10] Rubin MB, Yarin AL. On the relationship between phenomenological models for elastic}viscoplastic metals and polymeric liquids. J Non-Newton Fluid Mech 1993;50:79}88 and 1995;57:321.