CC BY
40
Computer simulation of target penetration on the basis of a combined discrete-continual approach
M.A. Chertov, A.Yu. Smolin, Yu.P. Stefanov, S.G. Psakhie, and Dewu Huan1
Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia 1 Shenyang Institute of Technology, Shenyang, 110015, China
In the present paper a combined discrete-continual approach proposed in the previous work [1] is applied to simulation of finite thickness steel plate penetration by long tungsten rod. Advantages of this approach for description of penetration, which go from the combination of best properties from a continual numerical method and the discrete numerical method of movable cellular automata (MCA), are discussed and demonstrated by computation of practical problem. In computations the influence of such parameters as impact velocities and plate thickness on penetration process pattern is studied. The dependences of final velocity and rod length as functions of impact velocity and plate thickness are calculated and verified against experimental data. The calculated curves and overall behavior of the simulated system have demonstrated sufficient compliance with experiment, which proves that the proposed approach is effective for description of penetration and other complex dynamic problems.
1. Introduction
It is no doubt that the penetration of target plates by high velocity impact is quite complicated physical phenomenon. One of the main reasons for the complexity of this process is that it simultaneously combines many different physical aspects such as plasticity, fracture, heating and heat transfer, plastification of target material due to heating, diffusion of penetrator and target material, etc. For the higher velocities (above 2.5-3 km/s) there is so-called hydrodynamic model that predicts the penetration depth of the impacting rod fairly well. This model is based on the assumption that at very high velocities spherical part of stress tensor or hydrodynamic pressure prevails components of stress deviator and, hence, the transversal stiffness of both the rod and target material can be neglected. At the same time, for the intermediate velocity range there are numerous phenomenological and dislocation models trying to describe the material response under dynamic loading. These models exhibit different properties concerning their precision at different strain rates, computational expenses required and to which extent underlying physical mechanisms of plasticity are accounted for. However, there is still a demand to build a model with reasonably small number of theory constants, which can adequately describe the gradual loosing of transverse stiffness by the material in the whole velocity range.
It is very important for the simulation of penetration process because the whole range of strain rates from zero to maximum at the channel is present in the specimen at the same time.
Computer simulation of high-rate deformation processes provides the invaluable possibility to have the insight of process dynamics with unlimited temporal resolution, which is impossible in the experiment; it considerably reduces the number of experiments as well as expenses allowing the extrapolation of experimental data. The existing grid-based computational codes using finite element or finite difference methods give numerical results [2, 3] that match reasonably well the experimental data. However, Lagrange grid codes used for simulation of solids have some well-known problems like high cell distortions, which demand complicated remeshing algorithms, resulting in accuracy losses. Many problems of grid codes has been overcome in meshless particle methods such as smoothed particle hydrodynamics (SPH) [4] and generalized particle algorithm (GPA) [5], etc. These methods, especially SPH, became very popular during the last decade.
One of the modern discrete methods used for simulating various aspects of mechanical behavior of heterogeneous materials and structures is the movable cellular automaton method (MCA) [6]. Unlike SPH, which is special approach
© M.A. Chertov, A.Yu. Smolin, Yu.P. Stefanov, S.G. Psakhie, and Dewu Huan, 2004
for discretization of differential equations of plastic flow, MCA is formulated in terms of particle interaction, each particle representing small but finite volume of material. The main advantage of the MCA method is its ability to naturally describe fracture of material up to its fragmentation and further mass mixing of fragments. The method has proved itself to be a very effective and prominent instrument for simulation of materials and structures under rather small elastoplastic deformations [6-8]. To enable computational study of high-rate processes such as shock loading, high-velocity impact, etc., the MCA-method has been improved by several modifications [9], description of them is also included here.
2. Modification of the MCA method for simulation of high-rate processes
Within the framework of the MCA method equations of motion for elements (automata) can be written as [6, 7]
J
dt2
i d^ dt2
= Fn +
E f\
=E K
(i)
where m1 is the mass of automaton 1; J1 is its moment of inertia; Ri is the position; e1 is the angle of the relative rotation; Fj = pj + O is the central component of the interaction force in the pair ij, t j is the tangential component of this force, FQ is the volume-dependent force acting on the automaton 1, Kj = qj (nj X O), qj is the distance from the center of the automaton 1 to the point of its contact with the neighbor j, nj is the unit vector defined as nj = = (Rj — R1 )/(qj' + qji).
This formulation explicitly accounts for the response caused by the volumetric changes. The volume-dependent force FQ acts due to the pressure from all the neighboring automata. The automaton pressure is determined by the change in the automation volume. In the simplest (linear) case this dependence takes the form
Pj = Kj £ j ,
(2)
where ej = (Qj -Q0)/Q0 is the bulk deformation; Q is the initial (equilibrium) volume of automaton j; Q1 is the current volume and KJ is the bulk modulus. The change in the automaton volume in a time At corresponds to computation based on the corresponding changes in the distances Aqj from the center of this automaton to the points of its contact with neighbors, i.e.
Nj
К Nj / ^
апj = njEAqj/ Eqjk Nj
k=0
k=0
(3)
Attempts to simulate high-rate processes such as shock loading, high-velocity impact, etc. using linear equation for pressure (2) has resulted in incorrect behavior due to large automaton overlapping. It is well known both from experimental data and from theory of interatomic interactions that response of the material is linear just only in the region of small deformations near equilibrium. When compression gets higher, interaction becomes essentially nonlinear because of non-linearity of interatomic potentials and, in case of shock loading, because of adiabatic heating of the material. To avoid non-physical automaton overlapping pressure in (2) must grow much faster than the linear dependence of bulk strain. It is well known that experimental data on thermodynamic properties of materials is often described by the shock adiabat, which in wide range of pressure can be represented by the linear relation [10]:
D = C0 + bU.
(4)
Here D is the shock-wave velocity and U is the mass velocity of the substance beyond the shock-wave front. Material constants c0 and b has been measured for a wide range of materials and can be easily found in the literature.
From (4) and from the Hugoniot relations for conservation of mass and momentum
e = — U/D, P = pDU, the following relation can be derived
P = Pc0
= K-
(5)
(6)
(1 + be)2 (1 + be)2
It is the simplest non-linear dependence of pressure on bulk strain. To use it in the MCA method it is natural to redefine the bulk elastic modulus K in (2) as
K ' = dP = K-1 - b£
d £
(1 + b£)3
(7)
where NJ is the number of neighbors of the automaton j.
Thus, as far as parameter c0 is assumed to be defined by the quasistatic bulk modulus K = pc0, the introduced non-linear equation of state requires only one additional parameter b that determines the inclination of D-U diagram. It can be found from the analysis of experimental shock adiabats or directly from the literature.
It is no doubt that for more accurate simulation of high rate deformations the introduction of nonlinear response function is not the only improvement required. It is also necessary to account for strain rate influence on the stress-strain diagram of the material. Heating due to adiabatic compression, which is now neglected, can have effect on mechanical properties as well. However, according to our understanding the introduction of nonlinear response function under high bulk compression is the most important step to be done. Note that the improvement discussed is also useful for simulation of quasistatic processes where high bulk compression occurs, for example, deformation under
£
£
Fig. 1. Rod and target plate structure at 0, 20 and 40 |xs
constrained boundary conditions, which is characteristic for processes in geological media.
The combined discrete-continual approach proposed in [1] allows to unite advantages provided by the grid and MCA methods. The MCA method is more suitable for description of intensive damage generation and mass mixing, while the grid method is several times faster than the particle method. So, it is natural to simulate the parts of structure, where deformations are small with the grid method to shorten calculation time and use the particle method in critical zones to increase space resolution. This approach makes easier calculation of full-size objects, without that it was necessary to simulate only fragments of the structure or increase the discrete element size. For simulation of penetration process the possibility to consider full-size object is quite important, because boundary conditions of the armor plate and properties of the fixers can have strong influence on the overall penetration pattern. It is also important that the alignment of grid with particles must be physically correct. To check this, in [1] authors considered propagation
of different wave types through the plane boundary. Although in [1] the grid method has been combined with the MCA method developed for small deformations, the proposed discrete-continual approach requires no special adjustment if the described MCA method for high-rate deformations is used, because deformations in the grid segment are still small.
3. Computation results
As an illustration of applicability of the combined MCA and grid method let us consider the results of simulation of steel plate penetration by long tungsten rod under normal impact in plane strain statement. Plate thickness is varied from 15 to 40 mm with step of 5 mm. The dimensions of the impacting rod are 10 X100 mm, impacting velocity is 1 520 m/s. Data on task geometry and loading conditions are taken to correspond to calculations described in the previous work [9]. The width of the area filled with particles is 80 mm, according to the calculation results such width is enough, because in this case the damaged zones don’t reach the interface boundary between grid and particles. Unlike [9], where special type of viscous boundary conditions is applied to the left and right side of the plate, here the plate ends are fixed.
Figure 1 shows initial, intermediate and final states of the simulated structure. One can see the intensive “mushrooming” of the impacting rod and formation of lips near the entrance of the penetration channel. The deformation and damage distribution presented here is not absolutely symmetric, because of stochastic distribution of strength properties of automata, which is added for more realistic fracture pattern.
Table 1 summarizes data on residual rod length and velocity after penetration in comparison with previous calculations [9] and other sources [2, 11]. The comparison indicates that rod deceleration is slightly decreased with regard to [9]. A possible reason for that is modified boundary conditions.
Table 1
Relative decrease of rod velocity and length as a function of plate thickness
Height, MCA + Grid MCA only
mm AU/U 0 Al/l0 AU/U 0 Al/l0
15 0.06 0.25 0.04 0.20
20 0.07 0.26 0.05 0.21
25 0.07 0.28 0.06 0.05* 0.02** 0.23 0.24* 0.19**
30 0.09 0.29 0.08 0.26
35 0.11 0.32 0.10 0.30
40 0.12 0.41 0.12 0.39
* — experimental data from [11]; ** — computer simulation results from [2]
4. Conclusion
The MCA method modified for description of high-rate deformation processes [9] in a combination with the continual grid method [12] is applied for simulation of thin plate penetration. The calculated dependences of residual velocity and erosion of the rod as functions of target plate thickness are compared with previous results obtained by MCA [9] and with results from other sources [2, 11]. It is shown that results of the combined discrete-continual approach are in good agreement with previous data. Overall penetration pattern seems quite realistic and there is no spurious behavior on the boundary between grid and particles. It can be concluded that the new approach is acceptable for simulation of penetration.
At present stage, the boundary between grid and particle segments can be only linear and have to be assigned before calculation. The main advantage of using mixed grid and particle approximation is the increased computation speed with the same accuracy. One of the goals for future development is dynamic transformation of highly distorted grid elements into particles. This would additionally increase calculation speed and make possible simulation of wider range of problems with more complicated and unknown in advance distribution of areas with small and intensive deformations.
References
[1] S.G. Psakhie, A.Yu. Smolin, Yu.P. Stefanov, P.V. Makarov, and M.A. Chertov, Simulation of complex medium behavior on the basis of combined discrete-continual approach, Pis’ma v ZhTF, 30, No. 17 (2004) 7.
[2] G.R. Johnson, R.A. Stryk, T.R. Holmquist et al., Int. J. Impact Engng., 10 (1990) 281.
[3] R.M. Summers, J.S. Peery, M.K. Wong et al., Int. J. Impact Engng., 20 (1997) 779.
[4] G.R. Johnson, R.A. Stryk, and S.R. Beissel, Comput. Meth. Appl. Mech. Engng., 139 (1996) 347.
[5] G.R. Johnson, S.R. Beissel, and R.A. Stryk, Comput. Mech., 25 (2000) 245.
[6] S.G. Psakhie, G. Ostermeyer, A.I. Dmitriev et al., Comput. Mater. Sci., 16, No. 4 (1999) 333.
[7] S.G. Psakhie, G. Ostermeyer, A.I. Dmitriev et al., Phys. Mesomech., 3, No. 2 (2000) 5.
[8] S.G. Psakhie, Y. Horie, G.P. Ostermeyer et al., Theor. App. Fract. Mech., 37 (2001) 311.
[9] M.A. Chertov, A.Yu. Smolin, E.V. Shilko et al., Simulation of Complex Plane Targets Penetration by Long Rod Using MCA Method, in Proc. of 6th Int. Conf. for Mesomechanics “Multiscaling in Applied Science and Emerging Technology”, Greece, Patras (2004) 307.
[10] W.H. Gust, J. Appl. Phys., 53, No. 5 (1982) 3566.
[11] V.M. Fomin, A.I. Gulidov, G.A. Sapozhnikov et al., High-velocity interaction of solids, Izdatelstvo SO RAN, Novosibirsk, 1999.
[12] M.L. Wilkins, Calculation of Elastic-Plastic Flow, in Methods in Computational Physics, Vol. 3, ed. by B. Alder, S. Fernbach and M. Rotenberg, Academic Press, New York (1964) 211.
