Investigation on Mie-Grüneisen type shock Hugoniot equation of state for concrete

This paper ascertains that the bilinear shock Hugoniot equation of state (EOS) can model the plasticizing process of the porous media like concrete material for high-velocity impact problems successfully. The negative slope of the bilinear Hugoniot for low particle velocity regime can simulate the process that the porosity of concrete may be compressed to form shock wave in concrete, through a series of numerical analyses over the investigation on the physical phenomena. The results of particle velocity for the concrete material are also discussed to be compared with those of nonporous aluminum alloy for 100 and 1000 m/s impact velocities. All the numerical simulations were carried out by applying the bilinear shock Hugoniot EOS to concrete which was linked to the binary object of a hydrocode: ANSYS Autodyn® [1−3] through a user’s subroutine.


INTRODUCTION
Through enormous numbers of experimental shock compression tests at Los Alamos National Laboratory, Lawrence Livermore National Laboratory, Russia etc. [4−7]. (written using their present organization and country names) where were conducted after World War II, it was clarified that the relationships between the particle velocity (up) and the shock velocity (Us) for most condensed matters sustain considerably good linearity. However, it is also known that, in the case of the material having significant porosity like geological materials, it shows unusual Hugoniot characteristics in the region of relatively lower particle velocity below 1000 m/s. Mainly for rocks, the group of Thomas J. Ahrens at Caltech conducted a series of studies to point out that such unusual phenomena are caused by phase transition [8]. For the concrete, there exists a compiled and plotted graph of Hugoniot properties by Werner Riedel et al. at Ernst-Mach-Institut (EMI) in Germany that are obtained through the experimental tests at Sandia National Laboratories, in Germany, in UK and in Japan; it is very useful. Figure 1 indicates nothing less than it [9 and 10−20]. In the same paper, they also reported the material properties of the concrete components containing some different aggregates by performing material characterization tests, furthermore discussed about the mesoscopic numerical analyses by using those obtained experimental data. The concretes of the analysis objet are consist of three types as follows: 1) Conventional compressive strength (Conv. B16): 35 MPa, density: 2.314 g/cm 3 2) Conventional compressive strength (Conv. B16): 35 MPa, density: 2.48 g/cm 3 3) High strength: 135 MPa, density: 2.477 g/cm 3 They estimated the relationships of the axial direction stress vs degree of compaction (1−V/V0) for these three concretes, where V is relative volume and suffix '0' means its reference value.
In our present investigation, for three types of concretes above mentioned, the relationships of the axial direction stress vs degree of compaction were evaluated assuming homogeneous continuum media by using simplified shock Hugoniot equations of state and constitutive models. Furthermore, in order to comprehend clearly the process in which the porosities of concretes are squashed up to form shock waves by exerting impulsive loading, a similar analysis was carried out to compare with their completely different phenomena for non-porous aluminum alloy (Aℓ 2024-T351), too.  Figure 2 shows the simplified shock Hugoniot relations (US = c0 + s.up) of three types of concretes and Aℓ 2024-T351 to be assumed in the present numerical analyses. In the Fig.2 three solid lines represent the bilinear Hugoniot relations simplified from the plotted points for the 1), 2) and 3) concretes mentioned above, a dashed straight line represents the usual Hugoniot relation for the Aℓ 2024-T351 from the LASL shock Hugoniot data [21] carried out for the purpose of reference. Note that the three bilinear lines for the concretes have negative gradients in the region of the particle velocity lower than around 100 m/s. The material properties applied to the three types of concretes in the numerical analyses are listed in Table 1 and 2. Since the impact pressure of GPa order is generated in the region of the particle velocity higher than 100 m/s, it is several times higher even than the compressive strength of high strength concrete, the constitutive model of the concrete subjected to such very high pressure is thought not to be so sensitive to the result, although its EOS is of great importance, so constitutive model of concrete was assumed to go along with the elastic-perfectly-plastic model for convenience sake.  The value of the shear modulus (G, rigidity) shown in Table 1 is calculated by the following eqns. (1) and (2):

Material Models
where K0 is a static or elastic bulk modulus, (i.e. when the particle velocity is 0) ρ0 is an initial mass density, ν is a Poisson ratio and c0 is a bulk (elastic) sound velocity.
In addition to an EOS and a constitutive equation, a failure criterion was applied to the concretes. Each concrete is assumed to be broken by the loaded high pressure and lose its yield stress and all the stress deviation components, when more than 3 % equivalent plastic strain is evaluated in a numerical discretization element.
In addition, a bilinear EOS can be utilized in the standard version of ANSYS Autodyn adopted in the present numerical simulations, but any negative gradients are not allowed to be input, so preliminarily to the main calculation, we wrote a user's subroutine to enable us to input the bilinear shock Hugoniot EOS having a negative gradient. Figure 3 depicts the sketch of the present numerical analysis model. Assumed is a beam-like concrete structure having 200 mm length in x-direction and a square cross section of 40 mm×40 mm in y-and z-directions. Applied is the impulsive boundary condition at the left end (x = 0) of the beam that sustains a constant velocity at the range of 100 to 2000 m/s. The opposite end of the beam is completely free. The surfaces of y-and z-directions are assumed to be subjected to free sliding and rigid wall boundaries. The 1/4 analysis system was adopted in consideration of its symmetrical property, i.e. mirror symmetry boundaries were applied to the planes of y = 0 and z = 0.

Numerical Analysis Model
In the meanwhile, numerical meshes seem to be unnecessary in y-and z-directions, differing from mesoscopic model by W. Riedel [9]. However, since the process of shock wave propagation is taken into account in y-direction and in z-direction, as well as in x-direction, the discretization method using homogeneous square meshes, 1 mm on a side, was adopted to the whole computational domain.  Figure 4 shows the comparison of the particle velocity profiles among four selected velocity boundary conditions at 20 µs. Pay attention to the fact that as the velocity of the boundary condition becomes higher, so the rear end of the profile stretches to more righthand side, since the horizontal axis of Fig.4 is the distance from the left-hand end of the initial numerical model.
Commonly for each case, the initial static velocity of longitudinal elastic wave (cl) is estimated by the following eqn. (3): Eventually, in the case of the concrete with conventional compressive strength and low density, i.e. the compressive yield strength is 35 MPa and the density is 2.314 g/cm 3 , the value of longitudinal elastic wave is estimated to be 6010 m/s, assuming that the particle velocity is 0 m/s. The right-hand tip of each profile reaches the position of over 120 mm (=6010 m/s×20 µs), because each longitudinal elastic wave propagates to such a position by the velocity of longitudinal elastic wave derived from eqn. (3), and this front of the wave is referred to as 'Elastic Precursor'. The particle velocity of left-hand follow-on profile becomes higher after the precursor because each concrete is applied to the overburdened velocity boundary condition. Accordingly, each concrete varied from elastic to plastic state without deviation components of stress and strain, this phenomena denotes that shock wave has been generated in each concrete media finally. We can see that it takes so much time to complete this plasticizing process, as the velocity of the boundary condition becomes lower.
In the case to which the 2000 m/s velocity boundary is applied, it is recognized that after some oscillation a constant particle velocity is attained as well as the rising time is very precipitous, in comparison with other cases.
Although similar results were observed in the other concretes, it is interesting to note that no shock wave was generated only in the case of the 100 m/s velocity boundary for the conventional strength and low density concrete. This fact suggests that there are critical points around these analytical conditions. Figure 5 indicates the graph in which the profiles of particle velocity and axial-(x-) direction stress (σxx) are plotted, respectively using left and right vertical axes, in the case of the 500 m/s velocity boundary for the conventional strength and low-density concrete at 20 µs. However, for this analysis condition, since the static hydrodynamic pressure is dominant in compared to the deviation components of stress, practically the σxx can be assumed to be equivalent to the pressure. The color-coded row above the graph depicts that each corresponding numerical mesh below sustains what material status. It is fairly difficult to judge from this legend, but in the boundary between the elastic (green) and the plastic (light blue) regions some disorder is generated in the y-direction, while no significant one in the boundary between the plastic and the failed regions. These facts seem to be caused by slight three-dimensional effect.
Observing these two profiles, the front of the elastic wave propagates to the distance of about 180 mm from the input velocity boundary, but subsequent waves sustain elastic state for a while, after that they reach yield point stress and mature plasticizing process. They finally generate material failure in the concrete and maintain constant particle velocity and axial stress after brief oscillation. The oscillation amplitude of axial stress is greater than that of particle velocity. The oscillation of particle velocity might be restrained by the inertia effect.  Figure 6 shows two kinds of graphs of maximum σxx and maximum effective plastic strain dependency on degree of compaction respectively (1−V/V0, where V is reference density, V=1/ρ, and ρ is density), obtained and plotted by totally 21 cases of numerical simulations for three sorts of concretes and 7 different constant velocity boundary conditions.

Stress and Effective Plastic Strain Dependency on Degree of Compaction
It should be noted that as the constant value of the velocity boundary becomes higher, so the degree of compaction becomes larger. Although a simple numerical model was adopted in these simulations, they successfully reproduce the W. Riedel et al.'s maximum xx graph obtained by using much more complicated mesoscopic model. On the other hand, the graph of effective plastic strain indicates very complicated variation. Especially, the case of conventional strength and high-density concrete seems to show fairly different variation from the other concretes. In the cases of conventional strength/low density and high strength concretes, as the value of velocity boundary becomes higher up to 1000 or 1500 m/s, so the effective plastic strain becomes larger monotonically, but in the impact region over 1500 m/s switches to decrease. It is coincident with the affirmation in this paper that constitutive model is of little importance and plastic flow is reduced comparatively in such a hypervelocity regime as concrete material, because deformation proceeds extremely rapidly. In contrast, in the case of conventional/high density concrete the variation of effective plastic strain shows multiple local maximum values, moreover the effective plastic strain still shows an upward tendency at the velocity boundary of 2000 m/s (about 0.35 degree of compaction). It appears quite strange, and it is necessary to make more future investigations.   Figure 7 indicates the particle velocity profiles for the constant velocity boundary conditions of 100 and 1000 m/s, when only the concrete material is replaced with Aℓ2024-T351 in Fig.3. The Mie-Grüneisen type linear shock Hugoniot equation of state [21]. (plotted by the dashed line in Fig.2) and Johnson-Cook constitutive model [22]. was applied to the Aℓ2024-T351.

Elastic Precursor Behavior of Aℓ2024-T351
As usually observed for other ductile materials, we can see clear elastic precursors and typical two-step rising profiles. The plasticizing regime shown in the corresponding results of Fig.4 and 5 cannot be seen in this graph, which is considered as squashing process of the porosity of concrete. Contrarily, we can see that clear shock wave fronts are generated in the solid or condensed material in the case of impact velocity of 1000 m/s.

CONCLUSIONS
As shown in Chapter 3, it was confirmed that the generating process of shock wave in concrete can be successfully simulated, subsequent to the plasticizing process in which the porosities of concrete are squashed, by applying Mie-Grüneisen type bilinear shock Hugoniot equation of state to concrete. The results were verified by comparing with those of more complicated mesoscopic numerical model by Werner Riedel et al. at Ernst-Mach-Institut [9]. It was also clarified that constitutive model is not a dominant factor for concrete, in the region of over several hundred m/s impact velocity. Except for the materials that have very high Hugoniot Elastic Limit (HEL) like fine ceramics, in the case of the brittle materials that have at most 0.1 GPa compressive strength like concrete, the plasticizing and shock wave generating process caused by the squash of the porosities might be adequately modeled using only simple material model of bilinear