A Dynamic Cavity Expansion Model for Rigid Projectile Penetration into Concrete considering the Compressibility and Nonlinear Constitutive Relations

The P-α equation of state (EOS) and a nonlinear yield criterion are utilized to characterize the dynamic constitutive behavior of concrete targets subjected to projectile normal penetration. A dynamic cavity expansion model considering the compressibility and nonlinear constitutive relations for concrete material is developed. Then, a theoretical model to calculate the depth of penetration (DOP) for rigid projectile is established. Furthermore, the proposed model is validated based on the available test data as well as the calculation results by the linear compressible EOS and linear yield criterion. This study shows that the proposed model derived using the P-α EOS and nonlinear yield criterion can effectively reflect the plastic mechanical properties of concrete and is also suitable for predicting the DOP of concrete targets. In addition, the influence law of concrete constitutive parameters such as the cohesion strength, shear strength, internal friction coefficient, and elastic limit pressure on the DOP is revealed.


Introduction
The dynamic response and destruction of concrete structures under explosive or impact loading have important theoretical and application value in the field of weapon damage and engineering protection [1]. Li et al. [2] systematically reviewed the empirical formulas, theoretical analyses, and numerical methods related to the concrete penetration problem. To evaluate the damage and failure behavior of concrete, a large number of penetration tests [3][4][5][6][7][8] have been conducted. The penetration process is driven not only by the mortar strength, aggregate size, aggregate sliding, and redistribution behavior but also by pore breakage and compaction. Therefore, the compressibility and nonlinear effect of concrete must be considered in the process of theoretical analysis.
Cavity expansion theory is an accepted analytical method for studying the penetration mechanism of projectiles into a target. Bishop et al. [9] proposed quasi-static cylindrical and spherical cavity expansion equations in a semi-infinite target.
In the current works that focus on reducing the complexity in the cavity expansion analysis, the simplified EOS is usually used such as incompressible [12], linearly compressible [12,22], and locked hydrostatic [11,13] models. Although the incompressibility hypothesis can reduce the difficulty of the analysis step, the DOP obtained by the incompressible EOS is much shallower than the corresponding test data [12]. The linearly compressible model is approximately available for the case of low striking velocities. Therefore, neither the linear model nor the locked model can accurately reflect the compression process of concrete.
After reviewing the yield criteria, it was found that the original and revised forms of the Tresca [11], Mohr-Coulomb [12,14,19,[22][23][24], and Drucker-Prager [15] models have been widely used in dynamic cavity expansion analyses for concrete-like material. However, all the above-mentioned yield criteria are in linear form. When considering the triaxial shear data in the present paper, the results show that the relationship between the pressure and shear strength is nonlinear, especially under high pressure. In summary, a model considering pore compaction, a nonlinear EOS, and nonlinear yield criterion are needed to calculate the DOP of projectiles into concrete targets.
This paper proposes a P-α EOS and nonlinear yield criterion for rigid projectiles penetrating concrete targets on the basis of a comparative analysis of existing concrete constitutive models. A spherical cavity expansion model considering the compressibility and nonlinear constitutive relations of concrete is established, and the cavity pressure is calculated. Furthermore, the calculation model for the DOP is derived, and the validity of the model is demonstrated by a comparison with the results of the widely recognized model as well as the available penetration test data [12]. Eventually, the influences of the constitutive parameters on the DOP are discussed.

Constitutive Model
In the following section, the hydrostatic behavior and strength model for concrete are treated separately and described by the P-α EOS and nonlinear yield criterion, respectively.

Equation of
State. The concrete exhibits pore compaction and permanent densification effect when subjected to highpressure loading [28]. The P-α EOS originally proposed by Herrmann [29] can be used to describe the compressibility of concrete. The P-α EOS for concrete can be expressed by Eq. (1).
where α and α 0 are the porosities of the current and initial material, respectively; ρ and ρ s are the densities of the current and fully compacted material, respectively; P is the current pressure, P e is the pressure at which pore compaction occurs, and P s is the fully compacted pressure.

Nonlinear Yield
Criterion. The mechanical properties and yield criterion for rock-like materials are always the con-cern of researchers in fields as mining science [30], geology [31], and rock engineering [32][33][34]. In the plastic range, to describe the confining pressure effect on yield strength for rock-like materials, the experimental data are smoothed and approximated according to the method first proposed by Lundborg [35], considering the similarities between concrete and rock in terms of mechanical response. For intact concrete, there is an initial tensile strength f t and a cohesion strength τ 0 . The local slope of the strength-pressure curve is called the internal friction coefficient μ; the shear strength reaches a constant τ m at a sufficiently high-pressure state. The nonlinear yield criterion in a spherically symmetric Euler coordinate system is where σ r and σ θ are the radial and hoop Cauchy stress components, which are taken as positive in compression. In summary, setting τ 0 = a 0 , 1/μ = a 1 , and 1/ðτ m − τ 0 Þ = a 2 , and then, Eq. (2a) can be expressed as σ r − σ θ = a 0 + P/ða 1 + a 2 P Þ. a 0 , a 1 , and a 2 are constants which obtained by a series of triaxial compression test data.

Development of Dynamic Spherical Cavity Expansion Model
The cavity expansion generates an elastic-crackedcompacted response (Figure 1(a)) or elastic-compacted response ( Figure 1(b)) which depends on the cavity expansion velocity V r , where r is the radial Eulerian coordinate, c and c 1 are the boundary velocities, c d is the dilatational velocity, and t is the time.
3.1. Compacted Region. The mass conservation and momentum conservation equations of the compressible medium in Euler coordinates are ∂σ r ∂r where r is the radial Eulerian coordinate and v is the particle velocity (outward is positive).

Geofluids
By defining the following dimensionless variables and introducing similar transformations, Transform Eq. (3) into the following dimensionless forms dU dξ h T ð Þf c dT dξ where f c and f t are the static uniaxial compressive and tensile strength of concrete. f ðTÞ, gðTÞ, hðTÞ, and iðTÞ can be derived with Eqs. (1) and (2).
From Eq. (2), the dimensionless relation between the pressure T and radial stress S obeys Eq. (5) is transformed by the Runge-Kutte method to evaluate it numerically.
When the Hugoniot jump condition is given, the initial values for U and T in Eq. (8) can be determined by the solutions in the cracked region (Figure 1(a)) or elastic region ( Figure 1(b)) at ξ = 1. These two types of target responses are presented separately below.

Elastic-Cracked-Compacted
Response. In the elastic region, the stress-strain relation is governed by Hooke's law.
Cavity Compacted Figure 1: Two types of target responses: (a) elastic-cracked-compacted response and (b) elastic-compacted response. 3 Geofluids where u is the outward radial displacement (positive direction).
Considering that elastic deformation is very small, ignoring the convective term in Eq. (3b) [14,22] and substituting Eq. (9) into Eq. (3b) obeys the wave equation By introducing a new dimensionless variable u = u/ct and adopting the similarity transformation, Eq. (10a) can be expressed in another form.
Integrating Eq. (11a) gives the following general solution.
where A 0 and B 0 are the unknown integral constants, which can be obtained using the following boundary conditions. u Combining Eqs. (12) and (13) yields Furthermore, the relation between the particle velocity and radial displacement is Substituting Eq. (12) into Eq. (15) and considering the strain in elastic region is small, Eq. (15) can be revised as Using Eqs. (9b), (12), (13), and (16), the dimensionless radial stress S 1 and particle velocity U 1 in the elastic region at ξ = β 1 /β are derived as Using the Hugoniot jump conditions in the cracked region at ξ = β 1 /β, the dimensionless radial stress S 2 and particle velocity U 2 can be obtained, i.e., Moreover, concrete is considered to be a linear compressible material in the cracked region when σ θ = 0, and then, the general solutions in the cracked region can be described as where C 0 and D 0 are the integral constants. Based on the boundary conditions in Eq. (18), it can be determined as follows: The above-mentioned nonlinear yield criterion is satisfied in the cracked region at ξ = β 1 /β; thus, By using Eq. (21), β can be expressed as In the cracked region at ξ = 1, the dimensionless radial stress S 3 , pressure T 3 , and particle velocity U 3 can be obtained as 3.3. Elastic-Compacted Response. The target response can be simplified to the elastic-compacted type when c > c 1 . The solution for u is consistent with Eq. (12), while the integral constants should be obtained from the following boundary conditions. u Substituting the boundary conditions into the displacement solution in Eq. (12), we have where x 0 , y 0 , and z 0 are constants, which can be written as In the elastic region at the elastic-compacted interface (ξ = 1), using Eqs. (9), (12), (24), (25) and (26), then S 3 , T 3 , and U 3 are obtained.
3.4. Hugoniot Jump Conditions. To solve the dimensionless conservation equation in Eq. (8) still requires a dimensionless pressure and particle velocity in the compacted region at the cracked-compacted interface or elastic-compacted interface; it can be obtained with the Hugoniot jump conditions. Noticing that subscript 3 represents the variable in the cracked or elastic region at ξ = 1, and subscript 4 represents the variable in the compacted region at ξ = 1. The concrete is treated as a linear compressible material in the elastic or cracked regions at ξ = 1, i.e., The P-α EOS is used for the compacted region; therefore, For the dimensionless Hugoniot jump conditions, it follows that By solving Eqs. (7) and (28)

Geofluids
To establish the penetration model, it is necessary to clarify the relationship between σ r and V r . The cavity expansion resistance of concrete is generally expressed as the following three components [15], i.e., where A, B, and C denote the dimensionless resistance constants.
Here, we can choose concrete with uniaxial compressive strength of 51 MPa to analyze the proposed dynamic cavity expansion model in detail. For the EOS, referring to the isotropic compression data of 51.2 MPa concrete [36], the fitting results using Eq. (1a) are shown in Figure 2(a). The initial and fully compacted densities of the concrete are 2350 kg/m 3 and 2767 kg/m 3 , respectively. Using the P-α EOS with α 0 = 1:177, P e = 17 MPa, and P s = 6370 MPa gives satisfying agreement with the test data. For the relation between the shear strength and pressure, the existing literature lacks the triaxial data of 51 MPa concrete, so the 48 MPa concrete triaxial data from Hanchak et al. [37] is used instead, and the fitted results by Eq. (2) are presented in Figure 2(b). Obviously, the nonlinear yield criterion presented in this paper can effectively reproduce the test data

Deep Penetration Analysis of an Ogive-Nosed Projectile.
The axial resistance acting on the projectile nose is given by the following expression.
If the tangential friction stress is neglected [38], the final axial penetration resistance F x can be expressed as  6 Geofluids where N 0 , N 1 , and N 2 are the nondimensional quantities related to the projectile nose shape, i.e., N 0 = 1, As Figure 3 shows, for ogive-nosed projectile, it follows that N 0 = 1, where ψ = s/d and s is defined in Figure 3. The final DOP H can be derived by integrating Eq. (33).
where D = ½4ACN 2 − ðBN 1 Þ 2 1/2 , and V 1 is projectile velocity at penetration depth 4a which was obtained from For Eqs. (35)- (37), m is the projectile mass, d = 2a is the shank diameter, and V s is the initial striking velocity.

Comparisons with the Penetration Test Data.
First, the proposed model is verified based on the test data of projectile penetration into 51 MPa concrete [4]. The P-α EOS and nonlinear yield criterion shown in Figure 2 are employed. Then, the predicted and experimental DOP are compared in Figure 4(a). It shows that the proposed model gives better agreement than the model in Ref. [12], especially at relatively high-impact velocities. Figure 4(b) illustrates the cavity surface radial stress versus the cavity expansion velocity of the elastic-cracked-compacted and elastic-cracked-plastic [12] models. Correspondingly, the cavity expansion resistance function parameters calibrated as A = 7:4, B = 2:5, and C = 1:1 and A = 5:8, B = 3:0, and C = 1:2, respectively. When V r /ð f c /ρ 0 Þ 1/2 < 1:8, the cavity surface radial stress σ r calculated by the elastic-cracked-compacted model is slightly larger than that of the elastic-cracked-plastic model. When V r / ð f c /ρ 0 Þ 1/2 > 1:8, the σ r in the elastic-cracked-compacted model is smaller than the elastic-cracked-plastic model, and the difference between the two models becomes greater with increasing of the cavity expansion velocity. This also explains why the elastic-cracked-plastic model cannot effectively predict the penetration depth under relatively high-impact velocities. Elastic-cracked-compacted Elastic-cracked-plastic [12] (b)
For 36 MPa concrete, using the hydrostatic pressure-shear strength test data from Ref. [3], the fitting results with nonlinear least squares are plotted in Figure 6(a). For the EOS, since the pressure-porosity data of 36 MPa concrete are absent, the concrete exhibits similar compaction characteristics when the strength is very close. Therefore, using the EOS parameters from 35 MPa [40] concrete seems reasonable. The constitutive parameters for 36 MPa concrete are summarized as follows: ρ 0 = 2:31 g/cm 3 , α 0 = 1:19, P e = 23:3 MPa, P s = 6000 MPa, a 0 = 22:89 MPa, a 1 = 0:72, and a 2 = 2:17 × 10 −3 MPa -1 . Figure 6(b) demonstrates that the proposed model provides better predictions than the classic model, and it can be inferred that the underestimation DOP of the classic model mainly comes from the use of the linear constitutive relations.

Further Discussions.
To obtain the parameters of the nonlinear constitutive in Eqs. (1) and (2), i.e., P e , P s , a 0 , a 1 , and a 2 , the pressure-porosity data from flyer-plate impact tests and the hydrostatic pressure-shear strength data from triaxial compression tests are needed. However, most DOP tests lack the necessary data required for analysis and calculation models. Under this circumstance, the parameters of a 0 , a 1 , and a 2 are taken from the method proposed by Forrestal et al. [6] Proposed model Model in Ref. [12] (c)    10 Geofluids Markovich et al. [41] as well as the value of P e suggested by Holmquist et al. [42]. The default values of P s for 35 MPa and 140 MPa concrete are both equal to 6 GPa in AUTODYN version 6.1 [43]. Here, the compacted pressure P s is also recommended to be 6 GPa for f c = 58:4 MPa and f c = 62:8 MPa concrete. Thus, the constitutive parameters of normal concrete are summarized as follows: a 0 = 2:442f c 0:4369 , P e = f c 3 , where f c is measured in MPa. Figure 7 gives the experimental data from Refs. [4,5] and the predicted DOP by Eq. (36), where the parameters of the P-α EOS and nonlinear yield criterion are given in Eq. (38). Most of the predictions show reasonable agreement with the test data except for two points, which are marked by a circle in Figure 7(b). The overestimations of the DOP for these two data points may lie in that if the striking velocity exceeds 987 m/s, the projectile mass loss will aggravate, and the mass abrasive penetration regime occurs.

The Influence of the Concrete Constitutive
Parameters on the DOP Taking 51 MPa concrete target as an example, the influences of the concrete constitutive parameters on the DOP are discussed below. Figure 8 illustrates the effects of one specific variable on the DOP while maintains other variables to be same. During the theoretical analysis, ρ 0 = 2:35 g/cm 3 and α 0 = 1:177. From Figure 8(a), it is observed that if a 0 increases by 50%, the DOP decreases by 3.8%, 3.4%, and 3.2% when the striking velocity is 400 m/s, 800 m/s, and 1200 m/s, respectively. Obviously, the DOP is relatively less affected by the cohesion in different velocity ranges.
As Figure 8(b) shows, if a 1 increases by 50%, the DOP increases by 6.8%, 7.1%, and 7.2% when the striking velocity is 400 m/s, 800 m/s, and 1200 m/s, respectively. In this respect, the coefficient of internal friction μ has a considerable influence on the DOP.
As can be seen in Figure 8(c), if a 2 increases by 50%, the DOP increases by 7%, 10.2%, and 12.5% when the striking velocity is 400 m/s, 800 m/s, and 1200 m/s, respectively. Therefore, we believe the von Mises plastic limit τ m has a substantial influence on the DOP.
As is observed in Figure 8(d), if P e increases by 50%, the DOP decreases by 2.5%, 1.9%, and 1.6% when the striking velocity is 400 m/s, 800 m/s, and 1200 m/s, respectively. In general, the influence of the elastic limit on the DOP decreases with increasing striking velocity and is negligibly small at high striking velocities.
From Figure 8(e), if P s increases by 50%, the DOP decreases by 6.8%, 12.0%, and 16.6% when the striking velocity is 400 m/s, 800 m/s, and 1200 m/s, respectively. It is observed that the influence of compacted pressure on the DOP increases with increasing striking velocities. For lower impact velocities, the effect of P s on the DOP is relatively small, and for higher striking velocities, it may become considerably large.

Conclusions
(1) The P-α EOS and nonlinear yield criterion are employed to characterize the plastic mechanical response of concrete under impact loading. Then, an improved dynamic spherical cavity expansion model considering the concrete compressibility and nonlinear constitutive relations is established (2) A penetration model to predict the DOP is obtained, and the proposed model seems to better predict the available test results compared to the classic model derived from the linear EOS and yield criterion (3) The parameters a 2 and P s have a great influence on the DOP while the internal friction coefficient μ has a considerable influence on the DOP. In contrast, a 1 and P e have a relatively small influence on the DOP and are weakly dependent on the striking velocities

Data Availability
The figures presenting the data analysis were all drawn in Origin 8.0. The formula derivation process and theoretical calculation data files used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.