1 Introduction

Modeling of normal and oblique penetration is of great importance in the defense industry which remains as an active research field currently. Repeated simulations and tests are common in the design and protection of important buildings where a fast and effective computational technique is urgently needed. Over the past several decades many meaningful research results have been done^[1-6] which build the foundation of research field on penetration. The review article by Ben-Dor et al.^[7] provides a detailed survey which summarizes the past and present analytical work on the penetration problem.

The method for penetration research includes experimental study, theoretical analysis, and numerical investigation. Engineering models including empirical and semi-empirical formulae are derived from the penetration test. Meaningful work with penetration by Forrestal et al.^[8-9] using a semi-empirical formula to predict the depth of penetration has been proven valid and effective for the normal penetration with a rigid projectile. Frew et al.^[10] have presented a semi-empirical formula to predict the depth of penetration into limestone target and obtained the empirical formula to calculate the resistance of the limestone target. For the oblique penetration, Warren and Poormon^[11] and Warren et al.^[12] have conducted a series of scaled tests on the penetration into aluminum and limestone targets using steel projectiles. Based on experimental data, Charles et al.^[13] have derived a penetration model for metallic targets with a least-squares regression analysis method. In general, engineering models and test data can be achieved through scaled tests. However, the tests are costing and inconvenient to conduct at the same time.

For theoretical analysis of penetration, the method of cavity expansion analysis (CEA) firstly proposed by Bishop et al.^[1] is widely used to derive the force between the projectile and target. Forrestal et al.^[14-15], Forrestal and Luk^[16], and Forrestal and Tzou^[17] have studied the cylindrical and spherical CEA in brittle materials. Zhen et al.^[18] have developed an analytical solution for the dynamic cylindrical cavity expansion with a finite radius considering the effect of the lateral free boundary. Besides, Chen and Li^[19] and Chen et al.^[20-21] have proposed general non-dimensional formulae for penetration and perforation of rigid projectiles which are consistent with experimental data. Generally, theoretical analysis is not adequate for complex oblique penetration problems which are commonly studied with experimental and simulation method.

Numerical investigation has been proven economical compared with scaled tests. Currently, the commercial finite element software is widely used for the simulation of penetration process. Considering the separation-reattachment behavior, Bernard and Craighton^[22] have developed a two-dimensional finite difference code for penetration of rigid projectile. Warren and Poormon^[11] and Warren et al.^[12] have developed an explicit transient dynamic finite element code PRONTO 3D based on a combined analytical and computational technique. The code is effective in the simulation of oblique penetration with deformable projectiles. Recently, Fang and Zhang^[23] have proposed a three-dimensional finite element analysis approach to simulate the penetration into rock-rubble overlays. The traditional finite element analysis is time-consuming especially for the deep penetration problems where the fast and repeated simulation cannot be realized. Besides, researchers mainly focus on the normal and two-dimensional penetration which is too simple for the real penetration. The three-dimensional gesture of projectile and penetration trajectory should be considered especially for the penetration into targets with randomly distributed stone and steel. Thus, a fast algorithm for three-dimensional penetration is needed for the current research and design activity.

In this paper, a fast algorithm is proposed for three-dimensional penetration of a rigid projectile into limestone target. The algorithm is designed on the basis of the isolation between the projectile and target. Thus, the analytical load on the projectile is used instead of the force between the projectile and target. The calculating time is saved by ignoring the response in the target which can realize the rapid prediction of penetration trajectory. The corresponding influence factors including the analytical load model, the cratering effect, the free surface effect, and the separation-reattachment phenomenon are considered in the algorithm. Besides, the method of the cavity ring is used to study the cavity expansion process. Furthermore, the description of the projectile's three-dimensional gesture is presented in this paper. As a result, the fast algorithm is programmed into a fast simulation code of three-dimensional penetration, named PENE3D. Besides, a series of selected normal and oblique penetration experiments are calculated using the PENE3D. The results are compared with the test data and the simulations in Ref. [10] and Refs. [11]-[12].

2 Design of fast algorithm

The theoretical analysis of penetration is mostly based on the CEA. For more specific analysis, the free surface effect and the separation-reattachment behavior should be considered to which the cavity expansion process is closely related. Besides, the cratering effect should be taken into account for the penetration into the brittle target like concrete and rock. With reference to the theory of ballistic trajectory, we propose a method to describe the three-dimensional gesture of the projectile which is common in the real penetration.

2.1 CEA 2.1.1 Analytical model of loading

The loading force on the projectile can be divided along two directions of the projectile's surface, the tangent direction and the normal one, as shown in Fig. 1. σ_n is used to describe the normal force on unit area^[8] as

(1)

where R_t is the static resistance of the concrete, and R_t is defined as an target strength constant related to the resistance to expand the cavity in the target. ρ_t is the density of the concrete, v_n is the normal velocity of the segment on projectile, and n is the normal unit vector of the segment on the projectile.

The static resistance R_t for limestone is related to the shank diameter of the projectile and can be calculated with the formula^[10] as

(2)

where K and k are coefficients achieved from data fitting. 2a₀ is the reference projectile diameter, and 2r_p is the diameter of the projectile. With K=607 MPa, k=86 MPa, and 2a₀=25.4 mm, the measured static resistance R_t is recovered for each data set in Frew's paper.

The friction force of the projectile in the penetration is calculated as

(3)

where τ is the unit vector in the tangent direction, and c_f is the coefficient of friction which is set as 0.01 in this paper.

2.1.2 Cavity expansion process

The target is considered as a uniform fluid medium for the analysis of the cavity expansion process. As shown in the diagram of Fig. 2, the origin of the cylinder coordinate system is set as the initial striking position of the penetrator. Then the target is divided into several discrete cavity segments with the width dz. Based on the conservation of mass and the assumption that the target is incompressible, the radial expansion velocity u of the target segment is expressed as the function of the cavity radius a_c, the cavity expansion velocity on the surface of the cavity that

(4)

where r is the distance from the segment to the axis of the cavity.

We assume that the work dE_s done by the projectile on the target segment can be totally transformed into the kinetics dE_k and the deformation energy dE_p of the target segment after the separation of the projectile and target. Then dE_k and dE_p are calculated by the formulae that

(5)

(6)

where R is the distance from the cavity axis to the boundary of the plastic region of the target, and P_t is the work needed to form the cavity of unit volume. For the energy conservation of the target segment, we can get the relationship of variables at time t_n and t_n+1 that

(7)

So we can achieve the cavity expanding velocity and the cavity radius a_{c, n+1} from and a_{c, n} using formulae (6) and (7) . The largest cavity radius is obtained if the cavity expanding velocity decreases to zero.

2.1.3 Target's free surface effect

In the oblique penetration, the assumption of infinite medium space is inappropriate for theoretical derivation of the loading model considering the influence from the free surface of the target. Thus, Warren and Poormon^[11] and Warren et al.^[12] proposed a decay function to fix the analytical loading. The distance from the axis of the projectile to the free surface of target d and the radius of the head of the projectile in the cross section perpendicular to the axis a are used as parameters in the derivation. In this paper, we use the same formula to analyze the free surface effect function that

(8)

(9)

where E is the Young's modulus of the target, τ is the shear strength of target, λ is a coefficient related to pressure, b=a_f(2E/(3τ))^1/3 is the radius of the plastic region in spherical cavity expansion theory, and α =6/(3+2λ) is an intermediate variable for derivation.

2.1.4 Cratering effect

A conical crater is always formed during the penetration into the brittle target, thus the penetration trajectory is divided into the crater region and the tunnel region^[10]. He^[25] considered an increasing stress on the projectile which is related to the depth of penetration. Thus, based on the same assumption with He^[25], a crater effect function is proposed here that

(10)

where h is the depth of penetration, and r_p is the radius of projectile.

The depth of crater is chosen as 5 times of the projectile's diameter which is in accord with the results of scaled test. The cratering effect can be neglected if the depth of crater is so little that the free surface effect plays a dominant role in the calculation of analytical model of loading.

2.1.5 Separation-reattachment phenomenon

The separation-reattachment-separation phenomenon always exists between the projectile and target at the local position. Based on Benard's work^[22], a theoretical analysis on the separation-reattachment phenomenon is put forward for three-dimensional penetration.

The approach angle is defined as the angle between the projectile and the cavity of the target, and can be expressed by the velocity of the mass center of projectile V and the velocity of the certain point on the projectile in the outer normal direction v_n.

(11)

While the certain point on the projectile separates from the inner surface of the tunnel, we assume that the approach angle for this point is equal to a certain angle which is defined as the separation angle a_s. For normal penetration with no yaw, the separation angle for different azimuth angles is the same, so the separation points are on the same circumference. For oblique penetration, the circular cross-section radius on the separation point is related to the azimuth angle. The projectile separates from the cavity of target if the approach angle is less than the separation angle.

For the reattachment phenomenon, the projectile reattaches the cavity if the distance from the axis to the surface of the projectile r_a is greater than the radius of the cavity r_c. The variables are expressed by the formulae

(12)

(13)

(14)

where r_o is the circular cross-section radius, ζ is the distance from the point on the axis of the projectile to the tip of the projectile, ζ_o is the distance from the separation point to the tip of the projectile, and δ is the lateral displacement from the certain point on the axis of the projectile to the axis of the cavity. ω_x and ω_y are the angular velocities of the projectile. V_x, V_y, and V_z are the translational velocities of the projectile in the different axis directions.

2.2 Analysis of projectile's three-dimensional movement 2.2.1 Description of projectile's three-dimensional gesture

The gesture of projectile is determined by the angles between the vectors such as the vector of the velocity, the unit vector of the projectile's axis and the normal vector of the target's surface. For the approximate two-dimensional penetration, oblique angle, pitch, and yaw are often used to describe the gesture of projectile which is inconvenient for the three-dimensional description. The oblique angle θ_y, θ_x and the yaw α_x, α_y are used as the variables to express the projectile's three-dimensional gesture in this paper, as shown in Fig. 3. The relation between these two kinds of description needs to be derived for the usage of the existing penetration results in the calculation with the fast algorithm in PENE3D.

As in Fig. 3, the description with an oblique angle, a pitch angle, and a yaw angle is expressed by two oblique angles and two yaw angles. If the downward pitch angle and the right yaw angle are defined as positive, the relationship of these two descriptions is expressed as follows:

(15)

where θ is the oblique angle, θ_P is the pitch angle, α_Y is the yaw angle, θ_y is the oblique angle in the OXZ plane, θ_x is the oblique angle in the OYZ plane, α_x is the yaw in the OXZ plane, and α_y is the yaw in the OYZ plane.

The local coordinate system which follows the projectile's movement is set with its origin at the mass center of the projectile. The z-axis is in the direction from the center of bottom to tip of the projectile, the x-axis and y-axis are determined by two points on the projectiles according to the projectile's initial gesture. When the projectile penetrates normally, the local coordinate system is the same as the global one at the impact moment. The relation of the local coordinate and the global coordinate is that

(16)

where M_t is the transformation matrix, and P_Global and P_Local are the same vector described in each coordinate system, respectively. The base vectors of the local coordinate system e_x^l, e_y^l, and e^zl are expressed with the global coordinates such as a_x, b_x, c_x that

(17)

Then the transformation matrix M_t is obtained by the base vectors that

(18)

2.2.2 Finite difference scheme for calculation

The movement of the projectile can be divided into two parts: the translational movement of the mass center of the projectile and the rotation about the instantaneous axis which passes through the mass center. The finite difference scheme is used to calculate the variables of movement such as the displacement and the angle. The forward difference scheme is used in the first two calculating time interval, and the central difference scheme is used in the following time steps. The computational precisions for these two difference schemes are first order and second order, respectively. Thus, the central difference scheme will be mainly discussed here.

The movement of the projectile is analyzed with the components in each axial direction of the Cartesian coordinate system for simplicity. The equation of projectile's translational movement is expressed as

(19)

where F_x, F_y, and F_z are the components of resultant force exerted on the projectile; X, Y, and Z are the global translational displacements of the projectile. At time t_n+2 the axial components of displacement X_n+2, Y_n+2, and Z_n+2 are expressed as

(20)

where F_{x, n+1}, F_{y, n+1}, and F_{z, n+1} are the axial components of resultant force at time t_n+1, dt is the time interval, and m is the mass of the projectile.

Three assumptions must be made before analyzing the rotation of projectiles. (i) The projectile is rigid during the rotation. (ii) The rotation of the projectile about its own axis is neglected for simplicity. (iii) The resultant force and torque of the projectile keep constant during the calculating time step.

The angular acceleration of the projectile is calculated with the moment of inertia and the resultant torque by the theorem of angular momentum that

(21)

where I_axis and I_tran are the axial and translational moments of inertia, respectively. L_x, L_y, and L_z are the axial components of the resultant torque. , and are the axial components of the angular acceleration. According to the second assumption, the angular acceleration about the z-axis is zero during the penetration. Then formula (21) is simplified to the uniform vector equation.

(22)

where is the angular acceleration vector, and L is the resultant torque vector.

The central difference scheme is applied to the angular acceleration that

(23)

where and are the angular velocities at time t_n and time t_n+2, and is the angular acceleration at time t_n+1. According to the third assumption, the instantaneous axis of rotation is defined using the average angular velocity as

(24)

The unit vector in the instantaneous axis n_axis is in the direction of the average angular velocity as follows:

(25)

where is the module of .

The rotation of one vector around another vector is analyzed to obtain the new transformation matrix M_T. As in Fig. 4, vector p rotates about the unit vector n with the angle θ to get the new vector p'. The relation of these two vectors is expressed as

(26)

where M is the transformation matrix for the rotation of a vector.

The derivation of M is described in three steps. Firstly, p and n are moved to p₁ and n₁ through the rotation with their relative positions fixed, and the transition vector n₁ is parallel to the z-axis in the positive direction. Secondly, p₁ rotates about the unit vector n₁ with the angle θ to get the new vector p'₁. Finally, the inverse operation is taken to move p₁, n₁ and p'₁ back to the original position to get the vector p'. In summary, the matrix M is obtained by the formula

(27)

where x_u, y_u, and z_u are the coordinates of the unit vector n.

3 Simulation results and discussion

In this section, a series of simulations are conducted with the fast algorithm in PENE3D on the normal and oblique penetration, respectively. The simulation results are compared with the results of simulations and tests in the published papers^[10-12]. The diameter and length of the steel projectile are 7.11 mm and 71.1 mm, respectively; the target is made of limestone with an average static resistance of 913 MPa.

3.1 Simulation of normal penetration

For the normal penetration, the pitch and yaw can be neglected since they have little influence on the depth of penetration. The tests in Ref. ^[10] are simulated with the algorithm in PENE3D. In the simulations, the parameters of the projectile and target are listed in Table 1 and Table 2, respectively. The projectiles are modeled as VAR 4340 (R_c=44.5) steel with three CRH ogive-noses, a diameter of 7.11 mm, a length-to-diameter ratio of 10, and a quasi-static yield stress of 1 462 MPa. The targets are limestones with a nominal unconfined compressive strength of 60 MPa, a density of 2 300 kg/m³, a shear strength of 40 MPa, and a static resistance of 913 MPa. The parameters of the projectile and target are listed in Table 1 and Table 2 in which λ and τ are the parameters used to define the pressure-dependent shear strength in the Mohr-Coulomb yield criterion.

The comparisons between the simulations and the test results are presented in Table 3 and Fig. 5. The constant static resistance R=913 MPa and the static resistance derived from the test are both used as parameters in our simulation. The penetration depth is greater than the results of the test when using the static resistance derived from the test in the simulation. For the static resistance R=913 MPa, the penetration depth is greater than the results of the test except for the cases with the initial penetrating velocity of 1 230 m/s and 1 340 m/s.

The error of penetration depth is probably caused by the consideration of the cratering effect in the simulation. When compared with the test results, the percentage error of the penetration depth in the simulation is less than 14% with the striking velocity greater than 787m/s. The comparison of the results shows that the simulation is acceptable in engineering for the penetration with the velocity greater than 600 m/s. Besides, these two static resistances are both acceptable for the simulation.

3.2 Simulation of two-dimensional penetration

Typical results from simulations are presented and compared with the results from experiments and simulations in Ref. [12]. Table 4 summarizes the initial conditions and the final projectile positions for both experiments and simulations. The pitch and yaw are equivalently changed into the oblique angle and yaw defined in simulations which are convenient for data input. The results are obtained from the impact point to the final position of the nose tip and the bottom center of the projectile.

3.2.1 30° angle of obliquity

As shown in Fig. 6 to Fig. 11, the post-test cast from the experiment is compared with the final position of the projectile obtained from simulations at several typical velocities. The unit for displacement is mm in the simulation. The final depth of the projectile in each case is greater than that in the test. The reason is the upper change of trajectory caused by the bending of the projectile. The simulation is more accurate on the horizontal direction than that in the vertical direction.

In the simulations for projectiles striking the target with velocities of 413 m/s and 510 m/s, the deformation of the projectile seems to happen easily because the free surface effect and the yaw have strong influence. The depths of penetration in the simulation increase and the oblique angles of the projectile at the end of penetration decrease when compared with the results from tests. The main reason for the difference between these two results is the deformation of the projectile in the test. The deformation is due to the free surface effect of the target, and it is not considered in the algorithm of PENE3D.

In the simulations for projectiles striking the target with middle and high velocities such as 610 m/s, 804 m/s, 863 m/s, and 1 023 m/s, the gestures of projectiles at the end of penetration in the simulations are close to those in the test. While the depth of penetration in the simulations is a little greater than that from the experimental results. The reason for the differences is the crater effect considered in the simulation and the slight deformation which makes the projectile more easily to rotate in the test.

3.2.2 15° angle of obliquity

As shown in Fig. 12 to Fig. 15, the tests in Warren's paper are simulated for the projectiles with the striking velocity of 411 m/s, 510 m/s, 773 m/s, 940 m/s, and 1 030 m/s, respectively. The unit for displacement is mm in the simulation.

It is observed that the simulations with the fast algorithm in PENE3D give similar results with those in Warren's paper. The depth of penetration in the simulation is slightly greater than that in the test due to the rigid property of the projectile. The displacements in our simulation in the horizontal direction are close to the results of tests and the simulations using PRONTO 3D. Compared with the penetration with the 30° angle of obliquity, the projectile which strikes the target with the 15° angle of obliquity tends to keep rigid after penetration. Besides, the trajectory is straighter than that in the penetration with the 30° angle of obliquity.

It is also observed that the simulation results with a smaller separation angle are closer to the test results. The reason is that the separation angle is small in the test due to the relatively low penetration velocity of the projectile, which is proven correct according to the post-test radiograph. In conclusion, the fast algorithm in PENE3D does well in the simulation of the two-dimensional penetration with the small oblique angle. Besides, it is necessary to consider the deformation of the projectile during the penetration to obtain a more accurate simulation.

3.3 Discussion of three-dimensional penetration and error analysis

In the simulation according to Table 4, the displacements in the x- and z-directions are considered for validation. Besides, the displacements in the y-direction (in the yaw direction) are also obtained. The final projectile position in OYZ plane is shown in Fig. 16 with the velocity of 1 249 m/s and the oblique angle of 30°. The displacement in the y-axial direction is 14.8 mm. For the final projectile position shown in Fig. 17, the initial velocity is 1 255 m/s, the oblique angle is 15°, and the displacement in the y-axial direction is 12.6 mm.

In these two situations, the displacements in the y-direction cannot be ignored, which means that it is inappropriate to consider those tests as two-dimensional. A three-dimensional analysis needs to be done to describe the three-dimensional gesture of the penetrator properly. Besides, the three-dimensional fast algorithm in this paper can be helpful and effective for the prediction of three-dimensional penetration.

In real penetration tests, the projectiles always strike the target with a certain yaw which makes the penetration process three-dimensional. Besides, three-dimensional penetration simulation is necessary to predict the trajectory reasonably for the targets with steel reinforcement or with the randomly distributed stones. In these cases, three-dimensional simulations using commercial finite element software usually need several hours on the workstation for calculation, and these are time-consuming. While for the same case, the time consumed on the simulation using the fast algorithm PENE3D is just several seconds. Within the acceptable calculation error, the fast algorithm here is more effective and convenient than the commercial software for the current research and design activity.

Though the algorithm is efficient to predict the penetration trajectory, shortages still exist in the fast algorithm for the simplified assumptions made in the theoretical analysis. Besides, systematical errors are existing in the finite difference scheme and the programming work. Furthermore, since the rigid assumption is made, errors also exist for the cases where the projectile transforms obviously. Thus, several developments are needed to decrease the errors. Firstly, the constitutive law which considers viscosity can be considered in the algorithm to obtain a more accurate loading model. Secondly, the finite element method can be added to the proposed algorithm to simulate the deformation of the projectile. Further work is still needed to obtain a more accurate description for the separation-reattachment phenomenon and the movement of the projectile with rotation about its own axis. The factors such as the mass loss and the deformation of the projectile should also be considered in our fast algorithm in the future research.

4 Conclusions

In this paper, a fast algorithm is proposed for three-dimensional penetration of rigid projectile into limestone target. The algorithm is designed on the basis of the isolation between the projectile and target. Thus, the analytical load on the projectile is used instead of the force between the projectile and target. The calculating time is saved by ignoring the response in the target which can realize the rapid prediction of penetration trajectory. The corresponding influence factors including the analytical load model, the cratering effect, the free surface effect, and the separation-reattachment phenomenon are considered in the algorithm. Besides, the method of the cavity ring is used to study the cavity expansion process. Furthermore, the description of the projectile's three-dimensional gesture is presented in this paper. As a result, the fast algorithm is programmed into a fast simulation code of three-dimensional penetration, named PENE3D.

A series of selected normal and oblique penetration experiments are calculated by the PENE3D. The results are compared with the test data and the simulations in the papers of Frew et al.^[10] and Warren and Poormon^[11] and Warren et al.^[12]. The predicted results are approximately in accord with the results from tests. The percentage error of the penetration depth is less than 16.2% with the striking velocity greater than 597 m/s in the simulation for normal penetration. Besides, the two-dimensional simulation of oblique penetration is adequate for the penetration with rigid projectiles. The simulation is accurate for the horizontal displacement and the oblique angle of the projectile. The general increase of the penetration depth could be caused by the crater effect and separation angle that considered in the simulation. Furthermore, three-dimensional penetration should be considered based on the simulations of the selected oblique tests.

In conclusion, the proposed algorithm in PENE3D in this paper is fast and effective in the simulation of the penetration process and the prediction of penetration trajectory. The algorithm can be easily integrated into a broader software which is of great help for the design of penetration weapons and the estimation before the penetration test. While errors still exist in the fast algorithm for the simplified theoretical assumptions such as the rigid assumption. Besides, systematical errors also exist in the finite difference scheme and the programming work. Thus, the finite element method and a better constitution law can be added to the proposed algorithm to decrease the errors. The factors such as the mass loss and the deformation of the projectile should also be considered in our fast algorithm in future research.