- Email: [email protected]
Prediction of fragmentation and experimentally inaccessible material properties of steel using finite element analysis
Finite Elements in Analysis and Design 104 (2015) 72–79
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
Prediction of fragmentation and experimentally inaccessible material properties of steel using finite element analysis Jeremy M. Schreiber a,b,n, Ivi Smid a, Timothy J. Eden b, David Jann c a
Department of Engineering Science and Mechanics, The Pennsylvania State University, University Park, PA 16802, USA Applied Research Laboratory, The Pennsylvania State University, University Park, PA 16802, USA c DE Technologies, King of Prussia, PA 19406, USA b
art ic l e i nf o
a b s t r a c t
Article history: Received 25 February 2015 Received in revised form 21 May 2015 Accepted 2 June 2015 Available online 26 June 2015
High strain-rate properties of materials are needed for predicting material behavior in extreme environments. The demand for high strain-rate properties continues to increase for commercial and military applications as the operating environments become more extreme, such as fragmentation, impact and explosions. To reduce time and expense, Finite Element Analysis (FEA) is being used to simulate these behaviors and reduce the number of experiments needed to characterize how a material performs at high-strain-rates. A finite element model for predicting fragmentation behavior of a high strength steel ring was developed using Abaqus Computer Aided Engineering (Abaqus) software. AISI 4340 steel, a low alloy Cr–Ni–Mo steel, was used in the analysis. The results of the finite element model were compared to the results from CTH, a two-dimensional Eulerian shock physics hydro-code. CTH was also used to develop a transient loading curve for the Abaqus model. The fracture strain in the model was adjusted to induce failure in the ring. Element deletion was used to model failure. A fracture strain less than 1 10 5 was needed to initiate fragmentation. The effects of mesh type and model defects were also investigated. & 2015 Elsevier B.V. All rights reserved.
Keywords: Fragmentation Finite element analysis Johnson–Cook High strain-rate Lagrangian analysis Eulerian analysis
1. Introduction Material fragmentation at high strain-rates is difficult to characterize experimentally. Typically, a test article or projectile is subjected to a high-strain-rate that results in fragmentation. The fragments are collected and sorted into bins. Several fragments are examined under a microscope to determine failure characteristics, making this process very time consuming and expensive. In the past 60 years, many empirical material models have been developed to predict failure at high strain-rates to correlate with experimental results [1–4]. Most of these models rely on high strain-rate data obtained using the Split Hopkinson Pressure Bar (SHPB) or other high strain-rate testing. The SHPB method measures stress pulse propagation through a metal bar to predict the stress–strain relationships of a material [5,6]. There are shortcomings to this method as it requires many assumptions to be made, and loses accuracy when materials undergo tension testing at strain-rates
n Corresponding author at: The Pennsylvania State University, 3075 Research Drive, State College, PA 16801, USA. Tel.: þ1 814 865 1096. E-mail address: [email protected] (J.M. Schreiber).
http://dx.doi.org/10.1016/j.finel.2015.06.001 0168-874X/& 2015 Elsevier B.V. All rights reserved.
above 103 s 1 [7,8]. Components can be designed to improve fragmentation behavior and identify failure initiation sites. The models help identify failure modes, and reduce the number of samples that need to be tested. For this effort, a model was developed that incorporates high strain-rate data from previous SHPB experiments in a constitutive material model to predict failure.
2. Experimental procedure and approach 2.1. Ring development A three dimensional solid model of a ring was developed using SolidWorks computer aided drafting (CAD) software [9]. The ring has a diameter of 81 mm, a wall thickness of 7 mm, and a notch depth of 3.5 mm. The notch is 5.7 mm tall with a 60 degree taper at the top. The ring was imported into Abaqus for analysis. The ring was partitioned into smaller sections to allow for the use of the automatic meshing function in Abaqus which is needed to utilize all mesh types when modeling high-strain-rate behavior. Partitioning in Abaqus was done by the three point method to develop a partitioning plane. The SolidWorks model of the ring is shown in Fig. 1.
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
73
Fig. 1. Three Dimensional SolidWorks model of the ring. The cross section view is shown on the right. The diameter of the ring is 81 mm. Wall thickness is 7 mm.
2.2. Material properties
Table 1 Strain-rate independent material properties of AISI 4340 steel [1,2,10].
The material properties used in the model of the ring were those of austenitized, quenched, and tempered AISI 4340 steel determined by Johnson and Cook [1,2]. The hardness of this steel was reported to be Rockwell C30, which would correlate to a mainly pearlitic microstructure. The standard, strain-rate independent material properties for this alloy are found in Table 1 [1,2,10]. High strain-rate properties were added to the model using the Johnson–Cook (J–C) plasticity model [1,2] This empirical model contains five experimentally determined constants and is commonly used to predict behavior of ductile materials under high strain-rate conditions such as explosive loading. This model is available in the Abaqus software by default [1,2,11]. It defines the von Mises flow stress of a material as a function of the power law relationship of plastic strain and strain hardening, strain-rate, and thermal softening. The J–C model is shown in Eq. (1) [1,2]. h ih i σ ¼ A þBεnp 1 þ Clnεnp ½1 T m ð1Þ H
Density (kg/m3)
200
7830
Poisson’s % ratio Elongation 0.29
Ultimate tensile strength (MPa)
Yield strength (MPa)
745
470
22
Table 2 Johnson–Cook material constants for AISI 4340 steel [1,2]. A (MPa)
B (MPa)
C
n
m
Tmelt (K)
800
510
0.014
0.26
1.03
1793
Predicted von Mises flow stress vs. true strain
1500
εp ¼ normalized effective plastic strain rate n
ðT T room Þ ðT melt T room Þ
The values of A, B, C, n, and m are experimental constants that are determined using uniaxial tension tests and the SHPB tests. Tmelt and Troom are the melting temperature of the alloy and the ambient air temperature, respectively. The material constants used in the J–C model and are found in Table 2 [1,2]. The J–C equation was used to predict the room temperature von Mises flow stress at six different strain-rates shown in Fig. 2. The results indicate that the von Mises flow stress increases as the strain-rate increases. Strain is only plotted to 0.4 strain due to the fact that the Johnson–Cook model only predicts a linear increase in strength at large amounts of strain. This linear increase at large strains is due to the effect of the exponential power law factor in the equation. The J–C strength model was not the only constitutive material model needed in this analysis. The J–C fracture model was used to initiate the failure of the elements through the use of a fracture strain value. Fracture strain can be measured easily under uniaxial tension at quasi-static strain rates using an extensometer, but at high-strain-rates this may not be possible. Modeling a complex structure at various fracture strain values maybe the only way to identify the proper fracture strain. In modeling, once an element reaches the fracture strain value, it is deleted from the equation. The J–C fracture model does not follow a typical nucleation and
von Mises Flow Stress (MPa)
1400
εp ¼ equivalent plastic strain
T H ¼ homologous temperature ¼
Elastic modulus (GPa)
1300 1200 1100 1000 900 800 Selected Strain Rates
700 600
0
0.05
0.1
0.15
0.2
0.00001
0.0001
1
100
0.25
0.3
0.01 10000
0.35
0.4
True Strain (mm/mm)
Fig. 2. von Mises Flow Stress vs. True Strain of AISI 4340 Steel at Room Temperature. Data is plotted at Room Temperature (300 K).
growth model. It is dependent on strain, strain-rate, temperature, and pressure [1,2]. The benefits are that it is less complicated to use, and most of the parameters can be found using adjusted quasi-steady state data. The J–C fracture model is also built directly into Abaqus, and is fully compatible with the J–C strength model. The J–C fracture model is defined in Eq. (2) [2], h i n εf ¼ D1 þ D2 expD3 nσ 1 þ D4 ln ε_ n ½1 þ D5 T n ð2Þ where εf is the equivalent strain to fracture, D1 –D5 are material constants, ε_ n is the dimensionless strain-rate, and T is temperature. Material constants for AISI 4340 from Johnson-and Cook are shown in Table 3 [2].
74
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
2.3. Model setup The first step was to start with an axisymmetric model of the ring to determine how large of a displacement was needed to initiate failure based on fracture strain. The fracture strain was introduced into the finite element solver using the J–C fracture model [11]. The values of fracture strain were arbitrarily chosen from very small strains to large strains, which necessitates multiple model runs. These values were input into the damage evolution option built into the Abaqus software. This option requires a displacement at failure input. This was set to the predicted fracture strain. Axisymmetric modeling allows for fast calculation times to determine rough values for overall von Mises stress, displacement and strain-rate before implementing a more calculation intensive 3D model. A dynamic, explicit analysis with an analysis time of 60 μs was chosen for all models to cover the entire high-strainrate event. An instantaneous pressure of 1.2 GPa was applied to the inner side of the ring to simulate explosive loading conditions in the axisymmetric models and initial 3D modeling [12,13]. This value was chosen to give both a baseline value for evaluation, but also a maximum pressure load found in some high explosives. The model was constrained on top and bottom faces of the ring to prevent model runaway. Runaway is a condition where the model is under-constrained and does not perform as predicted. This boundary condition was propagated to all subsequent models. Once the baseline behavior was determined using axisymmetric modeling, a 3D model was developed. A time dependent loading curve was investigated using 3D modeling only.
2.3.1. Eulerian and Lagrangian analysis Two different methods that are commonly used to analyze dynamic material systems are the Lagrangian and Eulerian [14]. The difference between these two methods is that in a Lagrangian method, the mesh is associated with a certain section of the part, whereas in the Eulerian method, the mesh is associated with a larger computational domain and the part is not assigned to the mesh, but moves freely throughout the meshed domain. The advantage of using a Lagrangian analysis is that the mesh is fixed to certain coordinates on the part, making adding boundary conditions and contact parameters easy [15]. In contrast, in the Eulerian analyses, the material is constantly moving through the computational domain of the system, leading to greater difficulty in assigning boundary and contact parameters [15]. Eulerian models excel when high strain is expected due to the retention of the mesh geometry. Mesh distortion is a concern in Lagrangian finite element models and there are re-meshing techniques that are typically used to address this issue. The downside of having a fixed mesh in an Eulerian analysis is that the overall model domain must be much larger to prevent material from leaving the computational domain, which increases calculation time [15]. In this work, the Eulerian CTH code is used to determine the loading profile for the Lagrangian Abaqus 3D model. The deformation behavior is also compared between the two models.
3. Results and discussion 3.1. Element overview
Table 3 Fracture parameters for AISI 4340 steel given by Johnson–Cook [1,2]. D1
D2
D3
D4
D5
0.05
3.44
2.12
0.002
0.61
The element types used in finite element modeling are chosen to provide specific information about a model. Each element has different properties, which give it certain characteristics for modeling such as hourglass control. This modeling effort utilized
Table 4 Summary of element types used in each modeling step. Geometry
Element
Abaqus designation
Settings
Comments
Axisymmetric
4-Node quadrilateral, linear
CAX4R
Stress/displacement element without twist
3D
8-Node brick, linear
C3D8R
3D
10-Node tetrahedron, quadratic 4-Node tetrahedron, linear
C3D10M
Reduced integration, hourglass control Reduced integration, hourglass control Modified element, hourglass control Default
3D
C3D4
Stress/displacement element Stress/displacement element, three additional displacement variables Stress/displacement element
Fig. 3. Predicted failure initiation values for uniform loading curve at t¼ 45 μs. Degree of failure is a qualitative estimate of where the model does not show any change in failure characteristics.
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
75
Fig. 4. Meshed 3D models before analysis. The left is meshed with an 8 node linear brick element, and the right side is meshed with a 10 node tetrahedral element.
Fig. 5. Results of both brick and tetrahedral models after 45 μs. The tetrahedral model gives a larger maximum von Mises stress in the surface regions. Stress units are in Pascal.
four different types of solid continuum elements. Table 4 gives an overview of the elements used in the modeling. 3.2. Axisymmetric modeling Analyses were performed using the axisymmetric models with fracture strain values ranging from 1 to 1 10 10. A fracture strain value of 1 10 5 produced fragmentation rather than extreme deformation. Values larger than 1 10 5 did not show fragmentation, while smaller fracture strain values produced the same fragmentation behavior as 1 10 5. The approach of this modeling effort was to determine what value would be best suited for the fracture strain in subsequent modeling. There was no discernable change in the element deletion at fracture strain values below 1 10 5, so the fracture strain for all subsequent 3D modeling was set to 1 10 5. The results of the axisymmetric modeling are shown in Fig. 3. The fracture strain determined using axisymmetric modeling is four orders of magnitude lower than the quasi-steady state elongation of 22% shown in Table 1 [16]. This is due to the very high strain-rates found during loading. At high strain-rates, metals change from a ductile failure to a brittle cleavage type fracture, which significantly decreases the elongation at failure [17]. This is quite evident in BCC metals. The element type used was a 4 node bilinear axisymmetric quadrilateral with reduced integration and hourglass control with the Abaqus designation CAX4R defined in Table 4. Element deletion signifying fragmentation first occurred
Fig. 6. Top down view of the inner diameter shift for eccentric models. Point a. shows the center of the outside diameter of the ring while point b. shows the shifted inner diameter. The total shift is 1 mm, and is not to scale. Wall thickness in the thinned section is approximately 6 mm.
at approximately 32 μs. Abaqus results show the peak von Mises stress to be 1.41 GPa. 3.3. Three dimensional modeling After determining the fracture strain value from the axisymmetric modeling, the 3D model was developed and evaluated. Two meshes, shown in Fig. 4, were used in modeling to check for convergence. Fig. 4 shows each mesh before calculation. The left mesh is an 8 node linear brick element with reduced integration and hourglass control with Abaqus designation C3D8R, while the right model uses 10 node quadratic tetrahedral elements with Abaqus designation C3D10M. These elements were selected since
76
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
Fig. 7. Fragmentation results for brick and tetrahedrally meshed rings with eccentricity defect. Notice that the quadratic elements in the tetrahedrally meshed model predict a higher von Mises stress. This is due to numerical noise of the surface elements in the quadratic tetrahedrally meshed model. Stress units are in Pascal.
Fig. 8. Ring refined with approximately 2 106 brick elements. It is evident that refining the mesh has increased the fragmentation resolution immensely. Stress units are in Pascal.
they are 3D stress elements that will provide accurate results under dynamic loading. Table 4 gives background information about these elements. The tetrahedral elements may need much more refinement to address any convergence issues associated with triangle-based elements. Results of the 3D calculation were similar to the axisymmetric model results. However, there appears to be a perfect rotational symmetry in the brick mesh of the 3D model that is giving misleading results that are shown in Fig. 5. This is due to how the solver determined a failure initiation site in the thick section of the brick model. Abaqus does not identify any critical flaw size for failure in this model, so the solver will assume that the part is perfect. The notch expanded first and elements were deleted, but none of the elements in the thick section were deleted as was found in the tetrahedral mesh. The tetrahedrally meshed part does not have any regularity in the mesh, i.e. the mesh does not have any rotational periodicity, so it shows failure throughout the part, not just at the notch section. Additionally, the quadratic tetrahedrally meshed model predicts a much higher maximum von Mises stress in the surface regions than the brick mesh. This is due to an over-prediction of stress formed from the design of the quadratic tetrahedral elements [18]. The nodal stresses in the quadratic elements are much higher due to the reduced averaging between gauss points in the element. Basic elements average the stresses across the entire element, while quadratic elements average between mid-nodes, leading to the change in stress approximation. In both models, a maximum strain-rate, ε_ of 1 104 was found. The peak stresses found before fracture at this
strain-rate are twice the ultimate tensile strength and three times the yield strength of the strain-rate independent material properties for AISI 4340. This correlates well with the theoretical estimate shown in Fig. 2. Fig. 5 shows the results of the calculations after 45 μs. 3.4. Eccentricity introduction After analyzing the previous results, it was thought that the approach used by Abaqus may be problematic in determining the failure origin of the ring due to the fully rotational symmetry and could lead to artificially inflated stress results. To address this issue, eccentricity was introduced into both the brick and tetrahedrally meshed Lagrangian models to simulate a defect in the metal, or a manufacturing induced eccentricity. The inner diameter of the ring was shifted 1 mm to remove concentricity which would destroy any rotational symmetry that was in the previous 3D models. Fig. 6 shows a schematic drawing of the shift. Both models were then re-meshed to add the eccentricity. Fragmentation in these models was more prevalent near the thinned section as expected. The quadratic tetrahedral mesh still predicted a larger maximum von Mises stress than the brick model. It appears that quadratic tetrahedral elements have a large amount of numerical noise at the surface elements caused by the extreme deformation. These elements over-predicted the surface stress by a factor of three. If the first three nodal stresses were disregarded, the overall bulk stress in the ring remains nearly the same as what was found in a linear tetrahedral or linear brick
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
elements. Another possible explanation is that stress concentrations exist at the notch area and are not detected by the linear elements. The results of each analysis presented in Fig. 7 show the eccentric model results that output excessively high stresses in the quadratic tetrahedral model. To avoid the skewed high surface stresses, subsequent data is reported at least three integration points away from the surface in the quadratic tetrahedral model.
Fig. 9. von Mises stress versus step time in AISI 4340 for three different element types used in eccentric uniformly loaded model. There is good agreement in data until approximately 40 μs due to the onset of element deletion. Note: the strain-rate independent yield strength of AISI 4340 is 470 MPa.
77
3.5. Mesh refinement and geometry effects Mesh refinement was done to ensure convergence of the stresses. The number of elements was increased from 6.9 103 to over 2 106 in the tetrahedrally meshed models. The increased mesh refinement predicts a similar maximum von Mises stress, but provided much better fragmentation resolution. At 60 μs, onset of a fragmentation pattern predicted by Mott et al. appeared [19,20]. Results of the mesh refinement modeling are shown in Fig. 8. The von Mises stress was plotted against step time using three different element types. A linear tetrahedral mesh C3D4, a quadratic tetrahedral mesh C3D10M, and a linear brick mesh C3D8R were used. These element types are shown in Fig. 8 and given in Table 4. The stress values were averaged across the entire model to reduce local effects due to element deletion. Stress relaxation in the elements that surround a deletion zone caused extreme fluctuations in the nodal stress outputs. To alleviate these issues, either the time stepping must be significantly decreased, or a bulk stress value must be utilized. These data showed good agreement until approximately 40 μs, where it was assumed the effects of element deletion start to significantly affect the stress output. The local stress values were heavily dependent on the distance from the location of element deletion. This may need further investigation. The numerical data are shown in Fig. 9. To ensure that the numerical noise in the quadratic tetrahedral model was due to surface effects and not model geometry issues such as stress concentrations, a shorter ring was developed and the average stress was compared against the standard concentric ring under uniform loading and linear brick elements. Brick elements were chosen since they did not show signs of artificially inflated stresses of the previous models, and demonstrate a change in the surface nodal stresses if geometry is mesh dependent. Both ends of the shorter ring are equidistant from the center of the ring notch. This is to help eliminate unequal amounts of torque applied to the notch area and should help indicate any differences in von Mises stress related to the model geometry. Fig. 10 shows the standard and new short ring geometry differences. The predicted fragmentation for both the original notched ring, left, and the new equal geometry ring, right, can be seen in Fig. 11. Both geometries predict similar fragmentation patterns. There is also no significant change in von Mises stress, which indicates that the mesh is not directly geometry dependent. 3.6. Introduction of non-uniform, time-dependent loading using CTH
Fig. 10. Equal length geometry, right, used to identify geometry related effects that may be associated with the notched ring.
A constant load was used in previous modeling steps in Abaqus to simulate the pressure experienced during explosive detonation.
Fig. 11. Comparison between the original ring, left, and the new equal geometry ring, right. It appears that the von Mises stress and the fragmentation patterns are very similar. Stress units are in Pascal.
78
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
The use of a constant load simplified the model to help eliminate any unexpected effects that may change the fracture strain value. However, this is not representative of explosive detonation. A time-dependent loading curve was introduced to the ring model to better understand the effect of non-steady state loading. To assist in developing a time-dependent loading curve, the fragmentation behavior of the symmetric ring was compared using a shock physics hydro-code developed by Sandia National
Fig. 12. Predicted deformation produced by Sandia CTH code. Notice that the thin section of the ring is not deleted as in the Abaqus model.
Laboratories named CTH [21]. The results of this modeling is shown in Fig. 12. There are few noticeable differences between the Abaqus results and the CTH results. CTH uses a plastic explosive loading curve rather than a fixed uniform pressure leading to significant variations in the model, which must be investigated. The plastic explosive, Composition B, was used in this model to simulate the explosive loading curve due to its well defined properties and wide use in commercial and military applications [22]. A tabular loading curve was developed using the time and pressure data determined by CTH. The tabular loading curve was used in subsequent Abaqus modeling. The tabular loading curve was used in place of a more sophisticated equation of state model such as the Jones–Wilkins–Lee model, so that any unnecessary contact definitions could be eliminated [23]. Eliminating contact definitions significantly reduces the model run time, as well as increases the stable time increment. It is important to note that each respective loading curve may have a significant effect on the behavior of the material. Fig. 13 plots the transient Composition B pressure to the wall used in CTH, against the constant load from previous modeling.
3.7. Time dependent loading using Abaqus
Fig. 13. Plot of the transient Composition B loading curve versus the static uniform applied load used in previous modeling. Notice the Composition B ultimate pressure occurs around 5.6 μs.
The loading curve developed using CTH was introduced into the fully symmetric ring using Abaqus, with an analysis time of 2.5 10 5 s and fixed stable time increments of 1 10 11 s. Fixing the stable time increment helps in outputting results at regular time intervals. The mesh was further refined in this model to increase convergence and eliminate any issues with exceeding the speed of sound through elements at the initial pressure face. Due to the excessive deformation of this model, it was also beneficial to have a brick element mesh that is oriented in such a manner that the brick is elongated in the direction of deformation. This allowed the elements to be squeezed into its optimal shape and then further deform into the inverse of its original shape. This method allows for re-meshing algorithms to be avoided, reducing model run time and increasing accuracy [11]. Utilizing these meshing features, a fine mesh of approximately 2 106 Abaqus C3D8R elements, mentioned in Table 4, can be used to effectively analyze this ring. The results of the refined mesh and non-steady-state loading are shown in Fig. 14. There is a noticeable change in the fragmentation behavior in the ring due to the change in loading.
Fig. 14. Results of shock loaded AISI 4340 ring after 5 10 5 s using brick elements. Notice fragmentation occurring throughout the entire model. Fragmentation pattern is noticeably different than that of the uniformly loaded model shown in Fig. 8. Stress units are in Pascal.
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
This is expected due to the significant increase in pressure that occurs in the transient Composition B loading curve. 4. Conclusions In this work, codes were used 4340 steel ring. conclusions have effort are:
axisymmetric, 3D Lagrangian, and 2D Eulerian to model fragmentation of a high strength AISI Several materials and finite element related made. The materials related conclusions for this
A fracture strain of 1 10 5 was found to initiate fragmenta
tion of AISI 4340 steel leading to element deletion around 30– 40 μs in all cases examined. Loading profile can greatly influence the fragmentation behavior as shown in the transient and steady-state cases The peak von Mises stress found in both loading cases was approximately 1.5 GPa. The Johnson–Cook model only predicts linear increases in plastic stress past the yield strength of the material, and does not account for material relaxation effects. Tensile high-strain-rate data is not readily available at these strain-rates, making this type of dynamic modeling very difficult. The finite element related conclusions for this effort are:
Meshing procedures have a significant effect on the fragmenta-
tion behavior and the von Mises stress. A linear brick mesh and linear tetrahedral mesh were found to provide the best agreement in stresses and fragmentation when reporting nodal stresses. Stress prediction during fragmentation is highly dependent on the distance from a newly deleted element. Brick elements provided the fastest convergence with the least amount of elements, decreasing calculation times. The outer three nodes of the quadratic tetrahedral elements greatly over-predicted the von Mises stress at the surface of the model compared to the linear tetrahedral elements and the linear brick elements. Stress over-prediction in quadratic elements is due to the addition of mid-nodes and reduced stress averaging across the element. There was good agreement with regards to the von Mises stress in axisymmetric and 3D Lagrangian models. Three dimensional fragmentation effects are lost in axisymmetric models. Axisymmetric modeling can be used to develop parameters for 3D models, reducing calculation times. Convergence in this modeling only was successful after mesh refinement and modifying time stepping to extremely small increments.
79
References [1] G.R. Johnson, W.H. Cook, A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures, in: Proc. 7th Int. Symp. Ballist. 547 (1983) 541–547. 10.1038/nrm3209. [2] G.R. Johnson, W.H. Cook, Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures, Eng. Fract. Mech. 21 (1985) 31–48. http://dx.doi.org/10.1016/0013-7944(85)90052-9. [3] F.J. Zerilli, R.W. Armstrong, Dislocation-mechanics-based constitutive relations for material dynamics calculations, J. Appl. Phys. 61 (1987) 1816. http://dx.doi. org/10.1063/1.338024. [4] W.S. De Rosset, Patterned Armor Performance Evaluation for Multiple Impacts. 〈http://oai.dtic.mil/oai/oai?verb=getRe cord&metadataPrefix=html&identifier=ADA417165〉, 2003 (accessed February 25, 2015). [5] B. Hopkinson, A method of measuring the pressure produced in the detonation of high explosives or by the impact of bullets, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 213 (1914) 437–456. [6] R.M. Davies, A critical study of the Hopkinson pressure bar, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci. 240 (1948) 375–457 〈http://www.jstor.org/ stable/91440〉. [7] B. a Gama, S.L. Lopatnikov, J.W. Gillespie, Hopkinson bar experimental technique: a critical review, Appl. Mech. Rev. 57 (2004) 223. http://dx.doi. org/10.1115/1.1704626. [8] S. Nemat-Nasser, Introduction to high strain rate testing, in:, ASM Handbook, vol. 8 Mech. Test. Eval., ASM Interational (2000) 427–428. [9] SolidWorks Premium Edition, (n.d.). [10] AISI 4340 Steel, Matweb Mater. Prop. Data. (n.d.). 〈http://matweb.com/search/ DataSheet.aspx?MatGUID=fd1b43a97a8a44129b32b9de0d7d6c1a〉. [11] ABAQUS Documentation, in: Version 6., Dassault Systemes, Providence, RI, USA, n.d. [12] J.J. Yoh, M.A. McClelland, J.L. Maienschein, C.M. Tarver, Test-based thermal explosion model for HMX, Proc. Combust. Inst. 31 (2007) 2353–2359. http: //dx.doi.org/10.1016/j.proci.2006.08.071. [13] J. Forbes, F. Garcia, C. Tarver, P. Urtiew, D. Greenwood, K. Vandersall, Pressure Wave Measurements During Thermal Explosion of HMX-Based High Explosives, in: Int. Detonation Symp. 〈http://www.intdetsymp.org/detsymp2002/ PaperSubmit/FinalManuscript/pdf/Forbes-168.pdf〉, 2002 (accessed January 9, 2015). [14] K.-J. Bathe, E. Ramm, E.L. Wilson, Finite element formulations for large deformation dynamic analysis, Int. J. Numer. Methods Eng. 9 (1975) 353–386. http://dx.doi.org/10.1002/nme.1620090207. [15] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, 2000 〈http://books.google.com/books/about/Nonlinear_ finite_elements_for_continua_a.html?id=C6goAQAAMAAJ&pgis=1〉, accessed January 9, 2015). [16] J.M. Schreiber, I. Smid, T.J. Eden, Fragmentation of a steel ring under explosive loading, in: M. De Graef, H.F. Poulsen, A. Lewis, J. Simmons, G. Spanos (Eds.), 1 Int. Conf. 3D Mater. Sci., John Wiley & Sons, Inc., Hoboken, NJ, USA., , , 2012: pp. 43–48. http://dx.doi.org/10.1002/9781118686768. [17] R.W. Armstrong, S.M. Walley, High strain rate properties of metals and alloys, Int. Mater. Rev. 53 (2008) 105–128. [18] E. Wang, T. Nelson, R. Rauch, C. Gmbh, Back to Elements—Tetrahedra vs. Hexahedra, 2004. [19] N.F. Mott, Fragmentation of Shell Cases, Proc. R. Soc. A Math. Phys. Eng. Sci. 189 (1947) 300–308. http://dx.doi.org/10.1098/rspa.1947.0042. [20] V.M. Gold, E.L. Baker, N.W. Koon, J.M. Hirlinger, A Method for Predicting Fragmentation Characteristics of Natural and Preformed Explosive Fragmentation Munitions, Picatinny Arsenal, 2001. [21] Sandia National Laboratory, CTH Capabilities, (n.d.). 〈http://www.sandia.gov/ CTH/Capabilities-leftnav-level1_js.html〉. [22] G. Security, Explosives-Compositions, (n.d.). 〈http://www.globalsecurity.org/ military/systems/munitions/explosives-compositions.htm〉. [23] G. Baudin, R. Serradeill, Review of Jones–Wilkins–Lee equation of state, EPJ Web Conf. 10 (2011) 00021. http://dx.doi.org/10.1051/epjconf/20101000021.
Contents lists available at ScienceDirect
Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel
Prediction of fragmentation and experimentally inaccessible material properties of steel using finite element analysis Jeremy M. Schreiber a,b,n, Ivi Smid a, Timothy J. Eden b, David Jann c a
Department of Engineering Science and Mechanics, The Pennsylvania State University, University Park, PA 16802, USA Applied Research Laboratory, The Pennsylvania State University, University Park, PA 16802, USA c DE Technologies, King of Prussia, PA 19406, USA b
art ic l e i nf o
a b s t r a c t
Article history: Received 25 February 2015 Received in revised form 21 May 2015 Accepted 2 June 2015 Available online 26 June 2015
High strain-rate properties of materials are needed for predicting material behavior in extreme environments. The demand for high strain-rate properties continues to increase for commercial and military applications as the operating environments become more extreme, such as fragmentation, impact and explosions. To reduce time and expense, Finite Element Analysis (FEA) is being used to simulate these behaviors and reduce the number of experiments needed to characterize how a material performs at high-strain-rates. A finite element model for predicting fragmentation behavior of a high strength steel ring was developed using Abaqus Computer Aided Engineering (Abaqus) software. AISI 4340 steel, a low alloy Cr–Ni–Mo steel, was used in the analysis. The results of the finite element model were compared to the results from CTH, a two-dimensional Eulerian shock physics hydro-code. CTH was also used to develop a transient loading curve for the Abaqus model. The fracture strain in the model was adjusted to induce failure in the ring. Element deletion was used to model failure. A fracture strain less than 1 10 5 was needed to initiate fragmentation. The effects of mesh type and model defects were also investigated. & 2015 Elsevier B.V. All rights reserved.
Keywords: Fragmentation Finite element analysis Johnson–Cook High strain-rate Lagrangian analysis Eulerian analysis
1. Introduction Material fragmentation at high strain-rates is difficult to characterize experimentally. Typically, a test article or projectile is subjected to a high-strain-rate that results in fragmentation. The fragments are collected and sorted into bins. Several fragments are examined under a microscope to determine failure characteristics, making this process very time consuming and expensive. In the past 60 years, many empirical material models have been developed to predict failure at high strain-rates to correlate with experimental results [1–4]. Most of these models rely on high strain-rate data obtained using the Split Hopkinson Pressure Bar (SHPB) or other high strain-rate testing. The SHPB method measures stress pulse propagation through a metal bar to predict the stress–strain relationships of a material [5,6]. There are shortcomings to this method as it requires many assumptions to be made, and loses accuracy when materials undergo tension testing at strain-rates
n Corresponding author at: The Pennsylvania State University, 3075 Research Drive, State College, PA 16801, USA. Tel.: þ1 814 865 1096. E-mail address: [email protected] (J.M. Schreiber).
http://dx.doi.org/10.1016/j.finel.2015.06.001 0168-874X/& 2015 Elsevier B.V. All rights reserved.
above 103 s 1 [7,8]. Components can be designed to improve fragmentation behavior and identify failure initiation sites. The models help identify failure modes, and reduce the number of samples that need to be tested. For this effort, a model was developed that incorporates high strain-rate data from previous SHPB experiments in a constitutive material model to predict failure.
2. Experimental procedure and approach 2.1. Ring development A three dimensional solid model of a ring was developed using SolidWorks computer aided drafting (CAD) software [9]. The ring has a diameter of 81 mm, a wall thickness of 7 mm, and a notch depth of 3.5 mm. The notch is 5.7 mm tall with a 60 degree taper at the top. The ring was imported into Abaqus for analysis. The ring was partitioned into smaller sections to allow for the use of the automatic meshing function in Abaqus which is needed to utilize all mesh types when modeling high-strain-rate behavior. Partitioning in Abaqus was done by the three point method to develop a partitioning plane. The SolidWorks model of the ring is shown in Fig. 1.
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
73
Fig. 1. Three Dimensional SolidWorks model of the ring. The cross section view is shown on the right. The diameter of the ring is 81 mm. Wall thickness is 7 mm.
2.2. Material properties
Table 1 Strain-rate independent material properties of AISI 4340 steel [1,2,10].
The material properties used in the model of the ring were those of austenitized, quenched, and tempered AISI 4340 steel determined by Johnson and Cook [1,2]. The hardness of this steel was reported to be Rockwell C30, which would correlate to a mainly pearlitic microstructure. The standard, strain-rate independent material properties for this alloy are found in Table 1 [1,2,10]. High strain-rate properties were added to the model using the Johnson–Cook (J–C) plasticity model [1,2] This empirical model contains five experimentally determined constants and is commonly used to predict behavior of ductile materials under high strain-rate conditions such as explosive loading. This model is available in the Abaqus software by default [1,2,11]. It defines the von Mises flow stress of a material as a function of the power law relationship of plastic strain and strain hardening, strain-rate, and thermal softening. The J–C model is shown in Eq. (1) [1,2]. h ih i σ ¼ A þBεnp 1 þ Clnεnp ½1 T m ð1Þ H
Density (kg/m3)
200
7830
Poisson’s % ratio Elongation 0.29
Ultimate tensile strength (MPa)
Yield strength (MPa)
745
470
22
Table 2 Johnson–Cook material constants for AISI 4340 steel [1,2]. A (MPa)
B (MPa)
C
n
m
Tmelt (K)
800
510
0.014
0.26
1.03
1793
Predicted von Mises flow stress vs. true strain
1500
εp ¼ normalized effective plastic strain rate n
ðT T room Þ ðT melt T room Þ
The values of A, B, C, n, and m are experimental constants that are determined using uniaxial tension tests and the SHPB tests. Tmelt and Troom are the melting temperature of the alloy and the ambient air temperature, respectively. The material constants used in the J–C model and are found in Table 2 [1,2]. The J–C equation was used to predict the room temperature von Mises flow stress at six different strain-rates shown in Fig. 2. The results indicate that the von Mises flow stress increases as the strain-rate increases. Strain is only plotted to 0.4 strain due to the fact that the Johnson–Cook model only predicts a linear increase in strength at large amounts of strain. This linear increase at large strains is due to the effect of the exponential power law factor in the equation. The J–C strength model was not the only constitutive material model needed in this analysis. The J–C fracture model was used to initiate the failure of the elements through the use of a fracture strain value. Fracture strain can be measured easily under uniaxial tension at quasi-static strain rates using an extensometer, but at high-strain-rates this may not be possible. Modeling a complex structure at various fracture strain values maybe the only way to identify the proper fracture strain. In modeling, once an element reaches the fracture strain value, it is deleted from the equation. The J–C fracture model does not follow a typical nucleation and
von Mises Flow Stress (MPa)
1400
εp ¼ equivalent plastic strain
T H ¼ homologous temperature ¼
Elastic modulus (GPa)
1300 1200 1100 1000 900 800 Selected Strain Rates
700 600
0
0.05
0.1
0.15
0.2
0.00001
0.0001
1
100
0.25
0.3
0.01 10000
0.35
0.4
True Strain (mm/mm)
Fig. 2. von Mises Flow Stress vs. True Strain of AISI 4340 Steel at Room Temperature. Data is plotted at Room Temperature (300 K).
growth model. It is dependent on strain, strain-rate, temperature, and pressure [1,2]. The benefits are that it is less complicated to use, and most of the parameters can be found using adjusted quasi-steady state data. The J–C fracture model is also built directly into Abaqus, and is fully compatible with the J–C strength model. The J–C fracture model is defined in Eq. (2) [2], h i n εf ¼ D1 þ D2 expD3 nσ 1 þ D4 ln ε_ n ½1 þ D5 T n ð2Þ where εf is the equivalent strain to fracture, D1 –D5 are material constants, ε_ n is the dimensionless strain-rate, and T is temperature. Material constants for AISI 4340 from Johnson-and Cook are shown in Table 3 [2].
74
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
2.3. Model setup The first step was to start with an axisymmetric model of the ring to determine how large of a displacement was needed to initiate failure based on fracture strain. The fracture strain was introduced into the finite element solver using the J–C fracture model [11]. The values of fracture strain were arbitrarily chosen from very small strains to large strains, which necessitates multiple model runs. These values were input into the damage evolution option built into the Abaqus software. This option requires a displacement at failure input. This was set to the predicted fracture strain. Axisymmetric modeling allows for fast calculation times to determine rough values for overall von Mises stress, displacement and strain-rate before implementing a more calculation intensive 3D model. A dynamic, explicit analysis with an analysis time of 60 μs was chosen for all models to cover the entire high-strainrate event. An instantaneous pressure of 1.2 GPa was applied to the inner side of the ring to simulate explosive loading conditions in the axisymmetric models and initial 3D modeling [12,13]. This value was chosen to give both a baseline value for evaluation, but also a maximum pressure load found in some high explosives. The model was constrained on top and bottom faces of the ring to prevent model runaway. Runaway is a condition where the model is under-constrained and does not perform as predicted. This boundary condition was propagated to all subsequent models. Once the baseline behavior was determined using axisymmetric modeling, a 3D model was developed. A time dependent loading curve was investigated using 3D modeling only.
2.3.1. Eulerian and Lagrangian analysis Two different methods that are commonly used to analyze dynamic material systems are the Lagrangian and Eulerian [14]. The difference between these two methods is that in a Lagrangian method, the mesh is associated with a certain section of the part, whereas in the Eulerian method, the mesh is associated with a larger computational domain and the part is not assigned to the mesh, but moves freely throughout the meshed domain. The advantage of using a Lagrangian analysis is that the mesh is fixed to certain coordinates on the part, making adding boundary conditions and contact parameters easy [15]. In contrast, in the Eulerian analyses, the material is constantly moving through the computational domain of the system, leading to greater difficulty in assigning boundary and contact parameters [15]. Eulerian models excel when high strain is expected due to the retention of the mesh geometry. Mesh distortion is a concern in Lagrangian finite element models and there are re-meshing techniques that are typically used to address this issue. The downside of having a fixed mesh in an Eulerian analysis is that the overall model domain must be much larger to prevent material from leaving the computational domain, which increases calculation time [15]. In this work, the Eulerian CTH code is used to determine the loading profile for the Lagrangian Abaqus 3D model. The deformation behavior is also compared between the two models.
3. Results and discussion 3.1. Element overview
Table 3 Fracture parameters for AISI 4340 steel given by Johnson–Cook [1,2]. D1
D2
D3
D4
D5
0.05
3.44
2.12
0.002
0.61
The element types used in finite element modeling are chosen to provide specific information about a model. Each element has different properties, which give it certain characteristics for modeling such as hourglass control. This modeling effort utilized
Table 4 Summary of element types used in each modeling step. Geometry
Element
Abaqus designation
Settings
Comments
Axisymmetric
4-Node quadrilateral, linear
CAX4R
Stress/displacement element without twist
3D
8-Node brick, linear
C3D8R
3D
10-Node tetrahedron, quadratic 4-Node tetrahedron, linear
C3D10M
Reduced integration, hourglass control Reduced integration, hourglass control Modified element, hourglass control Default
3D
C3D4
Stress/displacement element Stress/displacement element, three additional displacement variables Stress/displacement element
Fig. 3. Predicted failure initiation values for uniform loading curve at t¼ 45 μs. Degree of failure is a qualitative estimate of where the model does not show any change in failure characteristics.
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
75
Fig. 4. Meshed 3D models before analysis. The left is meshed with an 8 node linear brick element, and the right side is meshed with a 10 node tetrahedral element.
Fig. 5. Results of both brick and tetrahedral models after 45 μs. The tetrahedral model gives a larger maximum von Mises stress in the surface regions. Stress units are in Pascal.
four different types of solid continuum elements. Table 4 gives an overview of the elements used in the modeling. 3.2. Axisymmetric modeling Analyses were performed using the axisymmetric models with fracture strain values ranging from 1 to 1 10 10. A fracture strain value of 1 10 5 produced fragmentation rather than extreme deformation. Values larger than 1 10 5 did not show fragmentation, while smaller fracture strain values produced the same fragmentation behavior as 1 10 5. The approach of this modeling effort was to determine what value would be best suited for the fracture strain in subsequent modeling. There was no discernable change in the element deletion at fracture strain values below 1 10 5, so the fracture strain for all subsequent 3D modeling was set to 1 10 5. The results of the axisymmetric modeling are shown in Fig. 3. The fracture strain determined using axisymmetric modeling is four orders of magnitude lower than the quasi-steady state elongation of 22% shown in Table 1 [16]. This is due to the very high strain-rates found during loading. At high strain-rates, metals change from a ductile failure to a brittle cleavage type fracture, which significantly decreases the elongation at failure [17]. This is quite evident in BCC metals. The element type used was a 4 node bilinear axisymmetric quadrilateral with reduced integration and hourglass control with the Abaqus designation CAX4R defined in Table 4. Element deletion signifying fragmentation first occurred
Fig. 6. Top down view of the inner diameter shift for eccentric models. Point a. shows the center of the outside diameter of the ring while point b. shows the shifted inner diameter. The total shift is 1 mm, and is not to scale. Wall thickness in the thinned section is approximately 6 mm.
at approximately 32 μs. Abaqus results show the peak von Mises stress to be 1.41 GPa. 3.3. Three dimensional modeling After determining the fracture strain value from the axisymmetric modeling, the 3D model was developed and evaluated. Two meshes, shown in Fig. 4, were used in modeling to check for convergence. Fig. 4 shows each mesh before calculation. The left mesh is an 8 node linear brick element with reduced integration and hourglass control with Abaqus designation C3D8R, while the right model uses 10 node quadratic tetrahedral elements with Abaqus designation C3D10M. These elements were selected since
76
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
Fig. 7. Fragmentation results for brick and tetrahedrally meshed rings with eccentricity defect. Notice that the quadratic elements in the tetrahedrally meshed model predict a higher von Mises stress. This is due to numerical noise of the surface elements in the quadratic tetrahedrally meshed model. Stress units are in Pascal.
Fig. 8. Ring refined with approximately 2 106 brick elements. It is evident that refining the mesh has increased the fragmentation resolution immensely. Stress units are in Pascal.
they are 3D stress elements that will provide accurate results under dynamic loading. Table 4 gives background information about these elements. The tetrahedral elements may need much more refinement to address any convergence issues associated with triangle-based elements. Results of the 3D calculation were similar to the axisymmetric model results. However, there appears to be a perfect rotational symmetry in the brick mesh of the 3D model that is giving misleading results that are shown in Fig. 5. This is due to how the solver determined a failure initiation site in the thick section of the brick model. Abaqus does not identify any critical flaw size for failure in this model, so the solver will assume that the part is perfect. The notch expanded first and elements were deleted, but none of the elements in the thick section were deleted as was found in the tetrahedral mesh. The tetrahedrally meshed part does not have any regularity in the mesh, i.e. the mesh does not have any rotational periodicity, so it shows failure throughout the part, not just at the notch section. Additionally, the quadratic tetrahedrally meshed model predicts a much higher maximum von Mises stress in the surface regions than the brick mesh. This is due to an over-prediction of stress formed from the design of the quadratic tetrahedral elements [18]. The nodal stresses in the quadratic elements are much higher due to the reduced averaging between gauss points in the element. Basic elements average the stresses across the entire element, while quadratic elements average between mid-nodes, leading to the change in stress approximation. In both models, a maximum strain-rate, ε_ of 1 104 was found. The peak stresses found before fracture at this
strain-rate are twice the ultimate tensile strength and three times the yield strength of the strain-rate independent material properties for AISI 4340. This correlates well with the theoretical estimate shown in Fig. 2. Fig. 5 shows the results of the calculations after 45 μs. 3.4. Eccentricity introduction After analyzing the previous results, it was thought that the approach used by Abaqus may be problematic in determining the failure origin of the ring due to the fully rotational symmetry and could lead to artificially inflated stress results. To address this issue, eccentricity was introduced into both the brick and tetrahedrally meshed Lagrangian models to simulate a defect in the metal, or a manufacturing induced eccentricity. The inner diameter of the ring was shifted 1 mm to remove concentricity which would destroy any rotational symmetry that was in the previous 3D models. Fig. 6 shows a schematic drawing of the shift. Both models were then re-meshed to add the eccentricity. Fragmentation in these models was more prevalent near the thinned section as expected. The quadratic tetrahedral mesh still predicted a larger maximum von Mises stress than the brick model. It appears that quadratic tetrahedral elements have a large amount of numerical noise at the surface elements caused by the extreme deformation. These elements over-predicted the surface stress by a factor of three. If the first three nodal stresses were disregarded, the overall bulk stress in the ring remains nearly the same as what was found in a linear tetrahedral or linear brick
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
elements. Another possible explanation is that stress concentrations exist at the notch area and are not detected by the linear elements. The results of each analysis presented in Fig. 7 show the eccentric model results that output excessively high stresses in the quadratic tetrahedral model. To avoid the skewed high surface stresses, subsequent data is reported at least three integration points away from the surface in the quadratic tetrahedral model.
Fig. 9. von Mises stress versus step time in AISI 4340 for three different element types used in eccentric uniformly loaded model. There is good agreement in data until approximately 40 μs due to the onset of element deletion. Note: the strain-rate independent yield strength of AISI 4340 is 470 MPa.
77
3.5. Mesh refinement and geometry effects Mesh refinement was done to ensure convergence of the stresses. The number of elements was increased from 6.9 103 to over 2 106 in the tetrahedrally meshed models. The increased mesh refinement predicts a similar maximum von Mises stress, but provided much better fragmentation resolution. At 60 μs, onset of a fragmentation pattern predicted by Mott et al. appeared [19,20]. Results of the mesh refinement modeling are shown in Fig. 8. The von Mises stress was plotted against step time using three different element types. A linear tetrahedral mesh C3D4, a quadratic tetrahedral mesh C3D10M, and a linear brick mesh C3D8R were used. These element types are shown in Fig. 8 and given in Table 4. The stress values were averaged across the entire model to reduce local effects due to element deletion. Stress relaxation in the elements that surround a deletion zone caused extreme fluctuations in the nodal stress outputs. To alleviate these issues, either the time stepping must be significantly decreased, or a bulk stress value must be utilized. These data showed good agreement until approximately 40 μs, where it was assumed the effects of element deletion start to significantly affect the stress output. The local stress values were heavily dependent on the distance from the location of element deletion. This may need further investigation. The numerical data are shown in Fig. 9. To ensure that the numerical noise in the quadratic tetrahedral model was due to surface effects and not model geometry issues such as stress concentrations, a shorter ring was developed and the average stress was compared against the standard concentric ring under uniform loading and linear brick elements. Brick elements were chosen since they did not show signs of artificially inflated stresses of the previous models, and demonstrate a change in the surface nodal stresses if geometry is mesh dependent. Both ends of the shorter ring are equidistant from the center of the ring notch. This is to help eliminate unequal amounts of torque applied to the notch area and should help indicate any differences in von Mises stress related to the model geometry. Fig. 10 shows the standard and new short ring geometry differences. The predicted fragmentation for both the original notched ring, left, and the new equal geometry ring, right, can be seen in Fig. 11. Both geometries predict similar fragmentation patterns. There is also no significant change in von Mises stress, which indicates that the mesh is not directly geometry dependent. 3.6. Introduction of non-uniform, time-dependent loading using CTH
Fig. 10. Equal length geometry, right, used to identify geometry related effects that may be associated with the notched ring.
A constant load was used in previous modeling steps in Abaqus to simulate the pressure experienced during explosive detonation.
Fig. 11. Comparison between the original ring, left, and the new equal geometry ring, right. It appears that the von Mises stress and the fragmentation patterns are very similar. Stress units are in Pascal.
78
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
The use of a constant load simplified the model to help eliminate any unexpected effects that may change the fracture strain value. However, this is not representative of explosive detonation. A time-dependent loading curve was introduced to the ring model to better understand the effect of non-steady state loading. To assist in developing a time-dependent loading curve, the fragmentation behavior of the symmetric ring was compared using a shock physics hydro-code developed by Sandia National
Fig. 12. Predicted deformation produced by Sandia CTH code. Notice that the thin section of the ring is not deleted as in the Abaqus model.
Laboratories named CTH [21]. The results of this modeling is shown in Fig. 12. There are few noticeable differences between the Abaqus results and the CTH results. CTH uses a plastic explosive loading curve rather than a fixed uniform pressure leading to significant variations in the model, which must be investigated. The plastic explosive, Composition B, was used in this model to simulate the explosive loading curve due to its well defined properties and wide use in commercial and military applications [22]. A tabular loading curve was developed using the time and pressure data determined by CTH. The tabular loading curve was used in subsequent Abaqus modeling. The tabular loading curve was used in place of a more sophisticated equation of state model such as the Jones–Wilkins–Lee model, so that any unnecessary contact definitions could be eliminated [23]. Eliminating contact definitions significantly reduces the model run time, as well as increases the stable time increment. It is important to note that each respective loading curve may have a significant effect on the behavior of the material. Fig. 13 plots the transient Composition B pressure to the wall used in CTH, against the constant load from previous modeling.
3.7. Time dependent loading using Abaqus
Fig. 13. Plot of the transient Composition B loading curve versus the static uniform applied load used in previous modeling. Notice the Composition B ultimate pressure occurs around 5.6 μs.
The loading curve developed using CTH was introduced into the fully symmetric ring using Abaqus, with an analysis time of 2.5 10 5 s and fixed stable time increments of 1 10 11 s. Fixing the stable time increment helps in outputting results at regular time intervals. The mesh was further refined in this model to increase convergence and eliminate any issues with exceeding the speed of sound through elements at the initial pressure face. Due to the excessive deformation of this model, it was also beneficial to have a brick element mesh that is oriented in such a manner that the brick is elongated in the direction of deformation. This allowed the elements to be squeezed into its optimal shape and then further deform into the inverse of its original shape. This method allows for re-meshing algorithms to be avoided, reducing model run time and increasing accuracy [11]. Utilizing these meshing features, a fine mesh of approximately 2 106 Abaqus C3D8R elements, mentioned in Table 4, can be used to effectively analyze this ring. The results of the refined mesh and non-steady-state loading are shown in Fig. 14. There is a noticeable change in the fragmentation behavior in the ring due to the change in loading.
Fig. 14. Results of shock loaded AISI 4340 ring after 5 10 5 s using brick elements. Notice fragmentation occurring throughout the entire model. Fragmentation pattern is noticeably different than that of the uniformly loaded model shown in Fig. 8. Stress units are in Pascal.
J.M. Schreiber et al. / Finite Elements in Analysis and Design 104 (2015) 72–79
This is expected due to the significant increase in pressure that occurs in the transient Composition B loading curve. 4. Conclusions In this work, codes were used 4340 steel ring. conclusions have effort are:
axisymmetric, 3D Lagrangian, and 2D Eulerian to model fragmentation of a high strength AISI Several materials and finite element related made. The materials related conclusions for this
A fracture strain of 1 10 5 was found to initiate fragmenta
tion of AISI 4340 steel leading to element deletion around 30– 40 μs in all cases examined. Loading profile can greatly influence the fragmentation behavior as shown in the transient and steady-state cases The peak von Mises stress found in both loading cases was approximately 1.5 GPa. The Johnson–Cook model only predicts linear increases in plastic stress past the yield strength of the material, and does not account for material relaxation effects. Tensile high-strain-rate data is not readily available at these strain-rates, making this type of dynamic modeling very difficult. The finite element related conclusions for this effort are:
Meshing procedures have a significant effect on the fragmenta-
tion behavior and the von Mises stress. A linear brick mesh and linear tetrahedral mesh were found to provide the best agreement in stresses and fragmentation when reporting nodal stresses. Stress prediction during fragmentation is highly dependent on the distance from a newly deleted element. Brick elements provided the fastest convergence with the least amount of elements, decreasing calculation times. The outer three nodes of the quadratic tetrahedral elements greatly over-predicted the von Mises stress at the surface of the model compared to the linear tetrahedral elements and the linear brick elements. Stress over-prediction in quadratic elements is due to the addition of mid-nodes and reduced stress averaging across the element. There was good agreement with regards to the von Mises stress in axisymmetric and 3D Lagrangian models. Three dimensional fragmentation effects are lost in axisymmetric models. Axisymmetric modeling can be used to develop parameters for 3D models, reducing calculation times. Convergence in this modeling only was successful after mesh refinement and modifying time stepping to extremely small increments.
79
References [1] G.R. Johnson, W.H. Cook, A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures, in: Proc. 7th Int. Symp. Ballist. 547 (1983) 541–547. 10.1038/nrm3209. [2] G.R. Johnson, W.H. Cook, Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures, Eng. Fract. Mech. 21 (1985) 31–48. http://dx.doi.org/10.1016/0013-7944(85)90052-9. [3] F.J. Zerilli, R.W. Armstrong, Dislocation-mechanics-based constitutive relations for material dynamics calculations, J. Appl. Phys. 61 (1987) 1816. http://dx.doi. org/10.1063/1.338024. [4] W.S. De Rosset, Patterned Armor Performance Evaluation for Multiple Impacts. 〈http://oai.dtic.mil/oai/oai?verb=getRe cord&metadataPrefix=html&identifier=ADA417165〉, 2003 (accessed February 25, 2015). [5] B. Hopkinson, A method of measuring the pressure produced in the detonation of high explosives or by the impact of bullets, Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 213 (1914) 437–456. [6] R.M. Davies, A critical study of the Hopkinson pressure bar, Philos. Trans. R. Soc. London. Ser. A Math. Phys. Sci. 240 (1948) 375–457 〈http://www.jstor.org/ stable/91440〉. [7] B. a Gama, S.L. Lopatnikov, J.W. Gillespie, Hopkinson bar experimental technique: a critical review, Appl. Mech. Rev. 57 (2004) 223. http://dx.doi. org/10.1115/1.1704626. [8] S. Nemat-Nasser, Introduction to high strain rate testing, in:, ASM Handbook, vol. 8 Mech. Test. Eval., ASM Interational (2000) 427–428. [9] SolidWorks Premium Edition, (n.d.). [10] AISI 4340 Steel, Matweb Mater. Prop. Data. (n.d.). 〈http://matweb.com/search/ DataSheet.aspx?MatGUID=fd1b43a97a8a44129b32b9de0d7d6c1a〉. [11] ABAQUS Documentation, in: Version 6., Dassault Systemes, Providence, RI, USA, n.d. [12] J.J. Yoh, M.A. McClelland, J.L. Maienschein, C.M. Tarver, Test-based thermal explosion model for HMX, Proc. Combust. Inst. 31 (2007) 2353–2359. http: //dx.doi.org/10.1016/j.proci.2006.08.071. [13] J. Forbes, F. Garcia, C. Tarver, P. Urtiew, D. Greenwood, K. Vandersall, Pressure Wave Measurements During Thermal Explosion of HMX-Based High Explosives, in: Int. Detonation Symp. 〈http://www.intdetsymp.org/detsymp2002/ PaperSubmit/FinalManuscript/pdf/Forbes-168.pdf〉, 2002 (accessed January 9, 2015). [14] K.-J. Bathe, E. Ramm, E.L. Wilson, Finite element formulations for large deformation dynamic analysis, Int. J. Numer. Methods Eng. 9 (1975) 353–386. http://dx.doi.org/10.1002/nme.1620090207. [15] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, 2000 〈http://books.google.com/books/about/Nonlinear_ finite_elements_for_continua_a.html?id=C6goAQAAMAAJ&pgis=1〉, accessed January 9, 2015). [16] J.M. Schreiber, I. Smid, T.J. Eden, Fragmentation of a steel ring under explosive loading, in: M. De Graef, H.F. Poulsen, A. Lewis, J. Simmons, G. Spanos (Eds.), 1 Int. Conf. 3D Mater. Sci., John Wiley & Sons, Inc., Hoboken, NJ, USA., , , 2012: pp. 43–48. http://dx.doi.org/10.1002/9781118686768. [17] R.W. Armstrong, S.M. Walley, High strain rate properties of metals and alloys, Int. Mater. Rev. 53 (2008) 105–128. [18] E. Wang, T. Nelson, R. Rauch, C. Gmbh, Back to Elements—Tetrahedra vs. Hexahedra, 2004. [19] N.F. Mott, Fragmentation of Shell Cases, Proc. R. Soc. A Math. Phys. Eng. Sci. 189 (1947) 300–308. http://dx.doi.org/10.1098/rspa.1947.0042. [20] V.M. Gold, E.L. Baker, N.W. Koon, J.M. Hirlinger, A Method for Predicting Fragmentation Characteristics of Natural and Preformed Explosive Fragmentation Munitions, Picatinny Arsenal, 2001. [21] Sandia National Laboratory, CTH Capabilities, (n.d.). 〈http://www.sandia.gov/ CTH/Capabilities-leftnav-level1_js.html〉. [22] G. Security, Explosives-Compositions, (n.d.). 〈http://www.globalsecurity.org/ military/systems/munitions/explosives-compositions.htm〉. [23] G. Baudin, R. Serradeill, Review of Jones–Wilkins–Lee equation of state, EPJ Web Conf. 10 (2011) 00021. http://dx.doi.org/10.1051/epjconf/20101000021.